%!TEX root = ../dissertation.tex
%\begin{savequote}[75mm]
%Nulla facilisi. In vel sem. Morbi id urna in diam dignissim feugiat. Proin molestie tortor eu velit. Aliquam erat volutpat. Nullam ultrices, diam tempus vulputate egestas, eros pede varius leo.
%\qauthor{Quoteauthor Lastname}
%\end{savequote}

\chapter{Random-Subset Fitting of Digital Holograms for Fast 3D
  Particle Tracking}
\label{chap:random-subset}

\newthought{Digital holographic microscopy} (DHM) is a powerful tool for measuring
three-dimensional dynamics of microscale objects suspended in a
fluid. A holographic microscope records a two-dimensional (2D) image,
or hologram, resulting from interference between a reference light
field and the light scattered from a sample.  The 2D hologram encodes
phase information that can be used to infer the 3D structure of the
sample. Because holograms can be captured with a fast camera, DHM
enables measurements of 3D fluctuations and motion on short
timescales, making it particularly useful for measuring dynamics of
swimming organisms~\cite{lewis_swimming_2006, yeom_real-time_2006,
  sheng_digital_2007, moon_automated_2009}, living
cells~\cite{choi_three-dimensional_2009}, and colloidal
particles~\cite{xu_digital_2002, satake_three-dimensional_2005}.

The principal challenge in DHM is recovering the 3D information from
the 2D hologram. The most common method is numerical reconstruction,
the digital descendant of Gabor's optical
reconstruction~\cite{gabor_new_1948}. Numerical reconstruction
techniques use scalar diffraction theory to back-propagate light from
a measured hologram, yielding a 3D image of the
sample~\cite{schnars_digital_2002, fugal_practical_2009}.  However,
the 3D image can contain artifacts if the sizes of the objects are on
the order of the wavelength of light \cite{pu_intrinsic_2003,
  dixon_holographic_2011}, as is the case in colloidal suspensions.

An inverse-problem approach can yield more precise and accurate
measurements. Ovryn and Izen~\cite{ovryn_imaging_2000} and Lee
and coworkers~\cite{lee_characterizing_2007} measured positions and
optical properties of spherical colloidal particles by fitting a
Lorenz--Mie scattering solution to holograms, taking the particle
position, radius, and refractive index as fitting parameters. Fitting
holograms using scattering solutions takes advantage of known
information about the particle, such as its shape, and avoids
artifacts that lead to systematic
errors~\cite{moyses_robustness_2013}.  This inverse approach has been
used for high-precision measurements in fluid
dynamics~\cite{cheong_flow_2009} and colloid and interface
science~\cite{kaz_physical_2012}.  Recent work has extended the
fitting technique to clusters of spheres~\cite{fung_imaging_2012,
  perry_real-space_2012} and non-spherical
particles~\cite{wang_using_2014}.

However, the fitting technique is computationally intensive. For
example, Kaz, McGorty and coworkers~\cite{kaz_physical_2012} needed
about 100 CPU-hours to fit a single particle trajectory consisting of
14,000 frames. For a cluster of particles, the computational burden is
even larger: sequentially fitting 14,000 frames of a six-particle
cluster~\cite{perry_real-space_2012} would take several CPU-months.

How can we reduce the time to fit?  One strategy is to parallelize.
Although holograms can be processed more quickly using the parallel
processing inherent to a GPU~\cite{cheong_flow_2009}, it is still
difficult to process a set of such holograms in parallel.  This is
because the nonlinear optimization methods used to fit holograms
require initial guesses for the particle position, radius, and index
of refraction for each frame.  The initial guesses can be obtained
from the results of the previous frame, but only if the frames are fit
sequentially. For clusters containing many particles, automated
guesses are difficult because multiple scattering and near-field
effects make the hologram sensitive to small changes in particle
positions.

