{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fe944280-b8e8-48e5-88bc-ba3516c42b2a",
   "metadata": {},
   "source": [
    "# Image simulation with dask and generat\n",
    "\n",
    "The purpose here is to simulate images to identify the best methods for:\n",
    "\n",
    "- Determining the FLAT image\n",
    "- Segmenting cells from the background\n",
    "- Computing the ratio\n",
    "- Determine the minimal detectable gradient for a given error.\n",
    "\n",
    "Since subtracting the correct background value is crucial for accurate ratio imaging, we will test the distribution of background values with probplot for normality.\n",
    "\n",
    "External examples: https://cbia.fi.muni.cz/simulator/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73c15c9d-424d-4130-87b1-bc2880f74ca3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import dask.array as da\n",
    "import dask_image\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy\n",
    "import skimage\n",
    "import skimage.filters\n",
    "import skimage.io\n",
    "import tifffile as tff\n",
    "import zarr\n",
    "from dask_image import ndfilters, ndmorph\n",
    "from scipy import ndimage\n",
    "\n",
    "from nima import generat, nima, segmentation, utils\n",
    "\n",
    "nrg = np.random.default_rng(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9ee6532-ed1f-4674-81c9-918fa255d4ed",
   "metadata": {},
   "source": [
    "## FlatXY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db808977-e7e1-45c5-9fd3-3c505ff5e3a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fp_flatxy = \"/home/dati/dt-clop3/data/210920/flatxy.tf8\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da4b70f2-c47a-413a-80e9-3283c1a08bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "nima.read_tiffmd(fp_flatxy, channels=[\"G\", \"R\", \"C\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31e870b1-6e95-44b3-bd78-cebd85a33e59",
   "metadata": {},
   "outputs": [],
   "source": [
    "store = tff.imread(fp_flatxy, aszarr=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae1187b8-ce6f-41e1-94b3-be72f4fbb3ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "zc1a = zarr.open(store)\n",
    "zc1a.info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18c82e3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "dd = da.from_zarr(store)\n",
    "dd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59f2b3af",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = dd[0, 0]\n",
    "plt.imshow(img, vmax=60)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "327fccf3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from dask.diagnostics.progress import ProgressBar\n",
    "\n",
    "with ProgressBar():\n",
    "    bg_result = nima.bg(img.compute())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e2d809d-d769-4681-b889-123652dbeb99",
   "metadata": {},
   "outputs": [],
   "source": [
    "bg_result.figures[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "158cd5f6-aecd-4ea2-aecc-d010ef6ffd52",
   "metadata": {},
   "outputs": [],
   "source": [
    "im = img.compute()\n",
    "bgmax = segmentation._bgmax(im, bins=25, densityplot=True)\n",
    "bgmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df30cd47-cff7-43bd-9dfe-cf003e605826",
   "metadata": {},
   "outputs": [],
   "source": [
    "bg_result = segmentation.bg_refine_iteratively(im, bgmax=bgmax, probplot=1)\n",
    "bg, sd, iqr, figs = bg_result.bg, bg_result.sd, bg_result.iqr, bg_result.figures\n",
    "bg, sd, iqr, figs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "135cce78",
   "metadata": {},
   "outputs": [],
   "source": [
    "pp = da.mean(\n",
    "    dask_image.ndfilters.maximum_filter(dd[0:4000:20, 0], size=(100, 1, 1)), axis=0\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70ae847c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ppp = pp.compute()\n",
    "\n",
    "plt.imshow(skimage.filters.gaussian(ppp, 100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82524ece",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = img < skimage.filters.threshold_mean(img)\n",
    "skimage.filters.threshold_mean((img * m).clip(np.min(img))).compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0278c953",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img < skimage.filters.threshold_triangle(img))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a2feecc-b734-4278-9e94-a06272451dec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dabg(di):\n",
    "    m = di < skimage.filters.threshold_mean(di)\n",
    "    m1 = di < skimage.filters.threshold_mean((di * m).clip(np.min(di)))\n",
    "    # m2 = ndmorph.binary_dilation(~m1)\n",
    "    ndmorph.binary_dilation(~m1)\n",
    "    return da.ma.masked_array(di, mask=~m1)\n",
    "\n",
    "\n",
    "def bg(im):\n",
    "    m = im < skimage.filters.threshold_mean(im)\n",
    "    m1 = im < skimage.filters.threshold_mean((im * m).clip(np.min(im)))\n",
    "    m2 = skimage.morphology.binary_dilation(~m1, footprint=np.ones([1, 1]))\n",
    "    # m2 = im < skimage.filters.threshold_triangle(np.ma.masked_array(im, mask=~m))\n",
    "    return np.ma.masked_array(im, mask=m2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a721c405",
   "metadata": {},
   "outputs": [],
   "source": [
    "tff.imshow(dabg(img).compute())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4202c940-6bf8-416a-b088-38b995ca6abf",
   "metadata": {},
   "outputs": [],
   "source": [
    "flat = np.ma.mean(dabg(dd[333:500:1, 0]).compute(), axis=0)\n",
    "skimage.io.imshow(flat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12a13285",
   "metadata": {},
   "source": [
    "### threshold mean clipping to min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a010633",
   "metadata": {},
   "outputs": [],
   "source": [
    "tff.imshow(bg(dd[0, 0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "459199ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(bg(img).ravel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f833df46",
   "metadata": {},
   "outputs": [],
   "source": [
    "[np.ma.median(bg(dd[i, 0])) for i in range(10)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38c9290e",
   "metadata": {},
   "outputs": [],
   "source": [
    "segmentation.bg_refine_iteratively(img.compute(), bgmax=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa38e2c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "np.ma.mean(bg(img.compute()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c18f01f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "utils.bg2(img.compute(), bgmax=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f60d94f",
   "metadata": {},
   "source": [
    "### masked array (ma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19722a05",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.ma.masked_array([1, 4, 3], mask=[False, False, True])\n",
    "b = np.ma.masked_array([10, 2, 6], mask=[False, True, False])\n",
    "\n",
    "np.ma.median([a, b], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52f6d16c",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.ma.median(np.ma.stack([a, b]), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ba6ec27",
   "metadata": {},
   "outputs": [],
   "source": [
    "f3 = img > skimage.filters.threshold_local(img.compute(), 601)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31f0bf2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "f3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3326015",
   "metadata": {},
   "outputs": [],
   "source": [
    "img[~f3].mean().compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "680538ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "m1 = np.ma.masked_greater(img, 15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "264c51f7",
   "metadata": {},
   "source": [
    "## generat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48b08607-6e19-44b4-8a19-371cba20605c",
   "metadata": {},
   "outputs": [],
   "source": [
    "image = \"bias + noise + dark + flat * (sky + obj)\"\n",
    "image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f970ce93",
   "metadata": {},
   "source": [
    "- **bias**: generate a wave-like shape along x.\n",
    "- **noise**: random number will do.\n",
    "- **dark**: simply a scalar value.\n",
    "- **flat**: generate some 2D parabolic shape.\n",
    "- **obj**: circles-ellipsis. (MAYBE: like finite fractals to compare segmentation).\n",
    "- **sky**: None | some blurred circle-ellipsoid coincident and not with some obj.\n",
    "\n",
    "fg_prj := \n",
    "\n",
    "bg_prj := "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f79439da-ffec-4c80-bc65-8222ce9c7e67",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = Y = 128\n",
    "\n",
    "plt.figure(figsize=(12, 2.8))\n",
    "\n",
    "plt.subplot(1, 5, 1)\n",
    "plt.title(\"BIAS\")\n",
    "bias = generat.gen_bias(Y, X)\n",
    "skimage.io.imshow(bias)\n",
    "\n",
    "plt.subplot(1, 5, 2)\n",
    "plt.title(\"FLAT\")\n",
    "flat = generat.gen_flat(Y, X)\n",
    "skimage.io.imshow(flat)\n",
    "\n",
    "plt.subplot(1, 5, 3)\n",
    "plt.title(\"Object\")\n",
    "single_obj = generat.gen_object(Y, X, max_radius=7)\n",
    "skimage.io.imshow(single_obj)\n",
    "\n",
    "plt.subplot(1, 5, 4)\n",
    "plt.title(\"OBJS\")\n",
    "objects = generat.gen_objs(\n",
    "    generat.ImageObjsParams(\n",
    "        max_fluor=100, max_num_objects=13, max_radius=12, min_radius=6, ncols=Y, nrows=X\n",
    "    )\n",
    ")\n",
    "skimage.io.imshow(objects)\n",
    "\n",
    "plt.subplot(1, 5, 5)\n",
    "plt.title(\"Frame\")\n",
    "frame = generat.gen_frame(objects, bias=bias, flat=flat, sky=17, noise_sd=3)\n",
    "skimage.io.imshow(frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68049b10-909a-4333-b2d1-75c03b4b38e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "objects = generat.gen_objs(\n",
    "    generat.ImageObjsParams(\n",
    "        max_fluor=100,\n",
    "        max_num_objects=13,\n",
    "        max_radius=12,\n",
    "        min_radius=6,\n",
    "        ncols=Y,\n",
    "        nrows=X,\n",
    "    )\n",
    ")\n",
    "frame = generat.gen_frame(objects, bias=np.ones((Y, X)), flat=None, sky=7, noise_sd=6)\n",
    "segmentation.bg_refine_iteratively(frame, probplot=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8ea8996-7fdc-4249-8c00-9a8bb8315ee3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#         Method {'arcsinh', 'entropy', 'adaptive', 'li_adaptive', 'li_li'} used for the\n",
    "bg_result = segmentation.bg(frame, bg_params=segmentation.BgParams(kind=\"adaptive\"))\n",
    "bg_result.bg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8fc4f41-b387-4e63-b362-dda1da3a262d",
   "metadata": {},
   "outputs": [],
   "source": [
    "segmentation.bg_refine_iteratively(\n",
    "    frame, bgmax=segmentation._bgmax(frame, densityplot=True), probplot=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a44d75d8-1fbc-4a07-a8ba-cdf6f7ad01a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "r2 = utils.bg2(frame, bgmax=30)\n",
    "\n",
    "plt.plot(r2[3])\n",
    "r2[:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8cbd54f-9a57-4138-a290-c1be396037fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "lim = np.arcsinh(frame)\n",
    "plt.imshow(lim)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f1e2f34-10cf-49df-962c-1cd8a363733c",
   "metadata": {},
   "outputs": [],
   "source": [
    "lim = ndimage.percentile_filter(lim, 80, size=10)\n",
    "plt.imshow(lim)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d422ad13-2717-4aa1-bf02-24757295c220",
   "metadata": {},
   "outputs": [],
   "source": [
    "bg_result = segmentation.bg_refine_iteratively(frame)\n",
    "m0 = frame < bg_result.bg + 1 * bg_result.sd\n",
    "plt.imshow(m0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b29ce9c1-03cb-4988-b33a-c2e9d759b4df",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = segmentation.prob(frame, bg_result.bg, bg_result.sd)\n",
    "plt.imshow(p)\n",
    "plt.colorbar()\n",
    "plt.figure()\n",
    "plt.plot(p[100, :], \"o\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97491257-3911-45b2-95d6-9874ebab0be2",
   "metadata": {},
   "outputs": [],
   "source": [
    "arcsinh_perc = 80\n",
    "radius = 10\n",
    "perc = 0.1\n",
    "\n",
    "lim = np.arcsinh(frame)\n",
    "lim = ndimage.percentile_filter(lim, arcsinh_perc, size=radius)\n",
    "thr = (1 - perc) * lim.min() + perc * lim.max()\n",
    "m = lim < thr\n",
    "bgv = frame[m]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "460c5025-7d67-4eb8-96a7-178672938f00",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(131)\n",
    "plt.hist(p.ravel(), bins=99)\n",
    "# plt.semilogy()\n",
    "plt.subplot(132)\n",
    "skimage.io.imshow(p)\n",
    "plt.subplot(133)\n",
    "skimage.io.imshow(p > 0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2e9d3a9-735d-431d-8fe4-f06b0a68754e",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob_frame = segmentation.geometric_mean_filter(p, 1.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e0acadd-7dd7-4312-8c62-820775ab6967",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = prob_frame > 0.0029\n",
    "segmentation.fit_gaussian(frame[mask]), scipy.stats.distributions.norm.fit(frame[mask])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bfea458-7a05-4487-96d3-a0b9fa6c3d31",
   "metadata": {},
   "outputs": [],
   "source": [
    "bg, sd = bg_result.bg, bg_result.sd\n",
    "mgeo = segmentation.geometric_mean_filter(utils.prob(frame, bg, sd), 1) > 0.01\n",
    "# mgeo = skimage.filters.median(utils.prob(frame, bg,  sd)) > 0.01\n",
    "\n",
    "mgeo = (\n",
    "    ndimage.percentile_filter(utils.prob(frame, bg, sd), percentile=1, size=2) > 0.005\n",
    ")\n",
    "# mgeo = ndimage.uniform_filter(utils.prob(frame, bg222,  sd222), size=1) > 0.005\n",
    "# mgeo = ndimage.gaussian_filter(utils.prob(frame, bg222,  sd222), 0.25) > 0.005\n",
    "skimage.io.imshow(mgeo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f05146-19ee-4deb-bcad-77aa69549dd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "objs = generat.gen_objs(\n",
    "    generat.ImageObjsParams(\n",
    "        max_fluor=15,\n",
    "        max_num_objects=20,\n",
    "        max_radius=12,\n",
    "        min_radius=6,\n",
    "        ncols=264,\n",
    "        nrows=64,\n",
    "    )\n",
    ")\n",
    "frame = generat.gen_frame(objs, None, None, sky=10, noise_sd=6)\n",
    "\n",
    "bgres = nima.bg(frame.astype(\"float\"))\n",
    "# plt.hist(bgs, bins=8)\n",
    "bgres.bg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8ce2981-13d8-4049-adb0-c3f86de82393",
   "metadata": {},
   "outputs": [],
   "source": [
    "segmentation.bg_refine_iteratively(frame, probplot=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d31e6934-c55d-490d-a4fb-a105f9ef1dae",
   "metadata": {},
   "source": [
    "simulation for:\n",
    "- varying error\n",
    "- varying number of cells\n",
    "- varying intensity of cells\n",
    "- in the absence of bias correction\n",
    "- in the absence of flat correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d53cfde-3f5e-4003-98fa-fd6fa0d62cfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "objs_pars = generat.ImageObjsParams(max_num_objects=40, max_fluor=15, nrows=64)\n",
    "objs_pars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54df41f6-480d-45f3-8a0d-3ad6b0f423e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "utils.bg2(frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc4d574d-18ba-479c-82d8-4135bb389e12",
   "metadata": {},
   "outputs": [],
   "source": [
    "bg_arcsinh = []\n",
    "bg_entropy = []\n",
    "bg_adaptive = []\n",
    "bg_liadaptive = []\n",
    "bg_lili = []\n",
    "bg_invyen = []\n",
    "bg_ = []\n",
    "bg2_ = []\n",
    "\n",
    "\n",
    "def safe_call(func, *args, **kwargs):  # noqa: ANN002,ANN003\n",
    "    \"\"\"Wrap a function.\n",
    "\n",
    "    This wrapper attempts to call the given function with the provided arguments.\n",
    "    If the function call raises an exception, it returns NaN.\n",
    "    \"\"\"\n",
    "    try:\n",
    "        return func(*args, **kwargs)\n",
    "    except ZeroDivisionError:\n",
    "        print(\"Division by zero encountered.\")\n",
    "        return np.nan\n",
    "    except Exception as e:  # noqa: BLE001\n",
    "        print(f\"An error occurred: {e}\")\n",
    "        return np.nan\n",
    "\n",
    "\n",
    "objs_pars = generat.ImageObjsParams(max_num_objects=10, max_fluor=100, nrows=128)\n",
    "noise_sd = 4\n",
    "\n",
    "for _ in range(1):\n",
    "    print(_)\n",
    "    objs = generat.gen_objs(objs_pars)\n",
    "    frame = generat.gen_frame(objs, None, None, sky=10, noise_sd=noise_sd)\n",
    "    bg_arcsinh.append(\n",
    "        safe_call(nima.bg, frame, segmentation.BgParams(kind=\"arcsinh\")).bg\n",
    "    )\n",
    "    bg_entropy.append(\n",
    "        safe_call(nima.bg, frame, segmentation.BgParams(kind=\"entropy\")).bg\n",
    "    )\n",
    "    bg_adaptive.append(\n",
    "        safe_call(nima.bg, frame, segmentation.BgParams(kind=\"adaptive\")).bg\n",
    "    )\n",
    "    bg_liadaptive.append(\n",
    "        safe_call(nima.bg, frame, segmentation.BgParams(kind=\"li_adaptive\")).bg\n",
    "    )\n",
    "    bg_lili.append(safe_call(nima.bg, frame, segmentation.BgParams(kind=\"li_li\")).bg)\n",
    "    # bg_invyen.append(safe_call(nima.bg, frame, segmentation.BgParams(kind=\"inverse_yen\")).bg)\n",
    "    bg_.append(safe_call(segmentation.bg_refine_iteratively, frame).bg)\n",
    "    bg2_.append(safe_call(utils.bg2, frame, bgmax=500)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d28a78ec-928f-49f2-8937-9259b6c94b39",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sim2df(num_of_repeats, noise_sd, objs_pars):\n",
    "    bg_arcsinh = []\n",
    "    bg_entropy = []\n",
    "    bg_adaptive = []\n",
    "    bg_liadaptive = []\n",
    "    bg_lili = []\n",
    "    bg_invyen = []\n",
    "    bg_ = []\n",
    "    bg2_ = []\n",
    "\n",
    "    for _ in range(num_of_repeats):\n",
    "        # print(_)\n",
    "        objs = generat.gen_objs(objs_pars)\n",
    "        frame = generat.gen_frame(objs, None, None, sky=10, noise_sd=noise_sd)\n",
    "        bg_arcsinh.append(\n",
    "            safe_call(nima.bg, frame, segmentation.BgParams(kind=\"arcsinh\")).bg\n",
    "        )\n",
    "        bg_entropy.append(\n",
    "            safe_call(nima.bg, frame, segmentation.BgParams(kind=\"entropy\")).bg\n",
    "        )\n",
    "        bg_adaptive.append(\n",
    "            safe_call(nima.bg, frame, segmentation.BgParams(kind=\"adaptive\")).bg\n",
    "        )\n",
    "        bg_liadaptive.append(\n",
    "            safe_call(nima.bg, frame, segmentation.BgParams(kind=\"li_adaptive\")).bg\n",
    "        )\n",
    "        bg_lili.append(\n",
    "            safe_call(nima.bg, frame, segmentation.BgParams(kind=\"li_li\")).bg\n",
    "        )\n",
    "        bg_invyen.append(\n",
    "            safe_call(nima.bg, frame, segmentation.BgParams(kind=\"inverse_yen\")).bg\n",
    "        )\n",
    "        bg_.append(safe_call(segmentation.bg_refine_iteratively, frame).bg)\n",
    "        bg2_.append(safe_call(utils.bg2, frame, bgmax=800)[0])\n",
    "\n",
    "    df = pd.DataFrame(\n",
    "        np.column_stack(\n",
    "            (\n",
    "                bg_arcsinh,\n",
    "                bg_entropy,\n",
    "                bg_adaptive,\n",
    "                bg_liadaptive,\n",
    "                bg_lili,\n",
    "                bg_invyen,\n",
    "                bg_,\n",
    "                bg2_,\n",
    "            )\n",
    "        ),\n",
    "        columns=[\n",
    "            \"arcsinh\",\n",
    "            \"entropy\",\n",
    "            \"adaptive\",\n",
    "            \"li_adaptive\",\n",
    "            \"li li\",\n",
    "            \"inv_yen\",\n",
    "            \"bg\",\n",
    "            \"bg2\",\n",
    "        ],\n",
    "    )\n",
    "    df[\"sd\"] = noise_sd\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "517b0b4f-c685-4d35-a72d-91d1dd15d485",
   "metadata": {},
   "outputs": [],
   "source": [
    "objs_pars = generat.ImageObjsParams(max_num_objects=10, max_fluor=100, nrows=128)\n",
    "# df_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4b65815-3712-4e9d-b461-b3f35cb9d891",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_all = pd.DataFrame()\n",
    "for noise_sd in [1, 2, 3, 4]:\n",
    "    df = sim2df(num_of_repeats=49, objs_pars=objs_pars, noise_sd=noise_sd)\n",
    "    df_all = pd.concat([df_all, df])\n",
    "\n",
    "f = df_all.boxplot(vert=False, showfliers=0, by=\"sd\")\n",
    "plt.xlim((5, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72b66bf4-4dc1-4f69-96ee-55061a7263e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8c6065a-194b-4b1b-8524-74737542f893",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfmelted = df_all.melt(ignore_index=1, id_vars=\"sd\")\n",
    "dfmelted.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "446029fd-92d5-4154-8af0-962d788dd542",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50fc1ef8-a6fe-48f8-b3cf-ebe5d0a0d233",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.barplot(dfmelted, x=\"variable\", y=\"value\", hue=\"sd\", dodge=1, palette=\"Blues_d\")\n",
    "plt.ylim(4, 10.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eeabb350-3270-49af-91f2-7bf43d74c968",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.stripplot(\n",
    "    dfmelted,\n",
    "    x=\"variable\",\n",
    "    y=\"value\",\n",
    "    hue=\"sd\",\n",
    "    dodge=1,\n",
    "    jitter=0.2,\n",
    "    alpha=0.6,\n",
    "    palette=\"Reds_d\",\n",
    ")\n",
    "sns.boxenplot(\n",
    "    dfmelted, x=\"variable\", y=\"value\", hue=\"sd\", dodge=1, alpha=0.6, palette=\"Reds_d\"\n",
    ")\n",
    "\n",
    "plt.ylim(7.0, 11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9709840f-d940-42db-9278-d21e464e6ef4",
   "metadata": {},
   "outputs": [],
   "source": [
    "skimage.io.imshow(frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6811977b-58e8-465f-9294-a4cbbbf3d494",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "sns.regplot(pd.DataFrame({\"x\": bg_arcsinh, \"y\": bg_}), x=\"x\", y=\"y\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6c775ad-5ada-42cd-966e-e211925b4816",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate a single small object\n",
    "small_object = generat.gen_object(nrows=12, ncols=12, min_radius=3, max_radius=4)\n",
    "\n",
    "# Plot the generated object\n",
    "plt.imshow(small_object, cmap=\"gray\")\n",
    "plt.title(\"Single Small Object\")\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Y\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7c5ea2c-6b32-4d4d-89c6-8cba474c8400",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.signal\n",
    "\n",
    "# Convolve the small object with the flat image\n",
    "convolved_image = scipy.signal.convolve2d(flat, small_object, mode=\"same\")\n",
    "\n",
    "# Plot the convolved image\n",
    "plt.imshow(convolved_image, cmap=\"gray\")\n",
    "plt.colorbar()\n",
    "plt.title(\"Convolved Image\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7b8cc42-7b2e-4b94-bbd5-e71f511b0376",
   "metadata": {},
   "outputs": [],
   "source": [
    "flat.shape[1] - small_object.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71c40ff7-0284-4bae-b375-6ad6c632e1bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of frames in the stack\n",
    "num_frames = 10000\n",
    "\n",
    "# Initialize an empty stack to store the frames\n",
    "stack = np.zeros_like(flat)\n",
    "\n",
    "# Iterate over each frame in the stack\n",
    "for _ in range(num_frames):\n",
    "    # Generate random coordinates for the position of the small object within the flat image\n",
    "    x_pos = nrg.integers(0, flat.shape[1] - small_object.shape[1])\n",
    "    y_pos = nrg.integers(0, flat.shape[0] - small_object.shape[0])\n",
    "\n",
    "    # Add the small object to the flat image at the random position\n",
    "    flat_image_with_object = flat.copy()\n",
    "    flat_image_with_object[\n",
    "        y_pos : y_pos + small_object.shape[0], x_pos : x_pos + small_object.shape[1]\n",
    "    ] += small_object\n",
    "\n",
    "    # Add the frame with the small object to the stack\n",
    "    stack += flat_image_with_object\n",
    "\n",
    "# Plot the summed stack\n",
    "estimated = stack / stack.mean()\n",
    "plt.imshow(estimated, cmap=\"gray\")\n",
    "plt.colorbar()\n",
    "plt.title(\"Summed Stack with Small Object\")\n",
    "plt.show()\n",
    "# plt.imshow(estimated - flat, cmap='gray')\n",
    "skimage.io.imshow(ndimage.gaussian_filter(estimated, sigma=3) - flat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ba9db07-f612-49cd-a9de-fb2bcb2d146c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the Fourier transform of the small object\n",
    "fourier_transform_obj = np.fft.fft2(small_object)\n",
    "\n",
    "# Calculate the magnitude spectrum of the Fourier transform\n",
    "magnitude_spectrum = np.abs(np.fft.fftshift(fourier_transform_obj))\n",
    "\n",
    "# Plot the magnitude spectrum\n",
    "plt.imshow(magnitude_spectrum, cmap=\"gray\")\n",
    "plt.colorbar(label=\"Magnitude\")\n",
    "plt.title(\"Magnitude Spectrum of Fourier Transform\")\n",
    "plt.xlabel(\"Frequency (kx)\")\n",
    "plt.ylabel(\"Frequency (ky)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45ddd792-d280-4ce7-94ab-c77c5a979ba3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the convolution theorem\n",
    "flat_fft = np.fft.fft2(stack)\n",
    "\n",
    "# Calculate the Fourier transform of the small object\n",
    "fourier_transform_obj = np.fft.fft2(small_object)\n",
    "\n",
    "# Pad the small object to match the shape of flat\n",
    "padded_obj = np.pad(\n",
    "    small_object,\n",
    "    (\n",
    "        (0, flat.shape[0] - small_object.shape[0]),\n",
    "        (0, flat.shape[1] - small_object.shape[1]),\n",
    "    ),\n",
    "    mode=\"constant\",\n",
    ")\n",
    "\n",
    "# Calculate the Fourier transform of the padded small object\n",
    "fourier_transform_padded_obj = np.fft.fft2(padded_obj)\n",
    "\n",
    "# Calculate the Fourier transform of the flat image\n",
    "flat_fft = np.fft.fft2(flat)\n",
    "\n",
    "# Perform element-wise division\n",
    "result_fft = np.fft.ifftshift(\n",
    "    np.fft.ifft2(np.fft.fftshift(flat_fft / fourier_transform_padded_obj))\n",
    ")\n",
    "# result_fft = np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(flat_fft / fourier_transform_obj)))\n",
    "\n",
    "# Take the real part to get rid of any numerical artifacts\n",
    "result = np.real(result_fft)\n",
    "\n",
    "# Plot the resulting flat image\n",
    "plt.imshow(result, cmap=\"gray\")\n",
    "plt.colorbar(label=\"Intensity\")\n",
    "plt.title(\"Resulting Flat Image\")\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Y\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd3b608f-869f-4a89-8bdd-d33de2beb84a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from nima import generat\n",
    "\n",
    "flat = generat.gen_flat()\n",
    "bias = 7 * np.ones((128, 128))\n",
    "\n",
    "objs = generat.gen_objs(generat.ImageObjsParams(max_fluor=20, max_num_objects=80))\n",
    "frame = generat.gen_frame(objs, bias=bias, flat=flat, noise_sd=2, sky=7)\n",
    "\n",
    "# Plot the frame\n",
    "plt.imshow(frame, cmap=\"viridis\", origin=\"lower\")\n",
    "plt.colorbar(label=\"Intensity\")\n",
    "plt.title(\"Simulated Image Frame without Bias\")\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Y\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32dd6e21-f9af-4b1e-b6a8-c6ff1acc17ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm\n",
    "\n",
    "# Generate a stack of frames\n",
    "num_frames = 100\n",
    "frame_stack = []\n",
    "for _ in tqdm(range(num_frames), desc=\"Generating Frames\"):\n",
    "    objs = generat.gen_objs(generat.ImageObjsParams(max_fluor=20, max_num_objects=80))\n",
    "    frame = generat.gen_frame(objs, bias=bias, flat=flat, noise_sd=2, sky=7)\n",
    "    frame_stack.append(frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08d8a3f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "p999 = partial(np.percentile, q=99.7)\n",
    "p999.__name__ = \"percentile 99.9%\"\n",
    "\n",
    "\n",
    "def diff_plot(im, flat, title):\n",
    "    f, axs = plt.subplots(1, 2)\n",
    "    diff = im / im.mean() - flat\n",
    "    skimage.io.imshow(diff, ax=axs[0])\n",
    "    axs[1].hist(diff.ravel())\n",
    "    f.suptitle(title)\n",
    "    return diff.mean(), diff.std()\n",
    "\n",
    "\n",
    "def prj_plot(t_prj, title, sigma=128 / 11):\n",
    "    im = ndimage.gaussian_filter(t_prj, sigma=sigma)\n",
    "    return diff_plot(im, flat, title)\n",
    "\n",
    "\n",
    "def prj(stack, func):\n",
    "    t_prj = func(stack, axis=0)\n",
    "    return prj_plot(t_prj, func.__name__)\n",
    "\n",
    "\n",
    "prj(frame_stack, np.max)\n",
    "prj(frame_stack, p999)\n",
    "prj(frame_stack, np.mean)\n",
    "prj(frame_stack, np.median)\n",
    "prj(frame_stack, np.min)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d38fc303",
   "metadata": {},
   "outputs": [],
   "source": [
    "objs = generat.gen_objs(generat.ImageObjsParams(max_fluor=20, max_num_objects=8))\n",
    "frame = generat.gen_frame(objs, bias=bias, flat=flat)\n",
    "plt.imshow(frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5967d57b",
   "metadata": {},
   "outputs": [],
   "source": [
    "bias = np.zeros((128, 128))\n",
    "flat = np.ones((128, 128))\n",
    "\n",
    "stack = np.stack(\n",
    "    [\n",
    "        generat.gen_frame(\n",
    "            generat.gen_objs(generat.ImageObjsParams(max_fluor=20, max_num_objects=20)),\n",
    "            None,\n",
    "            None,\n",
    "            noise_sd=8,\n",
    "            sky=20,\n",
    "        )\n",
    "        for i in range(1000)\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0c2e9e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "stat_bg = []\n",
    "stat_sd = []\n",
    "for s in stack[:]:\n",
    "    bg, sd = segmentation.iteratively_refine_background(s)[:2]\n",
    "    stat_bg.append(bg)\n",
    "    stat_sd.append(sd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21c2d1e6-703a-4487-86c0-24ee835fe054",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(121)\n",
    "plt.hist(stat_bg)\n",
    "print(np.mean(stat_bg), np.std(stat_bg))\n",
    "plt.subplot(122)\n",
    "plt.hist(stat_sd)\n",
    "np.mean(stat_sd), np.std(stat_sd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0214305d-6f19-41dd-9159-5eb63dd5a935",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(stack[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42ad7ca8",
   "metadata": {},
   "source": [
    "## what is the best projection for flat calculation? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ff52162",
   "metadata": {},
   "outputs": [],
   "source": [
    "bias = generat.gen_bias()\n",
    "flat = generat.gen_flat()\n",
    "stack = np.stack(\n",
    "    [\n",
    "        generat.gen_frame(\n",
    "            generat.gen_objs(generat.ImageObjsParams(max_fluor=20)),\n",
    "            bias,\n",
    "            flat,\n",
    "            noise_sd=1,\n",
    "            sky=2,\n",
    "        )\n",
    "        for i in range(1000)\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cd1e9f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def splot(stack, num=4):\n",
    "    _f, axs = plt.subplots(1, num)\n",
    "    for i in range(num):\n",
    "        axs[i].imshow(stack[nrg.integers(0, len(stack))])\n",
    "\n",
    "\n",
    "splot(stack)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca669ce4",
   "metadata": {},
   "outputs": [],
   "source": [
    "prj(stack, np.max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72dfef8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "prj(stack, np.mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "433eeb3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "prj(stack, np.median)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f14fca2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "p999 = partial(np.percentile, q=99.9)\n",
    "p999.__name__ = \"percentile 99.9%\"\n",
    "\n",
    "prj(stack, p999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b151b841-d864-44fc-89f8-243b5373123b",
   "metadata": {},
   "outputs": [],
   "source": [
    "im = np.mean(\n",
    "    ndfilters.median_filter(\n",
    "        da.from_array(stack[:100] - bias), size=(32, 16, 16)\n",
    "    ).compute(),\n",
    "    axis=0,\n",
    ")\n",
    "prj_plot(im, \"dd\", sigma=7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a3b46a1",
   "metadata": {},
   "source": [
    "#### Knowing the Bias."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a54c4d1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "prj(stack - bias, p999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "429d6e5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "prj(stack - bias, np.mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "752fdda7",
   "metadata": {},
   "source": [
    "### Using fg and bg masks?\n",
    "\n",
    "And assuming we know the bias of the camera."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d824423",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mask_plane(plane, bg_ave=2, bg_std=1.19, erf_pvalue=0.01):\n",
    "    p = utils.prob(plane, bg_ave, bg_std)\n",
    "    p = ndimage.median_filter(p, size=2)\n",
    "    mask = p > erf_pvalue\n",
    "    mask = skimage.morphology.remove_small_holes(mask)\n",
    "    return np.ma.masked_array(plane, mask=~mask), np.ma.masked_array(plane, mask=mask)\n",
    "\n",
    "\n",
    "plt.imshow(mask_plane(stack[113], *utils.bg(stack[113]))[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9c8e6d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "bgs, fgs = list(\n",
    "    zip(*[mask_plane(s - bias, *utils.bg(s - bias)) for s in stack], strict=False)\n",
    ")\n",
    "\n",
    "splot(bgs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70ff0a9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_prj = np.ma.mean(np.ma.stack(bgs), axis=0)\n",
    "prj_plot(t_prj, \"Bg mean\", sigma=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d85c55fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_prj = np.ma.max(np.ma.stack(fgs), axis=0)\n",
    "prj_plot(t_prj, \"Fg max (bias known)\", sigma=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "981bea30",
   "metadata": {},
   "outputs": [],
   "source": [
    "bgs, fgs = list(zip(*[mask_plane(s, *utils.bg(s)) for s in stack], strict=False))\n",
    "\n",
    "bg_prj1 = np.ma.mean(np.ma.stack(bgs[:]), axis=0)\n",
    "fg_prj1 = np.ma.mean(np.ma.stack(fgs[:]), axis=0)\n",
    "im = fg_prj1 - bg_prj1\n",
    "diff_plot(ndimage.gaussian_filter(im, 1), flat, \"Bg mean - fg mean\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9f7998e",
   "metadata": {},
   "outputs": [],
   "source": [
    "bg_prj = np.ma.mean(bgs, axis=0)\n",
    "fg_prj = np.ma.max(fgs, axis=0)\n",
    "# im = ndimage.median_filter(bg_prj-fg_prj, size=60) #- 2 * flat\n",
    "im = ndimage.gaussian_filter(bg_prj - fg_prj, sigma=14)  # - 2 * flat\n",
    "\n",
    "diff_plot(im, flat, \"m\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2e198b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_prj = np.ma.max(fgs, axis=0)\n",
    "prj_plot(t_prj, \"Fg MAX\", sigma=13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ed3a92f",
   "metadata": {},
   "outputs": [],
   "source": [
    "eflat = bg_prj - fg_prj\n",
    "eflat /= eflat.mean()\n",
    "eflat = ndimage.gaussian_filter(eflat, sigma=13)\n",
    "\n",
    "diff_plot(eflat, flat, \"eflat\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "164d4d74",
   "metadata": {},
   "source": [
    "## When bias and flat are unknown...\n",
    "\n",
    "- bias = bg_prj - sky * flat\n",
    "- bias = fg_prj - flat\n",
    "\n",
    "sky * flat - flat = bg_prj - fg_prj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d100277a",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_plot((bg_prj1 - bias) / 2, flat, \"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d457e605",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow((im - bias) / (im - bias).mean() - flat)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fd25fb1",
   "metadata": {},
   "source": [
    "## cfr. nima.bg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "959a7db0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# r = nima.bg((stack[113] - bias) / flat)\n",
    "r = nima.bg(stack[111])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9104c54f",
   "metadata": {},
   "outputs": [],
   "source": [
    "r[1].mean(), r[1].std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "234f8826",
   "metadata": {},
   "outputs": [],
   "source": [
    "utils.bg(stack[111])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b920d5ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "bias.mean() + 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66362d94",
   "metadata": {},
   "source": [
    "## geometric mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6745872c",
   "metadata": {},
   "outputs": [],
   "source": [
    "vals = [0.8, 0.1, 0.3, 0.1, 0.8, 0.8, 0.8, 0.1, 0.8]\n",
    "\n",
    "np.median(vals), scipy.stats.gmean(vals), np.mean(vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af2b39c9-0bb3-4e82-a56b-c065bfea1c0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = vals[0]\n",
    "p *= vals[1]\n",
    "p *= vals[2]\n",
    "p *= vals[3]\n",
    "p *= vals[4]\n",
    "p *= vals[5]\n",
    "p *= vals[6]\n",
    "p *= vals[7]\n",
    "p *= vals[8]\n",
    "np.power(p, 1 / 9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ceecb8d3-959c-4f89-b27d-d99a62ae591e",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.exp(np.log(np.array(vals)).sum() / 9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18617c34-2438-469c-bd48-d99faf9721be",
   "metadata": {},
   "outputs": [],
   "source": [
    "vv = np.array(vals).reshape(3, 3)\n",
    "vv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eab3ed7a-0d07-40cf-9ddc-86112f045d30",
   "metadata": {},
   "outputs": [],
   "source": [
    "segmentation.geometric_mean_filter(vv, 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8440837-3adf-41f3-8ab1-8eb795e067ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = skimage.morphology.disk(1.0).astype(float)\n",
    "n = np.sum(kernel)  # Total weight, or number of ones in the kernel\n",
    "print(n)\n",
    "\n",
    "plt.imshow(kernel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fae06173",
   "metadata": {},
   "outputs": [],
   "source": [
    "(0.8 * 0.8 * 0.1 * 0.1 * 0.1) ** (1 / 5)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}