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
+}