A second strategy is to reduce the size of the hologram.  Fitting a
scattering solution to a hologram is a highly over-determined problem:
typically of order 10 parameters are extracted from a hologram
containing of order 10000 pixels, each of which represents an
intensity measurement. The premise of our approach is that a far
smaller number of measurements can sufficiently constrain a fit.  This
premise was explored by Seifi and coworkers, who used staged
resampling~\cite{seifi_three-dimensional_2012} and pixel
masks~\cite{seifi_accurate_2013} to reduce the number of pixels in a
hologram and thereby reduce the number of computations required to
evaluate the scattering solution. Here we pursue an alternative
approach that requires no prior knowledge about the spatial-frequency
spectrum of a hologram.

We show that randomly sampling a subset of pixels in the hologram
reduces the fitting time by an order of magnitude or more, with little
effect on accuracy or convergence.  The random-subset approach
requires only a single sampling step, is simple to implement, works
with any scattering model, and can be combined with GPU computing or
resampling methods to reduce computing time even further.  The method
also allows us to analyze a time series of holograms in parallel.  To
do this, we first fit a small random subset of each hologram in a time
series, then use the results as initial guesses to fit the full
holograms in parallel.

\section{Background}
\label{sec:background}

We first outline the method of fitting scattering solutions to in-line
digital holograms, following Lee and
coworkers~\cite{lee_characterizing_2007}. A digital hologram
$\mathbf{H}$ captures phase information through interference of a
reference field ($E_\mathrm{ref}$) with the field scattered from an
object ($E_\mathrm{scat}$):
\begin{eqnarray}
  \label{eq:interference}
  \mathbf{H} &=& |E_\mathrm{ref} + E_\mathrm{scat}|^2\\ \nonumber
  &=& |E_\mathrm{ref}|^2 + E_\mathrm{ref}^* E_\mathrm{scat} + E_\mathrm{ref} E_\mathrm{scat}^* + |E_\mathrm{scat}|^2
\end{eqnarray}
Before fitting, we record the background intensity, which we assume to
to be equal to the reference wave intensity $|E_\mathrm{ref}|^2$. We
then calculate a normalized hologram $\mathbf{h}$:
\begin{equation}
  \mathbf{h} = \frac{\mathbf{H}}{|E_\mathrm{ref}|^2} = 1 + \frac{E_\mathrm{scat}}{E_\mathrm{ref}} +
  \frac{E_\mathrm{scat}^*}{E_\mathrm{ref}^*} + \left|\frac{E_\mathrm{scat}}{E_\mathrm{ref}}\right|^2
\end{equation}

We extract information about a physical system from $\mathbf{h}$
through an inverse-problem approach. We fit a scattering model to
$\mathbf{h}$ using a nonlinear optimization algorithm; that is, we
minimize the objective function
\begin{equation}
  \label{eq:fitting}
  g\left([x], [y], [z], [n], [r], \alpha \in \Re\right) =  \left|\mathbf{\mathbf{h_\mathrm{model}}}([x], [y], [z], [n], [r], \alpha) - \mathbf{h_\mathrm{rec}}\right|^2
\end{equation}
where the bracketed variables indicate vectors of parameters (one per
particle), $\mathbf{h_\mathrm{rec}}$ is the recorded hologram,
$\mathbf{h}_\mathrm{model}$ is the calculated hologram given parameter
  values, and $\alpha$ is an additional fitting parameter needed for
  convergence~\cite{lee_characterizing_2007}.  Depending on the shape
  and structure of the particles, other parameters may be needed.  For
  example, homogeneous nonspherical particles are characterized by
  their orientations as well as positions.  In what follows we ignore
  these additional parameters with no loss of generality.

We compute $\mathbf{h_\mathrm{model}}$ from a light scattering
solution. For a single particle we use the Lorenz--Mie
solution~\cite{lee_characterizing_2007}. For clusters of multiple
particles, we use an adaptation~\cite{fung_imaging_2012} of
Mackowski's multisphere superposition
solution~\cite{mackowski_calculation_1994}, which takes into account
all near- and far-field coupling between the scattered fields. If the
numerical minimization of the objective function converges, the result
of this fitting procedure is a measurement of the parameters $[x]$,
$[y]$, $[z]$, $[n]$, and $[r]$.

