RoddierWaveFrontReconstruction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Roddier method for wavefront reconstruction\n",
    "\n",
    "## Introduction\n",
    "This notebook mostly aims at studying the Roddier method for wavefront reconstruction in optics. It was first introduced in the following paper:\n",
    "\"Wave-front reconstruction from defocused images and the testing of ground-based optical telescopes\" by Claude Roddier and Fransois Roddier, 1993"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"1200\"\n",
       "            height=\"800\"\n",
       "            src=\"doc/Optics/Wave-front_reconstruction_from_defocused_images_and_the_testing_of_ground-based_optical_telescopes.pdf\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7f2a0092ec70>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import IFrame\n",
    "IFrame(\"doc/Optics/Wave-front_reconstruction_from_defocused_images_and_the_testing_of_ground-based_optical_telescopes.pdf\", width=1200, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic optic concepts\n",
    "\n",
    "We will mostly consider stigmatic optical system here, however multiple shapes for pupil might be considered. Also, as the paper focuses on diffraction effects analysis, and wavefront analysis, we will mostly consider what happens on the optical axis\n",
    "\n",
    "### Object vs Image space\n",
    "From [wikipedia](https://en.wikipedia.org/wiki/Optical_space)\n",
    "Optical spaces are mathematical coordinate systems that facilitate the modelling of optical systems as mathematical transformations.\n",
    "An optical space is a mathematical coordinate system such as a Cartesian coordinate system associated with a refractive index. The analysis of optical systems is greatly simplified by the use of optical spaces which enable designers to place the origin of a coordinate system at any of several convenient locations. In the design of optical systems two optical spaces, **object space** and **image space**, are always employed\n",
    "\n",
    "Object space usually refers to the space between the source and the optical system (although it is not limited to it); image space is the space between the optical system and the image.\n",
    "\n",
    "It is common practice to designate the **horizontal axis** of an optical space as the **z-axis** with the positive direction left to right. Similarly, the **y-axis** is **vertical** with the positive direction upward.\n",
    "\n",
    "### Stigmatism\n",
    "From [wikipedia](https://en.wikipedia.org/wiki/Stigmatism)\n",
    "In geometric optics, stigmatism refers to the image-formation property of an optical system which focuses a single point source in object space into a single point in image space. Two such points are called a stigmatic pair of the optical system.\n",
    "\n",
    "\n",
    "### Conjugate focal plane\n",
    "From [wikipedia](https://en.wikipedia.org/wiki/Conjugate_focal_plane)\n",
    "In optics, a conjugate plane or conjugate focal plane of a given plane $P$, is the plane $P'$ such that points on $P$ are imaged on $P'$.\n",
    "\n",
    "![title](data/Conjugate_planes.jpg)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Roddier paper walkthrough\n",
    "\n",
    "### Basic experimental setup\n",
    "\n",
    "Irradiance data $I_1$ and $I_2$ are acquired at two different focal planes $P_1$ and $P_2$ that are symmetrically out of focus at a distance $l=l'$.\n",
    "\n",
    "Equivalently, in object space, the defocus images conjugate can be thought of as being at a distance $\\Delta z$ from the pupil plane $P$: one before the entrance pupil and one beyond the pupil. Here is a reproduction of the original paper:\n",
    "\n",
    "![title](data/roddier_scheme_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Irradiance transport equation / Transport of intensity equation (TIE)\n",
    "\n",
    "#### Introduction\n",
    "Although this equation is referred to as \"Irradiance transport equation\" in the original paper from Roddier, it is more often refered to as \"transport of intensity equation\", see wikipedia: https://en.m.wikipedia.org/wiki/Transport-of-intensity_equation\n",
    "\n",
    "Important note, TIE is different from ***radiative transport equation*** (you can find a nice introduction on [youtube](https://www.youtube.com/watch?v=I_90kOpj7J4&ab_channel=AaronParsons) )\n",
    "\n",
    "It has been extensively studied in the past, and still is a topic of interest for some physics researchers, we find interesting to cite two papers before going into more details for Roddier method.\n",
    "\n",
    "The first one, named ***Deterministic phase retrieval: a Green's function solution*** is an article from 1983 that aimed at studying the partial differential equation, and formulate a solution based on Green's function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#IFrame(\"doc/Optics/teague1983.pdf\", width=1200, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second one, named: ***On a universal solution to the transport-of-intensity equation*** presents an iterative approach to solve the pde, with stability considerations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#IFrame(\"doc/Optics/1912.07371.pdf\", width=1200, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The irradiance transport partial differential equation reads:\n",
    "\n",
    "\\begin{align*}\n",
    "  \\frac{\\partial I}{\\partial z} = -\\left( \\nabla I \\cdot \\nabla W + I \\nabla^2 W\\right)\n",
    "\\end{align*}\n",
    "\n",
    "Where\n",
    "* $I(x, y, z)$ is the distribution of the illumination along the beam, a scalar quantity usually expressed in TBD\n",
    "* $W(x, y)$ is the wave-front surface. A scalar quantity that can be expressed in nm or in number of wavelength when one study the behaviour of a specific monochromatic light source.\n",
    "* $z$ is a scalar that stands for the distance to pupil $P$ along axis $\\vec{z}$. origin of $z$ is a the pupil plane $P$\n",
    "* $\\nabla$ is a differential operator called gradient. its output is a vector. In our case it could be written $\\begin{pmatrix} \\frac{\\partial}{\\partial x} \\\\ \\frac{\\partial}{\\partial y}\\end{pmatrix}$\n",
    "* $\\nabla^2$ is also called Lapacian, it is a scalar and can be written $\\Delta = \\nabla \\cdot \\nabla = \\left( \\frac{\\partial}{\\partial x}\\right)^2 + \\left( \\frac{\\partial}{\\partial y}\\right)^2$\n",
    "\n",
    "We apply this equation to the pupil plane $(z = 0)$, where we assume the illumination to be fairly uniform and equal to $I_0$ inside the pupil and 0 outside. In this plane $\\nabla I = \\vec{0}$ everywhere but at the pupil edge, where we have the following expression:\n",
    "\\begin{align*}\n",
    "  \\nabla I = -I_0 \\vec{\\eta} \\delta_c\n",
    "\\end{align*}\n",
    "Where\n",
    "* $\\delta_c$, is a scalar field, described by a Dirac distribution around the pupil edge\n",
    "* $\\vec{\\eta}$ is a vector field normal to the edge and pointing outward, with unit length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Edge case\n",
    "Lets plugin the edge case description into original equation:\n",
    "\n",
    "\\begin{align*}\n",
    "  \\frac{\\partial I}{\\partial z} &= -\\left( -I_0 \\vec{\\eta} \\delta_c \\cdot \\nabla W + I_0 P \\nabla^2 W\\right) \\\\\n",
    "  &= -\\left( -I_0 \\frac{\\partial W}{\\partial \\eta} \\delta_c + I_0 P \\nabla^2 W\\right) \\\\\n",
    "  &= I_0 \\frac{\\partial W}{\\partial \\eta} \\delta_c - I_0 P \\Delta W \\\\\n",
    "\\end{align*}\n",
    "\n",
    "Where we define\n",
    "* $P(x, y)$ a scalar field that stands for the pupil indicator function, i.e it is equal to 1 inside the pupil and 0 outside"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Illumination at the defocused plane conjugates\n",
    "Lets simply express the illumination near the pupile plane (which edge case as been seen before), at the conjugate of the defocused planes, $I_1$ and $I_2$:\n",
    "\\begin{align*}\n",
    "  I_1 &= I_0 - \\frac{\\partial I}{\\partial z} \\Delta z\\\\\n",
    "  I_2 &= I_0 + \\frac{\\partial I}{\\partial z} \\Delta z\\\\\n",
    "\\end{align*}\n",
    "\n",
    "#### Sensor signal\n",
    "The quantity $S$, called the sensor signal, can be computed as:\n",
    "\\begin{align*}\n",
    "  &\\frac{I_1-I_2}{I_1+I_2}\\\\\n",
    "  &=\\frac{I_0 - \\frac{\\partial I}{\\partial z} \\Delta z - I_0 - \\frac{\\partial I}{\\partial z} \\Delta z}{I_0 - \\frac{\\partial I}{\\partial z} \\Delta z + I_0 + \\frac{\\partial I}{\\partial z} \\Delta z} \\\\\n",
    "  &= \\frac{-2\\frac{\\partial I}{\\partial z} \\Delta z}{2I_0}\\\\\n",
    "  &= -\\frac{1}{I_0}\\frac{\\partial I}{\\partial z} \\Delta z\n",
    "\\end{align*}\n",
    "It should be noted that, since in practice images are recorded in the image space, one has to invert (rotate by\n",
    "180 deg) the outside-focus image before computing $S$ (I guess this is why minus sign is later removed).\n",
    "\n",
    "We can now replace $\\frac{\\partial I}{\\partial z}$ in signal expression, with its expression that depends on wavefront:\n",
    "\\begin{align*}\n",
    "  S &= \\frac{1}{I_0} \\left( I_0 \\frac{\\partial W}{\\partial \\eta} \\delta_c - I_0 P \\Delta W \\right) \\Delta z\\\\\n",
    "  &= \\left( \\frac{\\partial W}{\\partial \\eta} \\delta_c - P \\Delta W \\right) \\Delta z\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using conjugate of distance to focus plane\n",
    "The telescope objective reimages the beam cross section that is beyond the pupil plane at a distance $l$ before the focal plane. According to Newton's law:\n",
    "\\begin{align*}\n",
    " (\\Delta z + f)l &= f^2\\\\\n",
    " \\Delta z &= \\frac{f(f-l)}{l}\\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### General expression of signal\n",
    "\\begin{align*}\n",
    "  S = \\frac{f(f-l)}{l} \\left( \\frac{\\partial W}{\\partial \\eta} \\delta_c - P \\Delta W \\right)\n",
    "\\end{align*}\n",
    "\n",
    "This last equation shows that the sensor signal consists of two terms.\n",
    "The first term is proportional to the wave-front radial slope at the pupil edge and is localized at the beam edge. The second term maps the wave-front Laplacian across the beam.\n",
    "Since these two terms do not overlap, one can measure them separately and reconstruct the wave-front surface by solving a Poisson equation, using the wave-front derivative normal to the edge as a Neumann-type boundary condition.\n",
    "However, this equation is only a first-order approximation valid for small $\\Delta z$ values,\n",
    "that is, highly defocused stellar images. The algorithm described in the following section uses the solution of the Poisson equation as a first-order solution that is further refined in an iterative process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A basic introduction to Laplace and Poisson partial differential equations\n",
    "This short paragraph is based on [Poisson equation wikipedia page](https://en.wikipedia.org/wiki/Poisson%27s_equation) and [this youtube video](https://www.youtube.com/watch?v=lsY7zYaezto)\n",
    "\n",
    "At the core, a Poisson equation, sought after a scalar valued fonction, and the equation gives you an expression for the Laplacian of this function, which is also a scalar field, example\n",
    "\n",
    "\\begin{align*}\n",
    "  \\Delta f = \\phi\n",
    "\\end{align*}\n",
    "\n",
    "#### Laplace equation / Elliptic PDE\n",
    "\\begin{align*}\n",
    "  \\Delta f = 0 \\\\\n",
    "  \\Delta_{xy} f = 0 \\\\\n",
    "\\end{align*}\n",
    "\n",
    "If the scalar field $\\phi$ is zero everywhere on your (spatial) domain, then when have a Laplace equation. Solutions to Laplace equation are called harmonic. It is pretty interesting to take a quick look at where the harmonic terms comes from: fourier series are expressed in terms of sum of complex harmonic functions (the $e^{2\\pi jnFt}$). It appears that the laplacian (in one dimension) of this family of function is $-(2\\pi nF)^2 e^{2\\pi jnFt}$ hence, one can say that those are eigenfunctions for the laplacian operator. This will be one of the core reason why Fourier series are useful for Laplace equations.\n",
    "\n",
    "Thanks to the linearity property of Laplace operator, one can also see that, the set of solutions of laplace equation (without bounday conditions) is equivalent to the nullspace of $\\Delta_{xy}$.\n",
    "\n",
    "#### Poisson equation / Elliptic PDE\n",
    "\\begin{align*}\n",
    "  \\Delta f = \\phi \\\\\n",
    "  \\Delta_{xy} f = \\phi \\\\\n",
    "\\end{align*}\n",
    "\n",
    "If the scalar field $\\phi$ is not zero everywhere on your (spatial) domain, then when have a Poisson equation. Interestingly, thanks to the linearity property of the Laplacian, one can see that (if we abstract boundary conditions) a set of solution for the Poisson equation can be obtained by finding a particular solution to $\\Delta f = \\phi$ and adding any element from the nullspace of the Laplacian operator. Which is why the two equation are strongly related.\n",
    "\n",
    "#### Heat equation / parabolic PDE\n",
    "\\begin{align*}\n",
    "  \\Delta f = \\phi \\\\\n",
    "  \\Delta_{xy} f = \\frac{\\partial f}{\\partial t} \\\\\n",
    "\\end{align*}\n",
    "\n",
    "In case there's another variable and gradient operator involved in $\\phi$, we are given with heat equation, for which a very famous solution was found by Joseph Fourier, we'll talk about that shortly after. Basically, the evolution of your scalar function $f$ defined over some (spatial) domain $x,y$ with respect to some variable $t$, is given by the laplacian of that very same scalar function itself. This allows so simulate the system by small time steps.\n",
    "\n",
    "#### Wave equation  / hyperbolic PDE\n",
    "\\begin{align*}\n",
    "  \\Delta f = \\phi \\\\\n",
    "  \\Delta f = \\frac{\\partial^2 f}{\\partial t^2} f \\\\\n",
    "\\end{align*}\n",
    "Wave equation is a bit like heat equation, but this time, it is the acceleration of the system (the scalar field) that is described by its own laplacian.\n",
    "\n",
    "#### A bit of math\n",
    "One can see, that, whatever solution you find for your Laplace equation, it will need to be expressed in terms of the given solution, that satisfies your boundary conditions plus the nullspace of the laplacian operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neural PDE\n",
    "\n",
    "### References\n",
    "* This whole section has been designed based on [an article](https://towardsdatascience.com/solving-differential-equations-with-neural-networks-afdcf7b8bcc4) from [Mario Dagrada](https://twitter.com/MarioDagrada)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import IFrame\n",
    "IFrame(\"doc/NeuralPDE/MarioDagradaNeuralPDE.pdf\", width=1200, height=800)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Also [this](https://docs.sciml.ai/NeuralPDE/stable/) amazing website/github project that I also forked. You will also find a basic introduction to neural PDE, and why normalizing flows are useful for those type of models in the NormalizingFlows notebook.\n",
    "\n",
    "* Nvidia also release modulus, an opensource framework for physics informed PDEs introduced [here](https://developer.nvidia.com/blog/develop-physics-informed-machine-learning-models-with-graph-neural-networks/) and available [here](https://catalog.ngc.nvidia.com/orgs/nvidia/teams/modulus/containers/modulus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}