From e487d92307b99882bf57b39ff9628c23c8014383 Mon Sep 17 00:00:00 2001 From: Guillaume Raffy Date: Thu, 10 Mar 2022 12:09:15 +0100 Subject: [PATCH] investigated a method to estimate the dark for images that have been captured with missing reference dark images. This jupyter notebook is a proof of concept that shows that the method works, at least on synthetic images --- doc/find_dark.ipynb | 405 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 405 insertions(+) create mode 100644 doc/find_dark.ipynb diff --git a/doc/find_dark.ipynb b/doc/find_dark.ipynb new file mode 100644 index 0000000..b1c98c0 --- /dev/null +++ b/doc/find_dark.ipynb @@ -0,0 +1,405 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# how to estimate the white and dark from fluorescent image sequences" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We suppose that the captured images $I(t,s)$ follow this model:\n", + "\n", + "\\begin{equation*}\n", + "I(t,s) = D + F(t)L(s)\n", + "\\label{eq:image_model} \\tag{1}\n", + "\\end{equation*}\n", + "\n", + "where:\n", + "- $t$ is the time\n", + "- $s$ is the site (pixel location in an image)\n", + "- $D$ is a constant dark value\n", + "- $F(t)$ is the function that represents the decay of the florenscence over time $t$\n", + "- $L(s)$ models the ratio of light received by the sensor for each pixel site $s$\n", + "\n", + "\n", + "\n", + "The goal of this study is to see if it's possible to estimate $D$, $F$ and $L$ from $I$. If this works, this will allow to recover the dark value for series such as `soleil2016/GGHL_rDGL_SGF55_lambda_Em_cinsuite_1_Pos0` for which the dark value hasn't been captured\n", + "\n", + "---\n", + "**NOTE**\n", + "\n", + "In this study, we assume that the image $I$ is 1-dimensional, which means that a pixel site $s$ is just its $x$ corrdinate instead of the couple $x,y$. We work with 1 dmensional images for the sake of easiness, as this method would extend to 2 dimensional images.\n", + "\n", + "---\n", + "\n", + "## building a synthetic test image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as pyplot\n", + "from scipy import optimize\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "real_dark = 42.0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "num_times = 5\n", + "t = np.arange(0, num_times)\n", + "print(t)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fluo_decrease_speed = 3.0\n", + "fluo0 = 1.0\n", + "fluo = np.exp(-t/num_times * fluo_decrease_speed) * fluo0\n", + "print(fluo)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pyplot.plot(fluo)\n", + "pyplot.xlabel('$t$')\n", + "pyplot.ylabel('$f(t)$')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$L(s)$ is the lens signal for each pixel site $s$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "num_sites = 6\n", + "s = np.arange(0, num_sites)\n", + "lens = (np.sin(((s - num_sites*0.5) / num_sites + 0.25) * 2.0 * 3.14159 ) + 1.0) * 0.4 + 0.2\n", + "print(lens)\n", + "pyplot.plot(lens)\n", + "pyplot.xlabel('$s$ (pixel position)')\n", + "pyplot.ylabel('$L(s)$')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we have a $L(s)$ signal, a $F(t)$ signal, and a dark value $D$, we are able to create a synthetic sequence of images $I(t,s)$ using the model defined by equation $\\eqref{eq:image_model}$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "i = lens.reshape(1, num_sites) * fluo.reshape(num_times, 1) + real_dark\n", + "print(i)\n", + "plots = pyplot.plot(i.transpose())\n", + "pyplot.xlabel('$s$ (pixel position)')\n", + "pyplot.ylabel('$I(t,s)$')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## estimation of D, L and F from I\n", + "\n", + "$$A_n(s) = I(t_{n+1},s) - I(t_n,s)$$\n", + "$$A_n(s) = [D + F(t_{n+1})L(s)] - [D + F(t_n)L(s)]$$\n", + "$$A_n(s) = F(t_{n+1})L(s)- F(t_n)L(s)$$\n", + "$$A_n(s) = L(s)(F(t_{n+1})- F(t_n))$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = i[1:num_times] - i[0:num_times-1]\n", + "print(A)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$R_n(s) = \\frac{A_{n+1}(s)}{A_n(s)}$$\n", + "$$R_n(s) = \\frac{F(t_{n+2})-F(t_{n+1})}{F(t_{n+1})-F(t_{n})}$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "R = A[1:num_times-1] / A[0:num_times-2]\n", + "print(R)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$P_n(s) = \\Pi_{k=1}^{n}R_k(s)$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "P = np.multiply.accumulate(R)\n", + "print(P)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$P_n(s) = \\frac{F(t_3)-F(t_2)}{F(t_2)-F(t_1)} \\times \\frac{F(t_4)-F(t_3)}{F(t_3)-F(t_2)} \\times \\cdots \\times\\frac{F(t_{n+2})-F(t_{n+1})}{F(t_{n+1})-F(t_{n})}$$\n", + "$$P_n(s) = \\frac{F(t_{n+2})-F(t_{n+1})}{F(t_2)-F(t_1)}$$\n", + "$$F(t_{n+2})-F(t_{n+1}) = P_n(s) . (F(t_2)-F(t_1))$$\n", + "$$A_n(s) = L(s)(P_{n-1}(s) . (F(t_2)-F(t_1)))$$\n", + "$$L(s) = \\frac{A_n(s)}{P_{n-1}(s) . (F(t_2)-F(t_1))}$$\n", + "Let's define $f_1=F(t_2)-F(t_1)$\n", + "\n", + "Then we have:\n", + "$$L(s) = \\frac{A_n(s)}{P_{n-1}(s)f_1}$$\n", + "At this point, the only unknown to estimate $L$ is $f_1$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def compute_lens(a, p, f_1):\n", + " return a / p / f_1\n", + "\n", + "print(compute_lens(a=A[1], p=P[0], f_1=1.0))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def compute_lenses(a, p, f_1):\n", + " return a[1:] / p / f_1\n", + "\n", + "print(compute_lenses(a=A, p=P, f_1=-1.0))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def estimate_fluos(f_1, F_1, R):\n", + " \"\"\"estimate the fluo signal for each site\n", + " Args:\n", + " F_1: the fluo signal for the first frame (fluo[0])\n", + " f_1: fluo[1] - fluo[0]\n", + " \"\"\"\n", + " fluos = np.zeros((R.shape[0]+2, R.shape[1]))\n", + " fluos[0] = F_1\n", + " F_2 = f_1 + F_1\n", + " fluos[1] = F_2\n", + " for i in range(R.shape[0]):\n", + " fi = fluos[i+1] - fluos[i]\n", + " fj = R[i] * fi\n", + " Fi = fluos[i+1]\n", + " Fj = Fi + fj\n", + " fluos[i+2] = Fj\n", + " return fluos\n", + "\n", + "print(estimate_fluos(f_1=fluo[1]-fluo[0], F_1=fluo0, R=R))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that `estimate_fluos` works properly, since it's able to reconstruct the $F$ signal from $R$ and the right values for $F(2)-F(1)$ and $F(1)$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def estimate_dark_lens_fluo(I, A, P, R, f_1, F_1):\n", + " lenses = compute_lenses(A, P, f_1)\n", + " fluos = estimate_fluos(f_1, F_1, R)\n", + " lens = lenses[0]\n", + " fluo = fluos[:,0]\n", + " dark = I - lens.reshape(1, len(lens)) * fluo.reshape(len(fluo), 1)\n", + " return dark, lens, fluo\n", + "\n", + "print(estimate_dark(i, A, P, R, f_1=fluo[1]-fluo[0], F_1=fluo0))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that with the proper values for `f_1` and `F_1`, the dark value is found at every site at every time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(estimate_dark(i, A, P, R, f_1=fluo[1]-fluo[0], F_1=fluo0+0.5))\n", + "print(estimate_dark(i, A, P, R, f_1=fluo[1]-fluo[0]+0.5, F_1=fluo0))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And we can see that this is no longer the case if the values used for `f_1` and `F_1` are wrong. So, this gives the idea of finding the correct values of `f_1` and `F_1` by finding the values that result in a constant dark value" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def find_dark(I):\n", + " (num_times, num_sites) = I.shape\n", + " A = I[1:num_times] - I[0:num_times-1]\n", + " R = A[1:num_times-1] / A[0:num_times-2]\n", + " P = np.multiply.accumulate(R)\n", + " \n", + " def function_to_minimize(x):\n", + " f_1 = x[0]\n", + " F_1 = x[1]\n", + " dark_matrix = estimate_dark(I, A, P, R, f_1, F_1)\n", + " error = np.var(dark_matrix)\n", + " return error\n", + " \n", + " minimize_result = optimize.minimize(fun=function_to_minimize, x0=(-0.7, 0.5), method = 'Nelder-Mead')\n", + " if not minimize_result.success:\n", + " print(type(minimize_result))\n", + " print(minimize_result)\n", + " raise Exception(minimize_result.message)\n", + " f_1 = minimize_result.x[0]\n", + " F_1 = minimize_result.x[1]\n", + " print('estimated value of f_1 : %f' % f_1)\n", + " print('estimated value of F_1 : %f' % F_1)\n", + " dark, lens, fluo = estimate_dark_lens_fluo(I, A, P, R, f_1, F_1)\n", + " return dark, lens, fluo\n", + "\n", + "print('exact value of f_1 : %f' % (fluo[1]-fluo[0]))\n", + "print('exact value of F_1 : %f' % fluo[0])\n", + "\n", + "estimated_dark, estimated_lens, estimated_fluo = find_dark(i)\n", + "print('estimated dark value : ', estimated_dark)\n", + "print('exact lens : ', lens)\n", + "print('estimated lens : ', estimated_lens)\n", + "pyplot.plot(lens)\n", + "pyplot.plot(estimated_lens)\n", + "pyplot.xlabel('$s$ (pixel position)')\n", + "pyplot.ylabel('$L(s)$')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print('exact fluo :', fluo)\n", + "print('estimated fluo', estimated_fluo)\n", + "pyplot.plot(fluo)\n", + "pyplot.plot(estimated_fluo)\n", + "pyplot.xlabel('$t$ (time)')\n", + "pyplot.ylabel('$F(s)$')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Even if the values of $f_1$ and $F_1$ are not correct, the exact value of for the dark $D$ has been found, which proves that this method is able to estimate the dark. We can also observe that $F$ and $L$ signals are recovered up to a scale factor." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}