For larger clusters, we generally use fixed values for the index and
radius of the particles, as determined from manufacturer
specifications or other experiments.  The objective function $g$
therefore simplifies to:
\begin{equation}
  \label{eq:fitting_simp}
  g\left([x], [y], [z], \alpha \in \Re\right) = \left|\mathbf{h_\mathrm{model}}([x], [y], [z], \alpha) - \mathbf{h_\mathrm{rec}}\right|^2
\end{equation}

We use local optimization based on the Levenberg--Marquardt
algorithm~\cite{marquardt_algorithm_1963} to minimize $g$.  Owing to
the large parameter space and narrow minimum, local optimization
algorithms fare better than global ones. However, local approaches
require an initial guess for the set of parameters and converge only
if this guess is close to the minimum of the objective
function. Obtaining a sufficiently good guess is difficult.  For
clusters of spheres we must manually generate an initial guess for the
particle positions, using numerical reconstruction as a
guide. However, when studying a time series, we can use the results
from fitting one frame as the initial guess for the next, provided
that the capture rate is fast relative to particle motion.  This
approach minimizes manual intervention but precludes processing the
frames in parallel.

Fitting is computationally expensive because it requires repeated
evaluation of $\mathbf{h_\mathrm{model}}$, which, like
$\mathbf{h_\mathrm{rec}}$, is an array of intensity values.  The
objective function $g$ in Equation~\ref{eq:fitting_simp} is therefore
a sum over a large number of pixels:
\begin{equation}
  \label{eq:sum}
  g\left([x], [y], [z], \alpha\right) = \sum_{i, j} \left(h_\mathrm{model}(i, j; [x], [y], [z], \alpha) - h_\mathrm{rec}(i, j)\right)^2
\end{equation}
where the subscripts $i,j$ index over pixels in the recorded
hologram. To evaluate $g$, we must calculate $h_\mathrm{model}$ at
each pixel $(i,j)$ on the detector.  This corresponds to calculating
the scattered field at a set of locations in space. Calculating the
scattered field at each pixel is computationally intensive; in both
the Lorenz--Mie and multisphere superposition solutions, the field at
each point is given by a series expansion in special functions that
themselves require many operations to calculate. On a single core of a
2.67 GHz Intel Xeon X5650 CPU, it takes approximately one second to
compute a 200$\times$200 pixel hologram of a single
1-\si{\micro\meter}-diameter particle lying 15~\si{\micro\meter} from
the focal plane.  Fitting a single hologram of such a particle takes
15 to 30 seconds.  It takes much longer to fit the multisphere
superposition solution to a hologram of a cluster: about 3 seconds to
compute a hologram of a six-particle cluster of 1.3-\si{\micro\meter}
spheres, and over 500 seconds to fit a single hologram.


\section{Approach}

\begin{figure}
  \centering
  \includegraphics{figures/random-subset-fitting/holo_reprs}

  \caption{Several representations of the information in a hologram of
    six 1.3-\si{\micro\meter}-diameter polystyrene spheres in an aqueous
    suspension. All of these representations are 200$\times$200
    pixels, except (c), which is 67$\times$67 pixels.}
  \label{fig:holos}
\end{figure}

Our approach reduces the fitting time by ignoring most of the pixels
of a hologram, thereby reducing the number of evaluations required to
calculate the objective function $g$.  Consider the hologram shown in
Figure~\ref{fig:holos}a.  It contains a large number of pixels, but
only a fraction of these are needed to determine the particle position
and properties.  For example, much of the information is at low
spatial frequencies, as shown in Figure~\ref{fig:holos}b. Downsampling
techniques like those of Seifi and
coworkers~\cite{seifi_three-dimensional_2012} take advantage of the
fact that, in many holograms, the largest spatial frequencies are the
ones most affected by noise and can be discarded with little effect on
the fitted parameters. But even after downsampling
(Figure~\ref{fig:holos}c), a hologram retains symmetries and large
patches of nearly uniform value, suggesting that it can be compressed
further.

