# Challenge 5
## Fundamental Limits of a Single Photon Detector

# Group Members:

Nick Polydorides: [e-mail]
Oliver Dorn: [e-mail]
Guiyun Tian: [e-mail]
Richard Claridge (PA): [e-mail]
Milo Edwards (DSTL): [e-mail]
Emily Russell (DSTL): [e-mail]
Sean Holman: [e-mail]
Kurt Langfeld: [e-mail]

# Slides for the Friday Presentation

<!-- Kurt: these slides are hust a suggestion; feel free to move around or insert stuff between the slides
-->

-------------------------
-------------------------
-------------------------
-------------------------
-------------------------
-------------------------
-------------------------
-------------------------
-------------------------
-------------------------
-------------------------
-------------------------

-------------------------

<!--
**OVERLEAF TEMPLATE https://pr~~~~otect-eu.mimecast.com/s/-DzEC48Wpc9nvYgtOrCgY?domain=scanmail.trustwave.com - please edit challenge5.tex**

Initial members of the Challenge group:
David Abrahams, Emily Russell (dstl), GuiYun Tian, Kurt Langfeld, Nick Polydorides, Sean Holman, Luca Sapienza, Richard Claridge (PA), Tom Border (dstl), Milo Edwards (dstl)

Possible helpful articles:
* [Overview of single photon detection technologies](https://ieeexplore.ieee.org/document/6358468)
* [Photon Detection Characteristics and Error Performance of SPAD Array Optical Receivers](https://www.research.ed.ac.uk/portal/files/21842579/ssh1501.pdf)
* [Three-dimensional single-photon imaging through obscurants](https://www.osapublishing.org/oe/abstract.cfm?URI=oe-27-4-4590)
* [Laser profiling for airborne target classification](https://www.spiedigitallibrary.org/conference-proceedings-of-spie/10636/1063602/Laser-profiling-for-airborne-target-classification/10.1117/12.2303965.short)
* [Real-time 3D reconstruction from single-photon lidar data using plug-and-play point cloud denoisers](https://www.nature.com/articles/s41467-019-12943-7)
* [Single-photon computational 3D imaging at 45 km](https://arxiv.org/pdf/1904.10341.pdf)
* [Resolution-enhanced quantum imaging by centroid estimation of biphotons](
https://www.osapublishing.org/optica/abstract.cfm?uri=optica-6-3-347)

Additional papers of interest referenced at bottom.


Scenario: the laser (L) emits photon packets (10,000s photons) which experience absorption through the atmosphere and an obscurant (O e.g. fog, mist, snow, dust, foliage). Some are transmitted through to the target (T). Some photons (the photons of interest) are reflected and return to the detector which can also record their time of arrival, wavelength and intensity. Some photons are scattered and do not return to the detector, some photons are scattered and DO return to the detector. The detector also experiences noise from background photons (primarily Sun) and dark current. The detector can be gated such that photons that arrive d/c seconds after emission are prioritised.

* Can polarisation help as filter on top of the time-of-flight?
* The change of polarisation by scattering from metal would be known. Can this be used to identify metalic parts of the target?
* What information could we extract from phase or directionality (e.g. a plenoptic variant)
* If the target is static, can we reconstruct the target information from a (large) background of noise from moving obstructions?
* How does multiple sampling help (e.g. moving target/moving source)
* Feature size limit
* Impact of atmospheric turbulence
* Is it enough to consider a single detector, or do we need to think about the array all together? - how does this differ to a standard pixel array beyond timing information?
* Do we need to model multiple photon scattering events to describe the measurements?
* How do we parameterise the "target" so that we can say that we have resolved this from low photon counts? How do we parameterise the domain? (absorbers+scatterers?)
* As the interim medium becomes "thicker" the measurements will drop to zero.
* Confirm we define limit as feature size as distinct from rate, range, illumination required
* How sensitive are the detectors to direction? Can we assume detected photons come from a certain direction?
* What is a reasonable limit on rate - for detector and illuminator (repitition rate and pulse length)
* In atmospheric lidar imaging (DIAL), they use 2 wavelengths, one that's sensitive to absorption of the medium and another that's not. They subtract and that gives them a signal that's not sensitive to the scatter. Wish there was something similar for this.

The supplementary doc of the Nature paper is quite informative - denoising Poisson data.

Here is the SIAM paper explaining the method in the nature paper
https://epubs.siam.org/doi/pdf/10.1137/18M1183972

RC: What can we learn from MRI/CT -> we have a series of 2D slices in time from which we want to create an image (or consider voxel-wise)

Additional papers:

Heriot-Watt:
* [Real-time 3D reconstruction from single-photon lidar data using plug-and-play point cloud denoisers](
https://www.nature.com/articles/s41467-019-12943-7)
* [Long-range depth imaging using a single-photon detector array and non-local data fusion](
https://www.nature.com/articles/s41598-019-44316-x)
* [Three-dimensional single-photon imaging through obscurants](https://www.osapublishing.org/oe/abstract.cfm?URI=oe-27-4-4590)
* [Single-photon three-dimensional imaging at up to 10 kilometers range](https://www.osapublishing.org/oe/abstract.cfm?&uri=oe-25-10-11919)
* [3D LIDAR imaging using Ge-on-Si single–photon avalanche diode detectors](https://www.osapublishing.org/oe/fulltext.cfm?uri=oe-28-2-1330&id=425784)
* [Resolving range ambiguity in a photon counting depth imager operating at kilometer distances](https://www.osapublishing.org/oe/fulltext.cfm?uri=oe-18-9-9192&id=198258)

Sweden (FOI):
* [Laser profiling for airborne target classification](https://www.spiedigitallibrary.org/conference-proceedings-of-spie/10636/1063602/Laser-profiling-for-airborne-target-classification/10.1117/12.2303965.short)
* [Imaging and laser profiling for airborne target classification](https://www.spiedigitallibrary.org/conference-proceedings-of-spie/10434/104340K/Imaging-and-laser-profiling-for-airborne-target-classification/10.1117/12.2277531.short)
* [Past and present laser sensing activities at the Swedish Defense Research Agency (FOI)](https://www.osapublishing.org/abstract.cfm?uri=lsc-2019-LM2B.1&origin=search)
* [Photon-counting panoramic three-dimensional imaging using a Geiger-mode avalanche photodiode array](https://www.spiedigitallibrary.org/journals/optical-engineering/volume-57/issue-9/093104/Photon-counting-panoramic-three-dimensional-imaging-using-a-Geiger-mode/10.1117/1.OE.57.9.093104.short)
* [Panoramic single‑photon counting 3D lidar](https://www.spiedigitallibrary.org/conference-proceedings-of-spie/10796/2325691/Panoramic-singlephoton-counting-3D-lidar/10.1117/12.2325691.short)

Single photon tech overivews:
* [Single photon detection system for visible and infrared spectrum range](https://www.osapublishing.org/ol/abstract.cfm?uri=ol-43-24-6085)
* [Single Photon Detection and Timing: Experiments and Techniques](https://www.sciencedirect.com/science/article/pii/S006525390860222X)

Potential to monte-carlo the effect of polarisation - ray tracing approach
Nick has looked at work from Herriott-Watt. Problem is a de-noising problem. Formulate the observation as a poisson process; measure a cloud of points x,y,time. Look for correlations that could be sitting on the same surface. Apply inhibition & form surfaces.
Inverse problem for radiative transform equation?
What can we lift from optical comms - poisson channel modelling
How can we exploit asymmetries

Potential Directions:
* Polarisation >> Kurt
* 2-wavelength approach (requires prior information)
* Moving target/moving noise
* Reverse emitter and reciever (reciprocity)
* Model as a wireless optical communication channel (bypass the physical characterisation of the "uncertain" channel - here the obsructing medium)
* Can the Doppler effect been used to detect moving targets (combination with time-of-flight filters)?

==11:20AM===
* Breakout to work up ideas seperately
* Reconvene at 15:30 to share work
* Produce 'working slides' to eventually form final presentation
* Question: How big is the surface of a single SPAD? Can I assume it to be a point?
* A: of order 5um per pixel. Probably safe to assume a point in the first instance. Note about 512x512 pixels though
* Question: considerations of speed of target/obscurant relative to frame rate

For discussions (Kurt):

unobscured square in the light of 1M photons.


and with a screen in front of:


I saw Oliver Dorn (Manchester) in the audience - he is an expert on the RTE..

@Kurt: how many scattering events did you take in the above simulation?
Relevant scattering only takes place on the screen; twice: outgoing to the target and a second time on the way back.
This is not a restriction; could simulate 'fog' as well (just takes more computational time).

If it is possible can we have the time-resolved photon counts for this?
Yes, this is the next step...

**Some slides perhaps useful for the pesentation**

One pulse now illuminates with 10M photons
Scattering characetristics is chosen adhoc



*Great. I will use the first 3. I made some hand drawings to show the effect of the scatter screen on the tail, but your simulations are much nicer, so I will add them.
- Did you use an arbitrary speed for light?
- The scattering characteristic of the screen is anisotropic - I wonder if that has an effect on the shape of the image? Surely explains the low returns.

Nice paper on simulation for multiple scattering from computer graphics people
http://vision.seas.harvard.edu/inverse_transient/inverse.pdf

[several days on clusters to compute]

Update on the slides (Kurt):
* more isotropic scattering
* inclduing initial reflectiosn form the screen
* scattering model explained (Heney Greesite could be implemented; not expensive)
* 10M rays take 3secs on a Desktop (Fortran; Matlab takes a factor 10 longer)


-->

# Popular model for radiative transfer
In many applications with scattering the propagation of photons (or particles) can be modeled by the scalar time-dependent linear transport equation
\begin{equation}
\frac{1}{c}\frac{\partial u}{\partial t} + \theta \cdot \nabla u(x,\theta,t) +
a(x)u(x,\theta,t)=b(x) \int_{S^{2}} \eta(\theta \cdot
\theta^{'})u(x,\theta^{'},t)
d\theta^{'} + q(x,\theta,t),
\end{equation}
with appropriate initial and possibly a boundary condition.
Here $S^{2}$ denotes the unit sphere in $R^{3}$ and we assume that
\begin{equation}
\int_{S^{2}}\eta(\theta\cdot\theta^{\prime}) d\theta^{\prime} =1.
\end{equation}
In many biomedical applications the scattering phase function $\eta(\theta\cdot\theta^{\prime})$ is highly forward peaked. A popular model for these applications is the Henyey-Greenstein (HG) scattering function in 3D
\begin{equation}
\eta (\theta\cdot\theta^{\prime} )=\frac{1-g^{2}}{4\pi (1+g^{2}-2g\theta\cdot\theta^{\prime} )^{\frac{3}{2}}}
\,=\,\sum_{n=0}^{\infty}\frac{2n+1}{4\pi}g^{n}P_n(\theta\cdot\theta^{\prime}),
\end{equation}
where $P_n$ is the $n$-th order Legendre polynomial.
The anisotropy factor $-1\leq g\leq 1$ in this function has the meaning of a mean scattering cosine.
$g=0$ indicates isotropic scattering, $g>0$ primarily forward scattering, and $g<0$ primarily backward scattering. In biomedical applications
we have approximately $0.9<g<0.95$ or higher which indicates highly peaked forward scattering.

It is popular to denote further
\begin{equation}
\sigma_{a}(x)=a(x)-b(x)
\end{equation}
which measures the probability at a given location that a particle is absorbed and disappears completely.

# Applications

- Biomedical Imaging (diffuse optical tomography, fluorescence tomography)
- Atmospheric Physics (understanding stellar atmospheres, scintilation in satellite comms)
- Computer Graphics (modeling light propagation in turbid media)
- Seismic imaging (multiple scattering of seismic waves from earthquakes)
- Nuclear Reactor Physics (calculating criticality conditions for reactors)
- Multiple scattering in communication (e.g underwater communication)
- and many more ...

# Modeling

- Monte Carlo Simulations
- Numerical models (Finite Differences, Finite Elements, Spectral Methods, ...)
- Analytical methods based on Green's functions
- Neumann series approaches (Born series, Scattering series)

# Visualization (here Monte Carlo in 2D)
In the following we show some snap shots from Monte Carlo simulations for light propagation in scattering tissue. Dimensions are in cm, but similar behaviour might be expected when modeling scattering of photons propagating in forests.
Notice that the determination of parameters $b$ and $g$ significantly affects the overall behaviour of photon propagation in the turbid medium. Absorption is assumed constant and small in the simulations shown here.
Above figure: $b=0.05\mbox{cm}^{-1}$, $g=0.0$. Here the ballistic (unscattered)
contribution is dominant (bottom left image) and only a small scattered side-lobe appears (bottom right image).
Above figure: $b=0.5\mbox{cm}^{-1}$, $g=0.0$. The scattered part becomes stronger but the ballistic contribution is still visible. A more spread-out wave-front becomes visible and increased isotropic scattering contributions give rise to more significant particle densities behind this wave front.
Above figure: $b=2.5\mbox{cm}^{-1}$, $g=0.9$. Now scattering becomes more dominant with a highly forward peaked scattering phase function. The corresponding photon density distribution resembles now more a damped wave propagating inside $\Omega$ which might be well approximated by a telegraph equation.
Above figure: $b=25.0\mbox{cm}^{-1}$, $g=0.9$. Now scattering dominates the particle propagation but due to the highly forward peaked scattering phase function the penetration depth of the particles is still good. The diffusion approximation seems a good model for this situation.
Above figure: $b=25.0\mbox{cm}^{-1}$, $g=0.0$. This situation is similar to the one shown before, but now, due to the isotropic scattering, the penetration depth of the particles is smaller. Also here the diffusion approximation seems a good model for this situation.
Above figure: $\sigma_a=0.025\mbox{cm}^{-1}$, $b=100.0\mbox{cm}^{-1}$, $g=0.9$. Shown are the time-resolved MC-simulated data at different points along the boundary of the modelled domain. In particular, two time series are highlighted, one close to the source (backscattering) and one at the side. Light is heavily scattered by the environment, but single scattering events are mainly forward peaked.

*Reference:* O. Dorn.
Distributed parameter estimation for the time-dependent radiative transfer equation.
In H. Antil, D Kouri, M Lacasse, and D Ridzal, editors, Frontiers in PDE constrained Optimization, volume 163 of IMA Volumes in Mathematics and its Applications, pages 341-375. Springer, 2018.
DOI: 10.1007/978-1-4939-8636-1_10

# One dimensional model

In one dimension there are only two directions and we will label the flux to the right as $u_R$ and the left as $u_L$. The RTE is then a coupled system of two equations
\begin{equation}
\frac{1}{c} \frac{\partial u_R}{\partial t} + \frac{\partial u_R}{\partial x} + \sigma_{a+s} u_R = \sigma_s u_L,
\end{equation}
\begin{equation}
\frac{1}{c} \frac{\partial u_L}{\partial t} - \frac{\partial u_L}{\partial x} + \sigma_{a+s} u_L = \sigma_s u_R.
\end{equation}
Boundary conditions are ($h$ models the shape of the pulse)
\begin{equation}
u_R(t,0) = h(t),
\end{equation}
\begin{equation}
u_L = 0 \mbox{ for $x \gg tc$}.
\end{equation}
I know this model won't capture spreading due to scattering out of the line, but I wonder if that can be incorporated into a 1D model in some way perhaps assuming rotational symmetry about the direction of the laser.

The solutions $u_R$ and $u_L$ are expanded in scattering series
\begin{equation}
u_R = \sum_{j=0}^\infty u_{R,j}, \quad u_L = \sum_{j=0}^\infty u_{L,j},
\end{equation}
where
\begin{equation}
\frac{1}{c} \frac{\partial u_{R,0}}{\partial t} + \frac{\partial u_{R,0}}{\partial x} + \sigma_{a+s} u_{R,0} = 0, \quad u_{R,0}(t,0) = h(t),
\end{equation}
\begin{equation}
\frac{1}{c} \frac{\partial u_{L,0}}{\partial t} - \frac{\partial u_{L,0}}{\partial x} + \sigma_{a+s} u_{L,0} = 0, \quad u_{L,0} = 0 \mbox{ for $x \gg tc$}.
\end{equation}
which results in
\begin{equation}
u_{R,0}(t,x) = h\left ( t-\frac{x}{c} \right ) e^{-\int_0^x \sigma_{a+s}(x')\ \mathrm{d} x'}, \quad u_{L,0} = 0.
\end{equation}
For $j \geq 1$ we require
\begin{equation}
\frac{1}{c} \frac{\partial u_{R,j}}{\partial t} + \frac{\partial u_{R,j}}{\partial x} + \sigma_{a+s} u_{R,j} = \sigma_s u_{L,j-1}, \quad u_{R,j}(t,0) = 0,
\end{equation}
\begin{equation}
\frac{1}{c} \frac{\partial u_{L,j}}{\partial t} - \frac{\partial u_{L,j}}{\partial x} + \sigma_{a+s} u_{L,j} = \sigma_s u_{L,j-1}, \quad u_{L,j} = 0 \mbox{ for $x \gg tc$}.
\end{equation}
Solving these gives
\begin{equation}
u_{R,j}(t,x) = \int_0^x e^{-\int_{\tilde{x}}^x \sigma_{a+s}(x')\ \mathrm{d} x'} \sigma_s(\tilde{x})\ u_{L,j-1}\left ( t + \frac{\tilde{x}-x}{c},\tilde{x} \right ) \ \mathrm{d} \tilde{x},
\end{equation}
\begin{equation}
u_{L,j}(t,x) = \int_x^\infty e^{-\int_x^{\tilde{x}} \sigma_{a+s}(x') \ \mathrm{d} x'} \sigma_s(\tilde{x}) u_{R,j-1}\left ( t + \frac{x-\tilde{x}}{c},\tilde{x} \right ) \ \mathrm{d} \tilde{x}.
\end{equation}
Setting h to be a delta function we could get an idea of spreading due to forward/backward scattering.

# Differential absorption lidar (DIAL) for gas imaging

A simple example in doing tomography of gases using the RTE to demonstrate the technique of DIAL. This uses data from one fixed in space light source and SPAD-type arrays as the backscatter from gases is tiny. In the conventional sense the inverse problem is massively under-determined, so we regularise by coupling the RTE with the dispertion equation that governs the evolution of the plumes. The problem is similar to the setting of challenge 5 in that we are trying to detect a certain gas within a gas mixture obstructed by fog, and the image needs to be done close to real time. Of course solving the RTE in 3D takes long and so effective reparameterisations of the image space are crucial to reduce the dimension and make the problem better posed. In 3D the time-dependent Radiative Transfer Equation (RTE)
\begin{equation}
\Bigl ( \frac{1}{c}\frac{\partial}{\partial t} + v \cdot \nabla_{x} + \sigma_{\mathsf{on/off}}(x) \Bigr ) H_{\mathsf{on/off}}(x,v,t) = \int_{\mathbb{S}^{2}} H_{\mathsf{on/off}}(x,v',t) k(x, v\cdot v') dv'
\end{equation}
is used in order to describe the propagation of light intensities $H_{\mathsf{on/off}}$ for both wavelengths. It is assumed that the target is visible in the difference of the extinction coefficients $\sigma_{\mathsf{on}}-\sigma_{\mathsf{off}}$. Measurements are assumed to be obtained from a scanning Lidar instrument which measures data
\begin{align*}
\mathbf{M}_{ij}&\sim \mathsf{Poisson}(L_{\mathsf{on}}(t_i,v_j)) \\ \mathbf{N}_{ij}&\sim \mathsf{Poisson}(L_{\mathsf{off}}(t_i,v_j))
\end{align*}
with intensities $L_{\mathsf{on/off}}$ that correspond to angularly averaged solutions $H_{\mathsf{on/off}}$ of the RTE at time $t_i$ and initial pulse direction $v_j$. Reconstruction of differential absorption profiles from a simulated $50 \times 50$ Lidar scan with $50$ time bins, assuming noise corrupted data with single digit photon counts in most bins. A conventional reconstruction at that spatial resolution is not possible from such measurements. Method will be described in our draft paper (arXiv by the end of September).
(i) Ground truth (view in wind direction)
(ii) Ground truth (view from light source)
(iii) Absolute error (view from light source)
(iv) Absolute error (view in wind direction)