# Microscopy Group

## Improving the Resolution of Microscopy Images

### Mentor: Bogdan Toader
### Participants
- Anvarbek Atayev, University of Warwick
- Nayef Shkeir, University of Warick
- Caelen Feller, Trinity College Dublin
- Daniel Mckinnell
- David Dalton
- Francesco Viganò, University College London

[Overleaf Link](https://www.overleaf.com/project/6058abf5bad63b5b1fac9b6e)

### Monday Summary

- Examined theory of problem
- Divised a provisional algorithm/workflow to use, with description of possible sources of error and limitations.
- Divided into groups to:
- Discretize/implement the theoretical form of the PFS
- Explore the deconvolution solution using proxylab and data loading
- Organized for tuesday (github, these, etc)

### Plan for Tuesday

Necessary installations/tasks for Monday night:
- [Proxylab](https://gitlab.com/jboulanger/proxylab) (download with git or directly, add by using addpath('your directory') in MATLAB terminal)
- Matlab + Image processing toolbox
- Git
- Recommend looking at proxylab implementation and trying to implement some test problems yourself this evening as preperation for tomorrow/to make sure it's all set up!

Tasks:

- Deconvolution inverse problem solution implementation with given PFS/OFS
- Any more theory work we need to do?
- Discretization scheme here? What space are we working in?
- Try using proxylab to implement problem
- What algorithm is best?
- other things that we should consider
- $h_{\text{data}}$ computation
- Compute corrected bead source image using $h_{\text{theory}}$, deconv implementation and a sphere finding algorithm/hand manipulation
- Solve the inverse problem where we have known source and image, but unknown h! (This needs to be reviewed from a theory perspective)
- Apply this to the full image
-
- Need metrics for seeing how well the method does
- Need to explain the possible
---
- Aim 1 Tasks:
- Figure out how to approxipate the pupil fuions with Zerkine polynomials
- Implement and test
- Aim 2 Tasks:
- Check proxylab code and determine how to generate test data
- Perform deconvolution on test data (using what h?)
---

# Resources (slideshow and links)
- Proxylab for Matlab inverse problems (toolbox) [Link](https://gitlab.com/jboulanger/proxylab)
- [A. Chambolle, T. Pock, An introduction to continuous optimization for
imaging](https://hal.archives-ouvertes.fr/hal-01346507/document)
- A. Stokseth. Properties of a Defocused Optical System (Paper, pdf in chat)
- [Wikipedia page for Richardson-Lucy algorithm](https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution)
- Writeboard Fox [Link](https://r8.whiteboardfox.com/8715944-0967-7262)
- Richardson-Lucy Deconvolution [Link](https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution)
- Overleaf [Link](https://www.overleaf.com/project/6058abf5bad63b5b1fac9b6e)
- Point Spread Function [Wikipedia](https://en.wikipedia.org/wiki/Point_spread_function)

# Introduction

Given an object of interest (be it beads/protein) that we want to image, depicted by $u$, its image through a microscope is given by $f$. Ideally, the image recorded through the microscope would equal the source, i.e. $f = u$. But, microscopes aren't perfect, and can record blurry images due to properties of the microscope used, i.e. the image received would be of the form
$$
\begin{align}
f[n_1,n_2,n_3] &= (h \circledast u)[n_1,n_2,n_3] \\
&:= \sum_{m_1=-\infty}^{\infty}\sum_{m_2=-\infty}^{\infty}\sum_{m_3=-\infty}^{\infty}h[m_1,m_2,m_3]u[n_1-m_1,n_2-m_2,n_3-m_3]
\end{align}
$$
where $\circledast$ denotes a 3D discrete convolution, and $h$ is defined as the **point spread function** of a microscope (this describes the various ways the microscope distorts/blurs any point it records of an image).

We assume that $u[n_1,n_2,n_3] = 0$ for $n_i < 0, n_i > 256$, for $i = 1,2$ and $n_3 < 0$ and $n > 32$.


- **Domains:**
- $f, u \in \mathbb{R}^{N_1\times N_2 \times N_3}$, for some $N_i \in \mathbb{N}$. For our purposes, $N_1 = N_2 = 256$ and $N_3 = 32$, i.e. number of slices of image (Note: number of slices of beads image is different)
- what is $h$ in? ($h \in \mathbb{R}^{(2N_1 + 1)\times (2N_2 + 1) \times (2N_3+1)}$?)

#### Introduction of Gaussian Error
In reality, our image data is not only distorted by the properties of the microscope, but also by some noise, which we'll consider to be Gaussian, i.e. in reality, the data that we receive will be
$$
f = h \circledast u + \eta,
$$
where $\eta$ is Gaussian with mean $\mu$ and variance $\sigma^2$.

### Ideal Situation
Ideally, $f = u$, for which one solution is $h[j_1,j_2,j_3] = \delta_{0j_1}\delta_{0j_2}\delta_{0j_3}$, where $\delta_{ij} = 1$ if $i = j$ and $0$ if $i \neq j$, and $\eta \equiv 0$, in which case for all $n_1,n_2,n_3$
$$\begin{align}
f[n_1,n_2,n_3] &= (h \circledast u) [n_1,n_2,n_3] \\
&= \sum_{m_1=-\infty}^{\infty}\sum_{m_2=-\infty}^{\infty}\sum_{m_3=-\infty}^{\infty}h[m_1,m_2,m_3]u[n_1-m_1,n_2-m_2,n_3-m_3] \\
&= \sum_{m_1=-\infty}^{\infty}\sum_{m_2=-\infty}^{\infty}\sum_{m_3=-\infty}^{\infty}\delta_{0m_1}\delta_{0m_2}\delta_{0m_3}u[n_1-m_1,n_2-m_2,n_3-m_3] \\
&= u[n_1,n_2,n_3]
\end{align}$$

### Computing $f$
Note that, since we have a convolution, in Fourier space this transform into a multiplication, i.e.
$$
\hat{f} = \widehat{h \circledast u} = \hat{h}\hat{u}
$$
which is a much more efficient process (numerically speaking).

Thus, when computing $f$ from an inputted $u$ and a known $h$, speedwise (since matrix multiplication + IDFT is faster than convolution), it is beneficial to convert to Fourier space. That is, given known $u$ and $h$, to compute $f$, first compute $\hat{h} \hat{u}$, and then set $f = \text{IDFT}(\hat{h} \hat{u})$, where $\text{IDFT}$ is the [Inverse Discrete Fourier Transform](https://en.wikipedia.org/wiki/Discrete_Fourier_transform).

Note that $\hat{h}$ is the deemed the [Optical Transfer Function (OFT)](https://en.wikipedia.org/wiki/Optical_transfer_function).

---
# Aims
The aim of this project is to deblur the img-gfp-tubuline_aber.tif file, i.e. by solving for $u$, via deconvolution, the equation
$$
f = h \circledast u + \eta.
$$
Although, we should note that we do not yet know what PSF, $h$, is and this must be determined using a calibration file, called beads.tif, which contains an image file of a number of beads of known size. Therefore, the aims of this project can be split into two:

### Aim 1: Calibration
Determine an approximation, $h_{approx}$, of the true value of the PSF $h_{true}$ using $f_{beads}$ (i.e. beads.tif), where we know what $u_{beads}$ should look like.

### Aim 2: Deblurring
Determine $u_{tubuline}$ from blurred image $f_{tubuline}$ (i.e. img-gpf-tubuline_aber.tif) using determined $h_{approx}$ from Aim 1 via a deconvolution technique.

---
# 1. Theory

## 1.1 Point Spread Function

As discussed in the introduciton, whenever an image is recorded using a microscope, the recorded image is not perfect/is not exact and contains abberations and blurring that form due to the structure of the microscope. The abberations introduced by the microscope can be described a a function called the **point spread function** (PSF).

Generally speaking, what is a PSF [point spread function](https://en.wikipedia.org/wiki/Point_spread_function)?


- Also known as impulse response of the optical system
- Represents the response to a single "point object" of the lens
- The spatial domain version of the Optical Transfer Function (OTF, defined on Fourier space as the FT of the PSF).
- PSF represents the effect of the lens on objects (it's a convolution operator), so reversing it's action can be seen as a "known" operator deconvolution problem.
- However, for more accuracy, you need to consider more unpredictable aberrations in the data (sample container shape light influence, and refractive index changes)

### Optical Transfer Function
The optical transfer function (OTF) is defined as the Fourier transform of the PSF. As a Fourier transform, it is generally complex-valued, although, if the PSF is radially symmetric, then the OTF is real valued.
- Specifies how different spatial frequencies are handled by microscope.

### Theoretical PSF Models

A simple theoretical PSF model is given by the following equation

$$h(x,y,z) = \left\vert\int\int p(\kappa_x,\kappa_y)e^{2i\pi z \sqrt{(n/\lambda)^2 - \kappa_x^2 - \kappa_y^2}} e^{2i\pi(\kappa_x x+\kappa_y y)}\mathrm{d}\kappa_x\mathrm{d}\kappa_y\right\vert^2$$

where
- $\lambda$ is the wavelength (here is 0.525 $\mu\mathrm{m}$),
- $n$ is the refractive index (typically ~1.35 here)
- $p(\kappa_x, \kappa_y)$ is called the **pupil function**

The theoretical PSF is useful for providing

- A baseline with which to test our deblurring algorithm without accounting for aberrations (just use a non-aberration accounting one to get some initial metrics/results)
- A part of a more complex theoretical model which includes some of the abberations we can observe.

### Pupil Function
*Pupil Function* - "describes how a light wave is affected upon transmission through an optical imaging system" [Wikipedia](https://en.wikipedia.org/wiki/Pupil_function)

- Complex function of position of microscope aperture indicating the relative change in amplitude and phase of the light wave. Or just indicates if light is being transmitted if less general.
- Can be seen as a function which constrains the integral to a (unit) disc here.
- The pupil function captures all optimal aberrations that occur between the image plane and the focal place in the sample.

In general, the pupil function can be defined as
$$
\begin{align}
p(k_1,k_2) = A(k_1,k_2) e^{i\Theta(k_1,k_2)}, \quad (k_1,k_2) \in \mathbb{R}^2,
\end{align}
$$ where $\Theta$ is the *phase change* introduced by the optics/surrounding medium, and $A$ denotes the amplitude of light.

#### Example: Ideal System
In an ideal system, the pupil function can be seen as being equal to one at every point within the pupil (aperture of camera) and zero outside it. In case of a circular pupil, the pupil function can be defined as
$$
p(k_1,k_2) :=
\begin{cases}
1 & \text{if } k_1^2 + k_2^2 \leq \frac{\text{NA}}{\lambda}, \\
0 & \text{otherwise},
\end{cases}
$$
where
- $\text{NA}$ is the [numerical aperture](https://en.wikipedia.org/wiki/Numerical_aperture)
- $\lambda$ is the wavelength of the light used
- $\sqrt{\frac{\text{NA}}{\lambda}}$ is the pupil radius

#### Example: Non-Ideal System

## 1.2 Discretisation of theoretical PSF
We choose to discritise the PSF using a standard Riemann sum approximation, although higher order approximations can be used and will be considered if necessary. In particular, for a discretisation fineness parameter $N$, denote the discritation approximation of the PSF by
$$h_N(x,y,z) := \left\vert\sum_{i=-N}^{N-1}\sum_{j=-N}^{N-1} (\Delta_{k_i}\Delta_{\omega_j})p(k_i,\omega_j)e^{2i\pi z \sqrt{(n/\lambda)^2 - k_i^2 - \omega_j^2}} e^{2i\pi(k_i x+ \omega_j y)}\right\vert$$
where
- $k_i = \frac{i}{N}\frac{1}{\lambda}$, thus $\Delta_{k_i} = \frac{1}{N\lambda}$
- $\omega_i = \frac{i}{N}\frac{1}{\lambda}$, thus $\Delta_{\omega_i} = \frac{1}{N\lambda}$

and $h_N \to h$ pointwise as $N\to\infty$.

#### Data-Driven PSF Models

- Purely data-driven: Estimate the PSF deconvolution matrix from a comparable calibration image of spherical beads. Points to consider:

- Are the aberrations space-invariant? This determines the size of the PSF? (We can estimate it for just a single segmented point source and average or break the overall image into segmented regions which )

-

For a single lens, an on-axis point source will produce an [Airy disk](https://en.wikipedia.org/wiki/Airy_disc), an off-axis one producing a warped version of this?

# 2. Aim 1: Calibrating
Recall that in this section we are aiming to determine an approximation, $h_{approx}$, of the true value of the PSF $h_{true}$ using known $f_{beads}$ (i.e. beads.tif).

## 2.1 Zerkine Polynomials
Zerkine polynomials are a sequence of polynomials that are orthogonal on the unit disk ([Wikipedia](https://en.wikipedia.org/wiki/Zernike_polynomials))

There are even and odd Zerkine polynomials. Respectively, the even and odd Zerkine polynomials are defined as
$$
Z_{n}^{m}(\rho,\varphi) = R_n^m(\rho)\cos(m\varphi) \qquad Z_{n}^{-m}(\rho,\varphi) = R_n^m(\rho)\sin(m\varphi),
$$
where $R_n^m$ are radial polynomials defined as
$$
R_n^m(\rho) :=
\begin{cases}
\sum_{k=0}^{\frac{n-m}{2}}\frac{(-1)^k(n-k)!}{k!\left(\frac{n+m}{2}-k\right)!\left(\frac{n-m}{2}-k\right)!}\rho^{n-2k}, & n-m \text{ is even} \\
0 & \text{otherwise}
\end{cases}\quad \rho \in [0,1].
$$
#### Zerkine Transform
Any sufficiently smooth real-valued phase field over the unit disk $G(\rho,\varphi)$ can be represented in terms of its Zerkine coefficients (odd and even), i.e.
$$
G(\rho, \varphi) = \sum_{m,n}[a_{m,n}Z_n^m(\rho,\varphi) + b_{m,n}Z_n^{-m}(\rho,\varphi)],
$$
where the coefficients $a_{m,n}$ and $b_{m,n}$ are defined as
$$
\begin{align}a_{{m,n}}&={\frac {2n+2}{\epsilon _{m}\pi }}\left\langle G(\rho ,\varphi ),Z_{n}^{{m}}(\rho ,\varphi )\right\rangle ,\\b_{{m,n}}&={\frac {2n+2}{\epsilon _{m}\pi }}\left\langle G(\rho ,\varphi ),Z_{n}^{{-m}}(\rho ,\varphi )\right\rangle, \end{align}
$$
where
$$
\epsilon_m := \frac{1}{\pi}\int_{0}^{2\pi}\cos(m\varphi)^2 d\varphi.
$$
$\epsilon_m$ is defined as being 2 if $m=0$ and 1 if $m \neq 0$.

#### Radial Zerkine Transform
If our desired function $G$ is radial, i.e. independent of $\varphi$, then
$$
G(\rho) = \sum_{n=1}^{\infty}a_{0,n}Z_{n}^0(\rho,0).
$$
For brevity, write $a_{0,n}$ as $a_n$ and $Z_n^0(\rho,0)$ as $Z_n(\rho)$, which yields
$$
G(\rho) = \sum_{n=0}^{\infty}a_{n}Z_{n}(\rho).
$$
Recalling that $Z_{n} \equiv 0$ for $n$ odd, we have that for a radial function $G(\rho)$
$$
G(\rho) = \sum_{n=0}^{\infty}a_{2n}Z_{2n}(\rho), \quad \rho \in [0,1].
$$

#### Zerkine Polynomials of Interest: Radial
For an approximation of a radial $G$ via Zerkine polynomials up to second order, the *Zernike polynomials* that we will need are
$$
\begin{align}
Z_0(\rho) &= Z_0^0(\rho,0) = R_{0}^{0}(\rho) = 1 \\
Z_2(\rho) &= R_2^0(\rho) = 2\rho ^{2}-1
\end{align}
$$
with coefficients
$$
\begin{align}
a_0 &= a_{0,0} = \frac{1}{\pi}\langle G(\rho), Z_0(\rho)\rangle \\
a_2 &= a_{0,2} = \frac{3}{\pi}\langle G(\rho), Z_2(\rho)\rangle
\end{align}
$$
where $\langle F,G\rangle$ is defined as:
$$
\begin{align}
\langle F,G \rangle &= \int_0^1 \int_0^{2\pi} F(\rho, \varphi) G(\rho, \varphi) \rho d\rho d\varphi
\end{align}
$$

and thus we can write
$$
\begin{align}
G(\rho) = a_0Z_0(\rho) + a_2Z_2(\rho)
\end{align}
$$

#### Zerkine Polynomials of Interest: Non-Radial
For an approximation up to second order, we have that
$$
G(\rho,\varphi) \approx \sum_{m=0,n=0}^{2,2}[a_{m,n}Z_n^m(\rho,\varphi) + b_{m,n}Z_n^{-m}(\rho,\varphi)],
$$
where
$$
\begin{align}
R_1^1(\rho) &= \rho \\
R_2^2(\rho) &= \rho^2
\end{align}
$$
which means that
$$
\begin{align}
Z_1^{-1}(\rho,\varphi) &= 2\rho \sin(\varphi) \\
Z_1^1(\rho,\varphi) &= 2\rho \cos(\varphi) \\
Z_2^{-2}(\rho,\varphi) &= \sqrt{6} \rho \sin(2 \varphi) \\
Z_2^2(\rho,\varphi) &= \sqrt{6} \rho \cos(2 \varphi)
\end{align}
$$
with coefficients
$$
\begin{align}
a_{1,1} &= \frac{4}{\pi}\langle G(\rho,\varphi), Z_1^1(\rho,\varphi)\rangle \\
a_{2,2} &= \frac{6}{\pi}\langle G(\rho,\varphi), Z_2^2(\rho,\varphi)\rangle \\
a_{2,0} &= \frac{1}{\pi}\langle G(\rho,\varphi), Z_0^2(\rho,\varphi)\rangle\\
b_{0,0} &= a_{0,0} \\
b_{1,1} &= \frac{4}{\pi}\langle G(\rho,\varphi), Z_1^{-1}(\rho,\varphi)\rangle \\
b_{0,2} &= a_{0,2} \\
b_{2,0} &= \frac{1}{\pi}\langle G(\rho,\varphi), Z_0^{-2}(\rho,\varphi)\rangle\\
b_{2,2} &= \frac{6}{\pi}\langle G(\rho,\varphi), Z_2^{-2}(\rho,\varphi)\rangle
\end{align}
$$

## 2.2 Method
We shall use Zerkine polynomials to approximate an appropriate PSF. In particular, we shall use Zerkine polynomials to approximate the phase part of the *pupil function*, given by $\Theta$.

Let (Z_n^m) denote the sequence of Zerkine polynomial, then assuming that $\Theta$ is smooth enough, we can write $\Theta$ as (after converting into polar coordinates)
$$
\begin{align}
\Theta(k_1,k_2) \equiv \Theta(\rho,\varphi) = \sum_{n=0}^{\infty}\sum_{m=0}^{\infty}[a_{m,n}Z_n^m(\rho,\varphi) + b_{m,n}Z_n^{-m}(\rho,\varphi)],
\end{align}
$$
where
$$
\begin{align}
k_1 = \rho\cos(\varphi), \quad k_2 = \rho\sin(\varphi).
\end{align}
$$
and $(a_{m,n})$ and $(b_{m,n})$ are coefficients defined as in Section 2.1. The above is an infinte series, although we shall approximate this up to second order, and write
$$
\begin{align}
\Theta(\rho,\varphi) \approx \Theta_{approx}(\rho,\varphi) &:=
\sum_{n=0}^{2}\sum_{m=0}^{2}[a_{m,n}Z_n^m(\rho,\varphi) + b_{m,n}Z_n^{-m}(\rho,\varphi)] \\
&= 2a_{0,0}Z_0^0(\rho,\varphi) + 2a_{0,2}Z_2^0(\rho,\varphi) + a_{2,2}Z_{2}^{2}(\rho,\varphi) + a_{1,1}Z_{1}^{1}(\rho,\varphi) \\
&\qquad + b_{2,2}Z_{2}^{-2}(\rho,\varphi) + b_{1,1}Z_{1}^{-1}(\rho,\varphi)
\end{align}
$$
Given the above expansion of the phase part of the pupil function, we may then approximate the pupil function as (after assuming $A \equiv 1$)
$$
p(\rho,\theta) \approx p_{approx}(\rho,\theta) := \exp(i\Theta_{approx}(\rho,\theta)),
$$
and thus, we can approximate the PSF as
$$
\begin{align}
h(x,y,z) \approx h_{approx}(x,y,z) := \left\vert\int\int p_{approx}(\kappa_x,\kappa_y)e^{2i\pi z \sqrt{(n/\lambda)^2 - \kappa_x^2 - \kappa_y^2}} e^{2i\pi(\kappa_x x+\kappa_y y)}\mathrm{d}\kappa_x\mathrm{d}\kappa_y\right\vert
\end{align}
$$

##### Method
- To determine an approximation for the true $h_{true}$, we perform an inversion on the beads file to determine $u_{theory}$ with the Zerkine expanded pupil function.
- We know the size of the beads and we know that they are spherical, thus it suffices to choose the coefficients of the Zerkine polynomial expansion that result in a best fit.

### Determining Optimal Coefficients
Our aim is to find the coefficients $\mathbf{x} := (a_{0,0},a_{0,2},a_{2,2},a_{1,1}, b_{1,1},b_{2,2})$
such that $h_{approx}$ well approximates $h_{true}$.

### The inverse problem
Given a sample $u$ , a PSF $h$ and an image $f$, we want to solve the following minimisation problem to give us the deconvoluted (deblurred) image.
$$
\begin{align}
\min_{u}({ \|u\|_{TV} + \frac{\lambda}{2} \|f - h \circledast u \|_{2}^{2}})
\end{align}
$$
where
$$
\begin{align}
\| x \|_{TV} = \sum_{l=1}^{N}(|x^{(l)} - x^{({n_{l,1}})}|^2 + |x^{(l)} - x^{({n_{l,2}})}|^2)^{1/2}
\end{align}
$$
where $(n_{l,1},n_{l,2}) \in \{1,\cdots,N\}^2$ denote the positions of the horizontal/vertical nearest neighbours of $x^{(l)}$.
To solve this problem, we apply the Primal-Dual Hybrid Gradient (PDHG) algorithm from the proxylab toolbox:


This algorithm solves problems of the form
$$
\begin{align}
\min_{x \in X}(f(Kx) + g(x))
\end{align}
$$

where $f,g$ are convex and $K : X \to Y$ is a bounded linear operator. We can then write the ptoblem as
$$
\begin{align}
\max_{y}\inf_{x} \langle y,Kx\rangle - f^* (y) + g(x)
\end{align}
$$
where we can alternate a descent in the variable $x$ and an ascent in the dual variable $y$
$$
\begin{align}
x^{k+1} &= prox_{\tau g} (x^k - \tau K^* y^k)\\
y^{k+1} &= prox_{\sigma f^* } (y^k - \sigma K x^{k+1}).
\end{align}
$$

# Questions
- Is $h$ lossless? (i.e. does it result in loss of information from source?)
- Not in an ideal case, as it is a linear operator. There may be loss of information due to discretisation error though, and the deconvolution is an ill-posed problem.

# Definitions
- *Optical aberration*
In optics, aberration is a property of optical systems such as lenses that causes light to be spread out over some region of space rather than focused to a point. Aberrations cause the image formed by a lens to be blurred or distorted.

# Rapid Talk Script

## Introduction
- Imaging an object with a microscope changes it due to effects from the lens, e.g. image can be blurred
- Mathematically, this blurring effect of a microscope can be described by a point spread function, here called $h$
- If an image can be described by a matrix $u$, then its recorded image $f$ can be assumed to be of the form
$f = h \circledast u + \eta$, where $\eta$ represents noise within the system
- In this project we are presented with a 3D microscope image of tubulin, which we aim to deblur using deconvolution techniques (i.e. obtain $u$)

## Current Progress
- To achieve our aim, we've split the task into two steps: calibrate and then deblur.

### Step 1:
- We do not currenlty know what the PSF of the microscope is, and so must obtain it via an approximation.
- We have approximated the PSF by expanding it's pupil function as a series of Zerkine polynomials and solve for coefficients that result in a best match with a calibration image file (provided by Bogdan)
- We are currently in the process of implementing this numerically

### Step 2: Deconvolution
- The next step is to perform the deblurring proceduce.
- We have implemented a deconvolution operation on a simulated dataset where the true form of the PSF is known.
- Below we show a test image that has been sharpened wih the deconvolution algorithm, the image from the microscope, the deconvoluted sharpened image and the true image that we are aiming for.