Instead of downsampling, we simply choose a random subset of the
pixels in the hologram (Figure~\ref{fig:holos}d), and fit a scattering
solution to these. This method requires little information about the
sample or the optical train, yet is surprisingly effective. The
inspiration for this technique comes from compressed sensing, where it
has been shown that random sampling is often nearly as efficient as
optimal sampling~\cite{candes_near-optimal_2006,
  donoho_compressed_2006, baraniuk_random_2009}.

\begin{figure}
  \centering
  \includegraphics{figures/random-subset-fitting/objects}
  \caption{(a, d) The holograms we use in our study (200$\times$200
    pixels, 0.12~\si{\micro\meter}/pixel). (b, e) Numerical reconstructions
    showing a plane cutting through the center of the scatterers. (c,
    f) Renderings of the structures, based on fitting a scattering
    model to the holograms. Incident light propagation is into the
    page in all cases. The inset ball and stick model is rotated to
    show structure.}
  \label{fig:objects}
\end{figure}

We examined the effectiveness of this technique on two experimentally
recorded holograms of polystyrene spheres suspended in water
(Figure~\ref{fig:objects}a,d). Both holograms we examined are
200$\times$200 pixel regions cropped from larger holograms. The first
hologram (Figure~\ref{fig:objects}a) was of a 1-\si{\micro\meter}-diameter
sphere and was originally recorded for an \emph{in situ} viscometry
experiment. We modeled this hologram using the Lorenz--Mie scattering
solution~\cite{lee_characterizing_2007}, and we fit the solution to
the data using five parameters: three spatial coordinates, the radius,
and the index of refraction. The second hologram
(Figure~\ref{fig:objects}d) is of a cluster of six 1.3-\si{\micro\meter}-diameter
particles that are nearly touching one another. It was originally
recorded for the experiments described by Perry and
coworkers~\cite{perry_real-space_2012}. We modeled this hologram using
the multisphere superposition
solution~\cite{mackowski_calculation_1994}, and we fit the solution to
the data using 18 parameters: three center coordinates for each
particle.  We fixed the index of each particle to 1.58 and the radius
to 0.645~\si{\micro\meter}. To give the fitting algorithm more freedom to
rearrange particles, we allowed the model to adopt configurations
where the particles overlapped by up to 10 percent of the radius.

We measured the convergence properties and accuracy of our method as a
function of the pixel fraction $f$ and initial guess quality.  The
accuracy we report for the fitted parameters is relative to a set of
standard values, which we calculated by fitting full frames and
averaging the fit results over different initial guesses.  The pixel
fractions varied from $f=1$, corresponding to the full hologram, to
$f=1 \times 10^{-3}$.  For the single-particle hologram, we chose 10
random subsets at each fraction, and we used an initial guess for $z$
that differed by 0.2~\si{\micro\meter} from the standard $z$, corresponding to a
displacement along the optical axis.  For the six-particle hologram,
we chose 600 random subsets at each fraction, and we used a range of
initial guesses created by adding a Gaussian random variable to each
standard center coordinate.  To vary the guess quality, we increased
the width of the Gaussian.  We generated 50 initial guesses for each
of 12 distributions whose standard deviations ranged from
$\sigma=1\times10^{-4}$~\si{\micro\meter} to $\sigma=1$~\si{\micro\meter}.  In all cases, we
rejected guesses in which particles overlapped by more than 10 percent
of their radius.

We used an implementation of the Levenberg--Marquardt algorithm from
nmpfit, a Python adaptation of LMFIT from MINPACK-1
\cite{markwardt_non-linear_2009}. We considered a fit to have
converged if the sum of squares decreased by less than $10^{-10}$
between iterations, the normalized solution vector changed by less
than $10^{-10}$ between iterations, or the cosine of the angle between
the solution and and all columns of the Jacobian was smaller than
$10^{-10}$. If none of these conditions were met within 100
Levenberg--Marquardt iterations, we halted the minimization and
considered the fit to have failed to converge.

