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
This commit is contained in:
parent
fcd9dde587
commit
e487d92307
|
|
@ -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
|
||||
}
|
||||
Loading…
Reference in New Issue