All computations were performed on Harvard's Odyssey computing
cluster. Each fit ran on a single core of a 2.67 GHz Intel Xeon X5650
processor. We used parallel processing only to run independent trials
concurrently. Our implementation of random-subset fitting is available
in our hologram processing code HoloPy~\cite{dimiduk_holopy_2016},
available as free and open-source software.

\section{Results}

We find that random-subset fitting extracts particle parameters from
holograms more than an order of magnitude faster than fitting the full
frame.  At pixel fractions near $f=1$, the speedups we obtain increase
linearly with decreasing pixel fraction.  For both holograms, the
speedups plateau at fractions smaller than $f=0.01$
(Figure~\ref{fig:speedup}): the plateau is at 100$\times$ speedup for
the single-particle hologram and 25$\times$ for the hologram of the
six-particle cluster. The plateau occurs because the time required to
compute scattering coefficients is independent of the pixel
fraction. This overhead is larger for the multisphere superposition
solution, explaining the lower plateau for the six-particle cluster.

\begin{figure}
  \centering
  \includegraphics{figures/random-subset-fitting/speedup}
  \caption{Fitting speedup as a function of subset fraction. The
    speedup is the time required to fit the full hologram divided by
    the time required to fit a random subset, using the same initial
    guess and convergence criteria. The single-particle fits used a
    guess offset by 0.2~\si{\micro\meter} along the axis perpendicular to the
    optical axis, and the six-particle cluster fits used guesses
    offset by a Gaussian random variable with $\sigma=0.01$~\si{\micro\meter}
    added to each center coordinate of each particle. Error bars
    represent 99\% confidence intervals.}
  \label{fig:speedup}
\end{figure}

We also find that fitting such small fractions does not require better
initial guesses. For the full hologram of the six-particle cluster,
the fitting algorithm converges for initial guesses within about 100
nm (per coordinate) from the standard values, as shown in Figure
\ref{fig:convergence}. Fits converge with similar reliability for
pixel fractions as low as $f=0.01$. Reducing the fraction further
reduces the probability of convergence, but there is little reason to
use such small subsets given the diminishing returns on speedup.

\begin{figure}
  \centering
  \includegraphics{figures/random-subset-fitting/convergence}
  \caption{Fraction of six-particle cluster fits that converge, shown as a
    function of initial guess quality for subsets with $f=1$ (black),
    $f=0.1$ (red), $f=0.01$ (blue), and $f=0.001$ (green). The $f=1$
    and $f=0.1$ lines overlap.}
  \label{fig:convergence}
\end{figure}

Finally, we examine accuracy. In a simulated case with a known correct
answer (Figure \ref{fig:accuracy} dotted line) we find that we can
recover the correct particle positions to within numerical roundoff
error, even at subset fractions as low as $f=1\times 10^{-4}$ (which
correspond to only 20 pixels). Thus, in the absence of noise, the fit
accuracy is limited only by numerical precision and minimizer
cutoff. With our experimentally-measured holograms, we find a larger
deviation between the fitted parameters and the standard values at
fractions smaller than the full frame $f=1$
(Figure~\ref{fig:accuracy}). However, the errors are comparable to the
intrinsic measurement error, which is approximately $10^{-3}$~\si{\micro\meter}
for single particles~\cite{kaz_physical_2012} and $10^{-2}$~\si{\micro\meter} for
clusters~\cite{fung_measuring_2011}. Thus, we expect that fitting to
subsets as small as $f=0.1$ will not impact accuracy.

\begin{figure}
  \centering
  \includegraphics{figures/random-subset-fitting/error_vs_fraction}
  \caption{Per-particle fit error as a function of subset
    fraction. The horizontal lines show the expected intrinsic rms
    position error of fits to the single-particle hologram (dashed) and to the
    hologram of the six-particle cluster (solid). Error bars represent
    99\% confidence intervals. The rms position error in the $f=1$ fit arises
    because we fit the same hologram multiple times using different
    initial guesses. Slightly different particle configurations yield
    model holograms that match the recorded one equally well, owing to
    noise in the recorded holograms.
    \label{fig:accuracy}}
\end{figure}

\section{Discussion}

Our results show that fitting scattering models to random subsets of
holograms offers large speedups with few tradeoffs, either in
convergence or accuracy.  Two questions arise: How should it be used
as part of a workflow for analyzing holograms?  And why is random
sampling so effective?

Given how simple it is to apply random sampling, there seems little
reason not to use it when fitting holograms.  Some experimentation may
be required to find the appropriate subset fraction, which varies with
the scatterer and experimental conditions.  However, we find that
using ten percent of the pixels generally gives good results.

We have found the method to be particularly useful for analyzing time
series of holograms quickly and accurately.  We use a two stage
procedure, the first stage of which is a serial ``rough fit'' of each
hologram in the time series, using a fraction of 0.01 or 0.1.  We also
limit the minimization algorithm to two iterations. The rough fit
completes in approximately 1/200th the time required to serially fit
the full frames: we gain a factor of twenty by fitting to small
subsets, and another factor of ten by limiting the number of
iterations.  In the second stage, we use the results from the first
stage as initial guesses to fit the full holograms.  Because these
fits can be done in parallel, the speedup scales linearly with the
number of nodes.

As to the second question, on the effectiveness of the technique, we
offer only an intuitive justification here. Consider a linear
transformation, operating on a hologram, that compresses all the
information about particle parameters into a small number of basis
elements.  We expect this transformation to depend on the shape, size,
and optical properties of the scatterer.  Hence we do not---and quite
possibly cannot---know either the transformation or the optimal basis.
However, for certain classes of signals, it has been shown that random
sampling allows efficient compression without knowledge of the optimal
basis~\cite{candes_near-optimal_2006, donoho_compressed_2006,
  baraniuk_random_2009}. We cannot prove that a hologram is in this
class of signal, but the success of our method suggests that it is
close enough for practical purposes, at least for holograms of
spherical particles or clusters of spheres.  Future work might yield a
more rigorous mathematical justification of this idea. Future work
could also compare in detail to downsampling, cropping, and other
sampling other sampling methods.

\section{Conclusions}

We have shown that the random-subset technique speeds up fitting by a
factor of twenty or more, with little effect on accuracy or
convergence compared to fitting to the full hologram.  The technique
is simple to implement and works for different scattering models. It
can be used in concert with GPU or CPU-based parallel processing, with
other algorithmic improvements, or with faster implementations of
scattering models.  In our calculations we sampled from a uniform
distribution of pixel locations, but non-uniform sampling may allow
even greater speedups by weighting areas of a hologram most useful for
constraining a fit.  In future work, it will be particularly
interesting to determine if there are non-uniform distributions that
are effective for a wide variety of scatterers.

Because recovering precise 3D information from measured data is the
most time-consuming step in holographic microscopy, the method we have
demonstrated should allow faster turnaround in experiments and may
enable new experiments that were previously infeasible. The ability to
process large quantities of data more quickly will be particularly
useful for studying rare events, such as nucleation, that can only be
captured with high frame rates over long times.

\section{Acknowledgments}
This work was funded by the National Science Foundation (NSF) through
grant no.~DMR-1306410 and by the Harvard MRSEC through NSF grant
no.~DMR-0820484.  Scattering calculations and hologram fitting were
performed on the Odyssey cluster, managed by the FAS Science Division
Research Computing Group at Harvard University.

We thank Ameya Agaskar for helpful discussion on the mathematical
background of this technique.

Our nmpfit routine comes from STSCI\_PYTHON, a product of the Space
Telescope Science Institute, which is operated by AURA for NASA.
