## Abstract

Plankton interact with the environment according to their size and three-dimensional (3D) structure. To study them outdoors, these translucent specimens are imaged *in situ*. Light projects through a specimen in each image. The specimen has a random scale, drawn from the population’s size distribution and random unknown pose. The specimen appears only once before drifting away. We achieve 3D tomography using such a random ensemble to statistically estimate an average volumetric distribution of the plankton type and specimen size. To counter errors due to non-rigid deformations, we weight the data, drawing from advanced models developed for cryo-electron microscopy. The weights convey the confidence in the quality of each datum. This confidence relies on a statistical error model. We demonstrate the approach on live plankton using an underwater field microscope.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

It is important to study plankton for their role in both local and global ecosystems. These organisms are vital to the ecology of the planet [1], as they play a critical role in producing breathable oxygen (50% due to marine phytoplankton) and serve as the base of the marine food web. More recently, an increase in harmful algal blooms has been devastating to marine and freshwater supplies [2,3], while having a negative impact on coastal ecosystems.

Over decades, the basic taxonomy of plankton has been established both by morphological observations and by DNA studies that allow unique attribution. Modern oceanographers increasingly rely on *in situ* imaging systems to study populations in their natural habitat. However, almost all available inferences about plankton morphology derived from these instruments rely on two-dimensional (2D) projections [4–8] with the third dimension being conjectured. The determination of three-dimensional (3D) structure is important as it defines the *surface-area-to-volume ratio*, which is fundamental to the absorption and diffusion of nutrients and waste products. Also, flow is linked to this ratio, through the Reynolds number [9]. In addition to the 3D structure, it is important to assess the size distribution of the population (Fig. 1).

Recent advances in 3D imaging focus on lab systems [10–12]. However, many plankton are easily deformed by collection and are, thus, difficult to study in the lab. There are, therefore, a host of ocean and lake systems [13–15] to observe these organisms *in situ*. These systems look outward from their housing into the surrounding water. In the majority of systems, imaging is accomplished via the presence of organisms that simply pass through their field of view.

In order to image in relatively high resolution small plankton ($\lt\!{100}\;{\rm micrometers}$), there are four kinds of common *in situ* imaging systems that can be submerged and operated [16,17]. (i) The Imaging FlowCytobot (IFCB) [15] uses either hydrodynamic or acoustic particle focusing. This results in particles being positioned in the volume in focus, flowing through a chamber. In cases of more spherical particles, IFCB captures a single shot of a randomly sized organism [18] at a random pose. Other systems of relevance to our study are based either on (ii) holography [19–21], (iii) shadowgraph illumination [7,22], or (iv) dark-field microscopy.

Holography and shadowgraph imaging are viable modalities. However, we found several advantages to dark-field microscopy [13]: (a) Dark-field imaging is common, low cost, and robust. (b) Our experience with on-axis holography indicates that the point spread function is range-dependent (a quadratic phase effect) that generally yields lower contrast than dark-field. We also found, in the case of holography, difficulty in yielding viable depth information at the thicknesses of the organisms that we study. A dark-field image is a simple projection rather than a phase measurement. Projected images can provide information that can yield 3D structure using tomography, as we describe in this paper. (c) Small scattering particles in turbid waters, as are commonly found in coastal regions, yield random localized specks in dark-field images. This randomness can be filtered out by tomography (which is based on multiple images). On-axis holography (the practical basis for digital holography) responds to scatterers by non-local speckled patterns, which reduce overall contrast. (d) Dark-field illumination enhances edges of objects whose refractive index is close to that of the surrounding medium. This leads to better resolved appendages and internal structures of plankton. We, thus, chose to work with dark-field optical microscopes to maximize taxonomic resolution, enhance edge contrast, and ensure high performance in scattering-rich environments.

Underwater computer vision has benefited from a dramatic increase in computing capabilities, addressing image quality [23–27] and underwater robotics [28–33]. It is fitting that the benefits of computational imaging should be applied to *in situ* study of plankton. Tomography is a computational method that we use. Usually, tomography is performed around a single, temporally fixed object [34–38]. Here we are interested in drawing structure and statistics about a *population*. The data are an ensemble of random views of random objects within a population, each having a different random and unknown scale, orientation, and location (see Fig. 1). A basis to our solution is *single particle reconstruction* (SPR) [39], which is used in cryo-electron microscopy (Cryo-EM) [40,41]. SPR handles random rotations but not scale variation in a population. An attempt to generalize SPR to handle variable scales was done in Ref. [42], where the approach termed 3D-POP is defined. It is based on pairwise constraints on rotations and scales. However, the work in Ref. [42] is prone to significant estimation errors. Hence, our paper employs a more robust SPR approach [43] and generalizes it to analyze data subject to scale variations. The proposed method uses a set of *triplets* of images to assess how indicative are constraints derived from image pairs. Then, this assessment is used to weight each constraint. The constraints are then efficiently and optimally fused using linear algebra.

We show the effectiveness of the method on simulated data and publically available bright-field images from the IFCB. We also show it using dark-field images of plankton acquired by our *in situ* system: the Scripps Plankton Camera (SPC). SPC is an *in situ* *dark-field* microscope developed to continuously monitor plankton populations [44] (Fig. 2). It is a free-space imaging system, with the light and camera focused on a plane in an open volume of fluid. The SPC uses a white LED illumination source with the central rays blocked by an opaque disk. Objects in the focal plane, thus, scatter light toward the image plane while the direct light continues past the lens, thus creating the dark-field effect. The paper uses images from two test SPC datasets. Both were taken by a version of the SPC outfitted with a ${0.5} \times$ microscope objective, yielding an object sampling interval of 7.4 µm/pixel with adequate contrast at 30 µm. All presented SPC plankton samples are about 0.5–2.5 mm long. The object sampling interval of the IFCB is 0.29 µm/pixel, enabling imaging of objects in an order of magnitude smaller than the SPC. All plankton samples were manually classified by an expert taxonomist.

## 2. BACKGROUND ON SINGLE PARTICLE RECONSTRUCTION

#### A. Tomography Using Random Projections

We denote 3D spatial coordinates by ${\textbf x} = ({x_1},{x_2},{x_3}{)^{T}}$, where ${(\cdot)^{T}}$ denotes transposition. A specimen has a volumetric optical extinction distribution $\beta ({\textbf x})$. In dark-field imaging, $\beta ({\textbf x})$ expresses the density of scattering particles in the object. In bright-field imaging, $\beta ({\textbf x})$ represents the optical density of attenuating media.

Computational tomography (CT) estimates $\beta ({\textbf x})$ using a set of $N$ multi-view orthographic projections. Relative to the camera, the sample orientation in the $n$th image ${{\textbf y}_n}$ is expressed by a 3D rotation matrix ${{\textbf R}_n}$. Orthographic projection through the translucent sample is expressed by linear integration along the ${x_3}$ axis,

There are strong and widely used methods for estimating a volumetric distribution $\beta ({\textbf x})$ using tomography based on such a set of projections $\{{{\textbf y}_n}\} _{n = 1}^N$. This paper, thus, does not focus on estimating $\beta ({\textbf x})$. We note, however, that common tomographic methods [47–54] assume that the set of $N$ orientation matrices $\{{{{\textbf R}_1}, \ldots ,{{\textbf R}_N}} \}$ is known. Moreover, common tomography assumes that the specimen itself does not change between projections.

A random biological microscopic specimen collected *in situ* does not comply with these assumptions. Its orientation is random and unknown. Furthermore, each projection image corresponds to a *different specimen*, hence having random, unknown scale and lateral displacement. This paper focuses on recovering the unknown scale, displacement, and rotation of each specimen, in order to facilitate tomographic recovery of the 3D structure.

#### B. 3D-POP

There is a prior method for estimating the unknown scale, displacement, and rotation of each specimen: 3D-POP [42]. Our paper presents methods that yield better efficiency, accuracy, and robustness. We briefly describe the main elements of 3D-POP, drawing significantly on its notations.

Analysis is done in the 3D spatial frequency domain ${\textbf k} = ({k_1},{k_2},{k_3})$. The 3D Fourier transform of the volumetric object is ${{\cal F}_{\textbf k}}\{\beta \}$. A *central slice* is a frequency-plane in the 3D Fourier space that passes through the origin, ${\textbf k} = 0$. Any central slice is a 2D domain, which can be written in terms of a *rotated* horizontal slice,

Now let us consider the 2D Fourier transform of projection image ${{\textbf y}_n}$. Denote the resulting transformed representation by ${{{\tilde {\textbf y}}}_n}(\rho ,\psi)$, where $\rho$, $\psi$ are polar coordinates in the 2D spatial frequency domain. According the *Central Slice Theorem*, the 2D function ${{{\tilde {\textbf y}}}_n}$ is equivalent to the central slice ${\textbf R}_n^{- 1}{[{k_1}, {k_2}, 0]^{T}}$ in 3D $ k $-space Eq. (2).

Any two central slices $ n $, $ m $ intersect at a *common line* in 3D k-space (Fig. 3). Correspondingly, this line in 3D $ k $-space is equivalent to a radial line in the 2D function ${{{\tilde {\textbf y}}}_n}$ at an angle denoted by $\psi = {\psi _{n \to m}}$. It is also equivalent to a radial line in the 2D function ${{{\tilde {\textbf y}}}_m}$, at an angle denoted by $\psi ^\prime = {\psi _{m \to n}}$.

Within these 2D functions, finding the angle pair ${\psi _{n \to m}}$, ${\psi _{m \to n}}$, which corresponds to a common line, yields constraints [41] about the 3D rotations ${{\textbf R}_n}$, ${{\textbf R}_m}$. The constraints are

- (a) Estimating all projection rotations $\{{{\textbf R}_n}\} _{n = 1}^N$ using constraints as in Eq. (3) is a key for subsequent tomographic recovery of $\beta ({\textbf x})$.
- (b) To yield constraint as Eq. (3), the angle pair ${\psi _{n \to m}}$, ${\psi _{m \to n}}$ must be estimated reliably and efficiently for many real-data image-pairs. This challenge is exacerbated by the fact that $\beta ({\textbf x})$ is randomly shifted and scaled between images of different organisms collected
*in situ*.

Let us start with challenge (b). In the 2D Fourier domain, the one-dimensional (1D) function on a radial line ${{{\tilde {\textbf y}}}_n}(\rho ,{\psi _{n \to m}})$ is expected to highly correlate with the radial function ${{{\tilde {\textbf y}}}_m}(\rho ,{\psi _{m \to n}})$ since both correspond to the same common line. The correlation should generalize to account for a relative scale ${M_{n \to m}}$ between the specimens viewed in images $n$ and $m$. Relative spatial scale is expressed as an inverse-scale in the Fourier domain. The correlation should also generalize to account for a lateral shift ${\tau _{n \to m}}$ of the specimen along the common line direction. This shift yields a linear phase $2\pi \rho {\tau _{n \to m}}$ in the Fourier domain. Consequently, in Ref. [42] the Fourier domain radial correlation is defined by an inner product of scale and phase-shifted 1D functions sampled along two directions, $\psi ,\psi ^\prime $, using variable scale $M$ and lateral shift $\tau$:

Consequently, a data term that penalizes scale inconsistencies is

### 1. Problems That Persist

There are several problems with Eqs. (4) and (5) and their implementation. In Ref. [42], Eq. (5) was done by exhaustive search over a four-dimensional (4D) domain of variables. To make the search space manageable, Ref. [42] used a coarse-to-fine approach starting from a sparse sampling of the domain of variables. Equation (5) is performed on this small sparse set. This results in a narrower 4D domain to be searched in finer resolution. However, noise and other factors may lead to errors in Eq. (5). An error occurring in a coarse-sampling stage can lead to intolerable errors that cannot be solved in later, finer-scale sampling. In effect, robustness and accuracy are sacrificed for numerical tractability.

The model in Eqs. (4) and (5) assumes a similarity relation between specimens. It does not account for non-rigid deformations, which are common in natural non-rigid objects in a population. A non-rigid deformation is a perturbation that may lead to deviation from the estimated pose variables. It may also lead to an inaccurate solution of Eq. (5): by chance, correlation of an arbitrary orientation scale and translation can yield better correlation [Eq. (4)] than the correct pose due to the deformation. Hence, a major error in pose propagates to major errors in tomographic recovery. In this paper, we use a more robust approach.

The biggest challenge is that of recovering the orientation. Based on the pairwise relations as in Eqs. (3) and (5), all rotations were estimated in Ref. [42] using a computationally complex, non-convex, constrained optimization as follows:

## 3. STATISTICAL MODEL FOR ROBUST PLANKTON RECONSTRUCTION

#### A. Estimating Rotations

To mitigate this complexity, we describe a numerically efficient method that solves $\{{{{\hat {\textbf R}}}_n}\} _{n = 1}^N$ given pairwise relations. Contrary to Eq. (8), the method here is not formulated as a non-convex problem or using approximate relaxation. Instead, it uses a linear algebra approach [43] to achieve accurate results quickly. Define

Use the angle pair ${\hat \psi _{n \to m}}$, ${\hat \psi _{m \to n}}$ and a matrix ${{\textbf R}_{m \to n}}$, which satisfies Eqs. (3) and (9) as found as in Ref. [55]. Then, define matricesUsing this formulation, we propose an algorithm for empirically estimating the rotations matrices $\{{{{\hat {\textbf R}}}_n}\} _{n = 1}^N$.

- 2. Perform singular value decomposition (SVD) on ${\textbf S}$. Select the three largest singular values and their corresponding eigenvectors. Denote these column vectors as ${{\textbf h}_1}$, ${{\textbf h}_2}$, ${{\textbf h}_3}$.
- 3. Construct the estimated rotations by concatenating these column vectors into a matrix:

We describe in Sections 3.B and 3.C an approach to yield a *confidence weight* ${w_{n,m}}$ to each ${\hat \psi _{n \to m}}$, ${\hat \psi _{m \to n}}$ pair. A low value of ${w_{n,m}}$ indicates a low confidence in the correctness of this pair of angles. The weights are normalized so as to satisfy

#### B. Statistical Model of Common Lines

To derive weights ${w_{n,m}}$, we need to have a measure of confidence or uncertainty for each common line. The uncertainty depends on a probabilistic model of common lines, the noise sources that affect this model and finite sampling effects. From our experience, errors in ${\hat \psi _{n \to m}}$, ${\hat \psi _{m \to n}}$ cause problems that are more severe than problems stemming from errors in scale or shift. Hence, we focus our attention on angular errors.

Consider the 2D Fourier representation ${{{\tilde {\textbf y}}}_n}$ (Fig. 4). The representation of the true common line with ${{{\tilde {\textbf y}}}_m}$ is a line at true angle ${\psi _{n \to m}}$. Its empirical estimation is ${\hat \psi _{n \to m}}$. The corresponding angular error is

The direction of any common line in each 2D Fourier domain is estimated only modulo 180°. The reason is that the correlation [Eq. (4)] between any two such Fourier domain lines is invariant to reflection around the zero-frequency. Thus, $\Delta {\psi _{nm}} \in [- {90^ \circ}{,90^ \circ})$.In Cryo-EM literature [43], a parametric model had been suggested to express the distribution of $\Delta {\psi _{nm}}$. The model dictates the following:

- • There is probability $ P $ that a pair ${\hat \psi _{n \to m}}$, ${\hat \psi _{m \to n}}$ is
*indicative*: the angles are random though roughly correct. The probability is assumed to be Gaussian, in the form $\Delta {\psi _{nm}}$, $\Delta {\psi _{mn}} \sim {\cal N}(0,{\sigma ^2})$. Here $\sigma$ is the standard deviation of angular errors in indicative common lines. - • There is probability $1 - {P}$ that a pair ${\hat \psi _{n \to m}}$, ${\hat \psi _{m \to n}}$ is
*arbitrary*: each of these angles is random, distributed uniformly around a circle.

Therefore, the probability density function (PDF) of $\Delta {\psi _{nm}}$ is modeled by

### 1. Effect of Radiometric Noise

We simulated a synthetic volume and projected it using $N = 100$ random orientations. The projections included photon (Poisson) noise, typical to optical cameras. The best quality had photon noise that corresponded to a full-well depth of 10,000 photoelectrons, i.e., ${\rm SNR} = {100}$ near camera saturation. The lower quality corresponded to a full-well depth of 100 photoelectrons, i.e., ${\rm SNR} = {10}$ near camera saturation. Equations (4), (5), and (17) were then run, after which the model of Eq. (18) was fitted to the histogram of errors. The fit was good yet *insensitive to the radiometric SNR*, yielding ${P} = 0.897 \pm 0.007$ and $\sigma = {1.25^\circ} \pm {0.01^ \circ}$ across the range. Similar results were obtained when introducing random lateral displacements and scales to the imaged samples. For this reason, we do not see radiometric noise as a major source of uncertainty in common line estimation, in our type of imaging.

### 2. Specimen Deformation as Noise

The objects we consider always have variance in sub-structures, i.e., non-rigid deformations. Such random morphological noise disrupts proper detection and matching of common line pairs. To induce such noise in a simulation, we followed a recipe by Ref. [42]. Let the longest dimension of a volumetric distribution $\beta ({\textbf x})$ be bounded by $G$ voxels. The domain size of $\beta ({\textbf x})$ is $G \times G \times G$ voxels. This distribution is deformed to $\beta ({\textbf x} + {{\textbf v}^{(n)}}[{\textbf x}])$. Here ${{\textbf v}^{(n)}} = [{v_1},{v_2},{v_3}]$ represents a random vector (deformation) field of the specimen corresponding to image $n$. There is a 3D vector of deformation at every point in the object.

The deformation direction ${{{\hat{\textbf v}}}^{(n)}} = [\hat v_1^{(n)},\hat v_2^{(n)},\hat v_3^{(n)}]$ is randomly sampled on the unit sphere. The deformation magnitude varies spatially as a sinusoid having a random amplitude ${V^{(n)}} \sim {\cal U}(0,{V_{\max}}G)$. Each spatial component $c \in \{1,2,3\}$ has random phase $\phi _c^{(n)} \sim {\cal U}(- \pi ,\pi)$. The distortion of component $c$ is

This deformation field model represents a basic Fourier component of distortion. We set $G = 100$. For each value of ${V_{\max}}$, we computed $N = 100$ random projections and deformations. Equations (4), (5), and (17) were run, and then the model of Eq. (18) was fitted to the histogram of angular errors (Fig. 5). Across a range of ${V_{{\max}}}$ values, the fit is good. The values of $ P $ and $\sigma$ are sensitive to the magnitude of distortion amplitude ${V_{{\max}}}$ (Fig. 6).#### C. Common Line Triplets

To establish the weights ${w_{n,m}}$, we need a statistical model. Section 3.B presented a generic model parameterized by $ P $ and $\sigma$. When using real-data, we do not know the values of these statistical parameters. Thus, we derive these parameters based on the actual real-data, which includes many (noisy) estimates of angles ${\{{\hat \psi _{n \to m}}\} _{n,m}}$.

We denote the $ n $, $ m $ pair of images as $\overline {n,m}$. For any $\overline {n,m}$, we should assess the likelihood that the estimated angle pair is indicative or arbitrary. This can be done by relying on *common line triplets*. Denote a different triplet of projection images (indexed $n,m,l$) of the same object as $\overline {n,m,l}$. It corresponds to true rotations ${{\textbf R}_n} ,{{\textbf R}_m} ,{{\textbf R}_l}$, and consequently to true three matrices ${{\textbf R}_{n \to m}}$, ${{\textbf R}_{m \to l}}$, ${{\textbf R}_{l \to n}}$, following Eq. (9). To be correct rotation matrices, they must satisfy

### 1. Triplet Score and Statistics

Reference [43] proposed a *score* to quantify how closely Eq. (20) is approximated using estimated matrices,

Because the empirical values of ${\hat \psi _{n \to m}},{\hat \psi _{m \to l}},{\hat \psi _{l \to n}}$ are random, then ${s_{n,m,l}}$ is random. If all three angles ${\hat \psi _{n \to m}},{\hat \psi _{m \to l}},{\hat \psi _{l \to n}}$ are indicative, then $\overline {n,m,l}$ is said to be an *indicative triplet*, and ${s_{n,m,l}}$ is likely to be high. However, if any of these angles are arbitrary, then $\overline {n,m,l}$ is said to be an *arbitrary triplet*, and ${s_{n,m,l}}$ is likely to be low.

These likelihoods are expressed by PDFs [43]. For *arbitrary triplets*, the PDF of $s$ is denoted by $f_{{\rm triplet}}^{{\rm arbitrary}}(s)$, which can be pre-calculated using a simulation. This simulation does not rely on any image data. Rather, it tests relations between random rotation matrices:

- (a) Draw a set of a large number of random uniformly distributed common line angles ${\{\hat \psi _{n \to m}^{{\rm simul}}\} _{m,n}}$.
- (b) Calculate $s \in [0,1]$ for triplets drawn from this set using Eq. (21).
- (c) Calculate a normalized histogram of $s$ and fit a simple continuous curve [43] to it, in the domain $s \in [0,1]$.

This yields $f_{{\rm triplet}}^{{\rm arbitrary}}(s)$, and it is independent of the particular mechanism underlying any arbitrary orientation errors. The score is maximal for $s = 0$ and decays to zero at $s = 1$ (Fig. 7).

For *indicative triplets*, the PDF of $s$ is analogously denoted by $f_{{\rm triplet}}^{{\rm indicative}}(s|\sigma)$. Here $\sigma$ is the standard deviation of angular errors in indicative common lines, as used in Section 3.B.

The PDF $f_{{\rm triplet}}^{{\rm indicative}}(s|\sigma)$ is likewise pre-calculated in simulations, irrespective of the underlying mechanism causing orientation errors:

- (a) Draw a set of a large number of arbitrary orientations ${\{{{\textbf R}_n}\} _n}$, randomly uniformly distributed.
- (b) Calculate the set of corresponding correct common line angles ${\{\psi _{n \to m}^{{\rm simul}}\} _{m,n}}$.
- (c) Perturb these angles by random Gaussian errors having standard deviation $\sigma$ [as in Eq. (18)] for ${P} = 1$. This yields a set of indicative ${\{\hat \psi _{n \to m}^{{\rm simul}}\} _{m,n}}$.
- (d) Calculate $s \in [0,1]$ for triplets drawn from this set, using Eq. (21).
- (e) Calculate a normalized histogram of $s$ and fit a simple continuous curve [43] to it, in the domain $s \in [0,1]$. This fitted curve is parameterized by $\sigma$.

For the special case $\sigma = 6$, $f_{{\rm triplet}}^{{\rm indicative}}(s|\sigma)$ is indeed maximal close to $s = 1$ and decays to zero at $s = 0$ (Fig. 7).

### 2. Estimating the Statistical Model Parameters

To estimate $P$ and $\sigma$, we use the pre-calculated PDFs $f_{{\rm triplet}}^{{\rm arbitrary}}(s)$, $f_{{\rm triplet}}^{{\rm indicative}}(s|\sigma)$ together with the set ${\{{\hat \psi _{n \to m}}\} _{m,n}}$ empirically derived from the real-data. Consequently the weights ${w_{n,m}}$ are derived as in Section 3.C.3. Recall that for each $ m $, $ n $, the probability that $\{{\hat \psi _{n \to m}}\}$ is indicative is $ P $. Assume that the estimation of each common line is approximately independent of the estimation of the other lines. Hence, the probability that a triplet is indicative is $\approx {{P}^3}$, and the probability that a triplet is arbitrary is approximately $1 - {{P}^3}$.

From the data-based estimation of ${\{{\hat \psi _{n \to m}}\} _{m,n}}$, the matrices ${\{{{{\hat {\textbf R}}}_{n \to m}}\} _{m,n}}$ are derived. Then Eq. (21) yields a set of scores $\{{s_{n,m,l}}\} _{n \ne m \ne l = 1}^N$. The values in the set of scores yield a normalized histogram ${\cal H}(s)$. According to the model,

### 3. Confidence Weights

To derive the weights ${w_{n,m}}$, we can now rely on the estimated ${{\hat P}},\hat \sigma$ and the set of empirical scores $\{{s_{n,m,l}}\} _{n \ne m \ne l = 1}^N$. The confidence in the angles of $\overline {n,m}$, ${\hat \psi _{n \to m}},{\hat \psi _{m \to n}}$ is assessed by the set of triplet scores associated with $\overline {n,m}$:

The likelihoods that this set is drawn from indicative or arbitrary common lines are $\ell _{n,m}^{{\rm indicative}}(\hat \sigma)$ and $\ell _{n,m}^{{\rm arbitrary}}$, respectively. Details of estimates of $\ell _{n,m}^{{\rm indicative}}(\hat \sigma)$ and $\ell _{n,m}^{{\rm arbitrary}}$ are given in Appendix A. Given $\ell _{n,m}^{{\rm indicative}}(\hat \sigma)$ and $\ell _{n,m}^{{\rm arbitrary}}$, the likelihood of the scores ${{\cal S}_{n,m}}$ is ${{\hat P}}\ell _{n,m}^{{\rm indicative}}(\hat \sigma) + (1 - {{\hat P}})\ell _{n,m}^{{\rm arbitrary}}$. The probability density of $ n $, $ m $ being indicative, given the set ${{\cal S}_{n,m}}$, is, thus [43],From this result, a higher weight is given to a common line pair that yields a higher likelihood of indicative triplet scores. This means that the particular common line being weighted has high support across the entire set of common lines. A common line pair, which has little support from the other common lines for being indicative, has a low probability for being indicative and thus has a lower weight. The weight [Eq. (25)] is used for constructing the matrix $\widetilde {\textbf S}$ (Section 3.A). In turn, this enhances the estimation of the rotation matrices $\{{{{\hat {\textbf R}}}_n}\} _{n = 1}^N$.

#### D. Constraint for Scale Estimation

Recall from Eq. (5) that pairwise rotation angles, scales, and translations are searched in a 4D domain. This is a computationally expensive process that is prone to errors in estimated pairwise rotation angles, which can in turn induce errors in scale. Our procedure (Sections 3.A–3.C) reduces these problems, yielding a robust estimation of the rotations, despite uncertainties and significant outliers. The estimated rotations $\{{{{\hat {\textbf R}}}_n}\} _{n = 1}^N$ consequently *constrain the pairwise angles*. The constrained angles then lead to a simpler search, whose results proved to be more accurate. We follow these steps:

- (a) Run Eq. (5) for initial values of ${\{{\hat \psi _{n \to m}}\} _{m,n}}$, as well as initial estimates of pairwise scales and shifts. The search is executed on a coarse 4D grid for the sake of speed. The results can have many outliers with large uncertainties. Discard the initial estimates of pairwise scales and shifts.
- (c) Use $\{{{{\hat {\textbf R}}}_n}\} _{n = 1}^N$ to derive the common line between any pair of central slices $\overline {n,m}$. This is done in closed-form and yields a new set of common line pairs, denoted by ${\{{\tilde \psi _{n \to m}}\} _{m,n}}$. The new set of common lines ${\{{\tilde \psi _{n \to m}}\} _{m,n}}$ is, thus, a result of robust estimation, which accounts for rotations of all projections in the image set. This is contrary to the set ${\{{\hat \psi _{n \to m}}\} _{m,n}}$, which was crudely sampled based on pairs of images.
- (d) Constrain Eq. (5) to run on the common line pair ${\tilde \psi _{n \to m}},{\tilde \psi _{m \to n}}$:

The resulting pairwise scales are then used in Eqs. (6) and (7), to yield all scales, globally. Similarly, pairwise translations are fed into a linear optimization process. This yields all translations in a global frame, akin to [42]. Once rotations, scales, and translations are estimated, the pairwise search is repeated. In the new iteration, the search domain is narrowed and sampled more finely. This step is integrated into a broader coarse-to-fine approach that follows [42].

Eventually, an estimate of all orientation, scales, and shifts is obtained. Consequently, tomographic 3D recovery of the objects is performed, yielding $\hat \beta ({\textbf x})$. We used an algebraic reconstruction technique [56,57] for this purpose.

## 4. RESULTS

#### A. Quantitative Criteria

We performed simulations as well as experiments using real-data. For quantitative criteria, we follow Ref. [42] by noting that all rotations can only be estimated up to an unknown global rotation ${\textbf O}$. Similarly, all scales are defined up to a global unknown ${M_{{\rm rel}}}$. Hence [42], when ground-truth is available during simulations, error measures for rotations and scales are, respectively,

*in situ*plankton images have no ground-truth, quality is assessed via cross-validation. Rotations, scales, and translations are estimated using the entire ensemble. Then a single image $\tilde n$ is excluded from the $N$-sample ensemble. The 3D estimation $\hat \beta$, thus, uses $N - 1$ images. Then, $\hat \beta$ is re-projected using ${{{\hat {\textbf R}}}_{\tilde n}},{\hat M_{\tilde n}},{\hat \tau _{\tilde n}}$, yielding a synthetic image ${{{\hat {\textbf y}}}_{\tilde n}}$. It is compared to the true image ${{\textbf y}_{\tilde n}}$ using relative root mean square error. This is repeated $\forall \tilde n \in [1,N]$. There are ${N_{{\rm pixel}}}$ pixels. The re-projection error is

#### B. Simulation Based on Empirical Data

The simulation is based on a real object whose 3D volumetric distribution is empirically known [59]. In our simulations, this object underwent random deformations as described in Section 3.B.2, with ${V_{{\max}}} = 6\%$. It also underwent random scaling distributed as $\mathop {\log}\nolimits_{e} ({\rm scale}) \sim {\cal U}({- 0.7,0.7})$. This results in a scale factor of $\times 4$ between the lengths of the largest and smallest specimen. The object was projected to $N = 100$ random directions, each yielding an image of size $126 \times 126$ pixels. The projection then underwent random lateral shifts in the image domain, distributed as ${\rm shift} \sim {\cal U}({- 10,10}) [{\rm pixels}]$ in both image axes. The prototype object and sample projections are shown in Fig. 8.

The initial search domain in Eq. (5) used sampled angles with one-degree resolution. The domain had 10 relative scale samples uniformly spaced in the logarithm domain $[- \log (2) \log (2)]$ and 10 samples of relative offsets in a $[- 15 15]\; {\rm pixel}$ domain. As in Section 3.D, the domains become narrow through an iterative coarse-to-fine search. The initial pairwise matches were processed using the approach of this paper and compared to the result of Ref. [42]. Results are shown in Fig. 9 and Table 1.

#### C. Microscopy Experiments

### 1. Dark-Field Data

The first dataset contains images of *Daphnia* sp., a species of freshwater zooplankton (Fig. 10). Outliers are first pruned out [42]: images whose estimated common lines pairs have low correlation values (under a manually chosen threshold) are considered outliers. For this fixed threshold, pruning resulted in a cluster of $N = 147$ samples for [42]. Correlations on the constrained lines described in Section 3.D yielded $N = 145$ inlier samples for our approach. The initial search domain in Eq. (5) is as in the simulations. Comparing the resulting re-projections on sample frames from both approaches indicates that our approach is able to resolve finer details (Fig. 10). Quantitatively, our approach yields ${\epsilon_{{\rm re} {\text -} {\rm proj}}} = 1.38 \pm 0.002$. In comparison, the method of Ref. [42] yields ${\epsilon_{{\rm re} {\text -} {\rm proj}}} = 3.065 \pm 0.004$.

The second dataset consists of 1000 images of the cosmopolitan zooplankton *Oithona* sp. (Fig. 11). To keep runtime short in this large database, we used 50 equally spaced angular samples in the search domain of Eq. (5). We used the other domain sample parameters as those in the simulations. Pruning yielded $N = 921$ inlier samples for our approach and $N = 877$ for the method of Ref. [42]. Figure 11 compares re-projections while Fig. 12 compares 3D recoveries.

The proposed method yields ${\epsilon_{{\rm re} {\text -} {\rm proj}}} = 4.174 \pm 8 \cdot {10^{- 4}}$. In comparison, the method of Ref. [42] yields ${\epsilon_{{\rm re} {\text -} {\rm proj}}} = 6.559 \pm 0.002$.

Recall that our method yields an estimation of the size distribution. Figures 13(a) and 13(b) show the estimated size histograms of *Oithona* sp. and *Daphnia* sp. and the corresponding PDFs fitted using Eq. (31). The estimated size histograms and fitted PDF by Ref. [42] are very similar to ours.

### 2. Bright-Field Data

For direct comparison to the results of Ref. [42], we used the publicly available plankton image database provided by Ref. [60], specifically focusing on *Pyramimonas longicauda*. We note that this algal species is an order of magnitude smaller than the zooplankton considered in the other datasets. Pruning of poorly correlated data yielded a set of $N = 136$ images for our analysis and $N = 132$ for running [42]. Figure 14 visualizes the respective 3D recoveries.

Figure 15 shows several re-projections and compares them to the real images. Both figures indicate that the proposed approach yields more complete results, of higher contrast.

Quantitatively, our approach yielded ${\epsilon_{{\rm re} {\text -} {\rm proj}}} = 3.446 \pm 0.004$. In comparison, the method of Ref. [42] yielded ${\epsilon_{{\rm re} {\text -} {\rm proj}}} = 3.942 \pm 0.005$.

## 5. DISCUSSION

Using our method, massive 2D image databases collected globally by ocean and lake *in situ* observation systems can lead to 3D structure and population statistics estimates. Tackling this problem at scale necessitates more efficient processing. For example, using 3D volumetric models of plankton, learning-based methods can be trained to infer species and distributions, based on image projections taken at random poses.

Our experiments indicate that the method is robust to differences in organism size: the SPC images are of plankton in the 0.5–2.0 mm range, while those captured by the IFCB are tens of microns across. The method seems to be robust to the body plan: *Daphnia* and *Oithona* are both zooplankton while *Pyramimonas* is a photosynthetic algae. We believe the approach will yield reasonable reconstructions for a variety of plankton, though domain experts should be consulted when approaching new taxa. Gelatinous organisms will be particularly challenging since much of their bodies are transparent in seawater.

The approach presented here is not limited to plankton. We believe it can be used in population analysis of a wide variety of small, translucent suspended objects. Beyond other organism classes, the method might also be applied to non-biological targets. These may include submerged suspended sediments, airborne aerosols and hydrosols, ice crystals in clouds, and possibly small meteoroids.

## APPENDIX A

We now describe how estimates of $\ell _{n,m}^{{\rm indicative}}(\hat \sigma)$ and $\ell _{n,m}^{{\rm arbitrary}}$ are obtained, for use in Eq. (24).

- • Assume for a moment that the $ n $, $ m $ common line is arbitrary. Then, by definition, all elements in ${{\cal S}_{n,m}}$ are scores of arbitrary triplets. Then, the likelihood of this set of scores being drawn from arbitrary common lines is$$\ell _{n,m}^{{\rm arbitrary}} \equiv {\ell ^{{\rm arbitrary}}}({{\cal S}_{n,m}}) = \prod\limits_{l \notin \{n,m\}} f_{{\rm triplet}}^{{\rm arbitrary}}({s_{n,m,l}}).$$Here, it is assumed that ${s_{n,m,l}},{s_{n,m,l^\prime}}$ are approximately statistically independent $\forall l \ne l^\prime $,
*given*that the common line of $\overline {n,m}$ is arbitrary. - • Assume for a moment that the common line of $\overline {n,m}$ is indicative. Then, $\overline {n,m,l}$ is indicative only if both the common lines of $\overline {n,l}$ and $\overline {m,l}$ are indicative as well. For such an event, the likelihood of the score ${s_{n,m,l}}$ is, thus, ${{{\hat P}}^2}f_{{\rm triplet}}^{{\rm indicative}}({s_{n,m,l}}|\hat \sigma)$.

The complementing case is that either or both common lines of $\overline {n,l}$ and $\overline {m,l}$ are arbitrary. Then, the triplet $\overline {n,m,l}$ is arbitrary, despite the common line of $\overline {n,m}$ being indicative. In this case, the likelihood of the score ${s_{n,m,l}}$ is, thus, $(1 - {{{\hat P}}^2})f_{{\rm triplet}}^{{\rm arbitrary}}({s_{n,m,l}})$. Consequently, in total, if the common line of $\overline {n,m}$ is indicative, the likelihood of the score ${s_{n,m,l}}$ is ${{{\hat P}}^2}f_{{\rm triplet}}^{{\rm indicative}}({s_{n,m,l}}|\hat \sigma) + (1 - {{{\hat P}}^2})f_{{\rm triplet}}^{{\rm arbitrary}}({s_{n,m,l}})$. Therefore, compounding all scores in ${{\cal S}_{n,m}}$, the likelihood of this set being drawn while the common line of $\overline {n,m}$ is indicative is

## Funding

European Research Council (810370).

## Acknowledgment

We thank Aviad Levis and Tali Treibitz for their advice, and Johanan Erez, Ina Talmon, and Daniel Yagodin for technical support. Yoav Schechner is the Mark and Diane Seiden Chair in Science at the Technion. He is a Landau Fellow—supported by the Taub Foundation. His work was conducted in the Ollendorff Minerva Center. Minvera is funded through the BMBF. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 810370-ERC-CloudC.

## Disclosures

The authors declare no conflicts of interest.

## Data Availability

Data underlying the results presented in this paper are available in Refs. [4,60].

## REFERENCES

**1. **G. C. Hays, A. J. Richardson, and C. Robinson, “Climate change and marine plankton,” Trends Ecol. Evol. **20**, 337–344 (2005). [CrossRef]

**2. **J. Ryan, R. Kudela, J. Birch, M. Blum, H. Bowers, F. Chavez, G. Doucette, K. Hayashi, R. Marin, C. Mikulski, and J. Pennington, “Causality of an extreme harmful algal bloom in Monterey Bay, California, during the 2014–2016 northeast Pacific warm anomaly,” Geophys. Res. Lett. **44**, 5571–5579 (2017). [CrossRef]

**3. **C. Anderson, E. Berdalet, R. Kudela, C. Cusack, J. Silke, E. O’Rourke, D. Dugan, M. McCammon, J. Newton, S. Moore, and K. Paige, “Scaling up from regional case studies to a global harmful algal bloom observing system,” Front. Mar. Sci. **6**, 250 (2019). [CrossRef]

**4. **“Scripps plankton camera website,” http://spc.ucsd.edu.

**5. **G. Gorsky, M. Picheral, and L. Stemmann, “Use of the underwater video profiler for the study of aggregate dynamics in the North Mediterranean,” Estuarine Coastal Shelf Sci. **50**, 121–128 (2000). [CrossRef]

**6. **C. S. Davis, S. M. Gallager, M. Marra, and W. K. Stewart, “Rapid visualization of plankton abundance and taxonomic composition using the video plankton recorder,” Deep Sea Res. II **43**, 1947–1970 (1996). [CrossRef]

**7. **M. D. Ohman, R. E. Davis, J. T. Sherman, K. R. Grindley, B. M. Whitmore, C. F. Nickels, and J. S. Ellen, “Zooglider: an autonomous vehicle for optical and acoustic sensing of zooplankton,” Limnol. Oceanogr. Methods **17**, 69–86 (2019). [CrossRef]

**8. **R. K. Cowen and C. M. Guigand, “In situ ichthyoplankton imaging system (*I* SIIS): system design and preliminary results,” Limnol. Oceanogr. Methods **6**, 126–132 (2008). [CrossRef]

**9. **S. Vogel, *Life in Moving Fluids: The Physical Biology of Flow-Revised and Expanded Second Edition* (Princeton University, 2020).

**10. **S. Lee, K. Kim, A. Mubarok, A. Panduwirawan, K. Lee, S. Lee, H. Park, and Y. Park, “High-resolution 3D refractive index tomography and 2D synthetic aperture imaging of live phytoplankton,” J. Opt. Soc. Korea **18**, 691–697 (2014). [CrossRef]

**11. **S. Johnsen, *The Optics of Life: A Biologist’s Guide to Light in Nature* (Princeton University, 2012).

**12. **S. Colin, L. P. Coelho, S. Sunagawa, C. Bowler, E. Karsenti, P. Bork, R. Pepperkok, and C. De Vargas, “Quantitative 3D-imaging for cell biology and ecology of environmental microbial eukaryotes,” Elife **6**, e26066 (2017). [CrossRef]

**13. **F. Lombard, E. Boss, A. M. Waite, J. Uitz, L. Stemmann, H. M. Sosik, J. Schulz, J.-B. Romagnan, M. Picheral, and J. Pearlman, “Globally consistent quantitative observations of planktonic ecosystems,” Front. Mar. Sci. **6**, 196 (2019). [CrossRef]

**14. **M. Picheral, L. Guidi, L. Stemmann, D. M. Karl, G. Iddaoud, and G. Gorsky, “The underwater vision profiler 5: an advanced instrument for high spatial resolution studies of particle size spectra and zooplankton,” Limnol. Oceanogr. Methods **8**, 462–473 (2010). [CrossRef]

**15. **R. J. Olson and H. M. Sosik, “A submersible imaging-in-flow instrument to analyze nano-and microplankton: imaging FlowCytobot,” Limnol. Oceanogr. Methods **5**, 195–203 (2007). [CrossRef]

**16. **J. S. Jaffe, “Underwater optical imaging: the past, the present, and the prospects,” IEEE J. Ocean. Eng. **40**, 683–700 (2014). [CrossRef]

**17. **A. D. Mullen, T. Treibitz, P. L. Roberts, and J. S. Jaffe, “An underwater microscope for *in situ* imaging of seafloor organism,” in *Novel Techniques in Microscopy* (OSA, 2017), paper NTu1C–1.

**18. **E. Acevedo-Trejos, G. Brandt, J. Bruggeman, and A. Merico, “Mechanisms shaping size structure and functional diversity of phytoplankton communities in the ocean,” Sci. Rep. **5**, 8918 (2015). [CrossRef]

**19. **H. Sun, P. W. Benzie, N. Burns, D. C. Hendry, M. A. Player, and J. Watson, “Underwater digital holography for studies of marine plankton,” Philos. Trans. R. Soc. A **366**, 1789–1806 (2008). [CrossRef]

**20. **E. J. Davies and R. Nepstad, “In situ characterisation of complex suspended particulates surrounding an active submarine tailings placement site in a Norwegian fjord,” Regional Stud. Mar. Sci. **16**, 198–207 (2017). [CrossRef]

**21. **J. Sheng, E. Malkiel, and J. Katz, “Digital holographic microscope for measuring three-dimensional particle distributions and motions,” Appl. Opt. **45**, 3893–3901 (2006). [CrossRef]

**22. **S. Samson, T. Hopkins, A. Remsen, L. Langebrake, T. Sutton, and J. Patten, “A system for high-resolution zooplankton imaging,” IEEE J. Ocean. Eng. **26**, 671–676 (2001). [CrossRef]

**23. **A. Zweifler, D. Akkaynak, and T. Mass, and T. Treibitz, “*In situ* analysis of coral recruits using fluorescence imaging,” Front. Mar. Sci. **4**, 273 (2017). [CrossRef]

**24. **D. Akkaynak, T. Treibitz, B. Xiao, U. A. Gürkan, J. J. Allen, U. Demirci, and R. T. Hanlon, “Use of commercial off-the-shelf digital cameras for scientific data acquisition and scene-specific color calibration,” J. Opt. Soc. Am. A **31**, 312–321 (2014). [CrossRef]

**25. **D. M. Kocak, F. R. Dalgleish, F. M. Caimi, and Y. Y. Schechner, “A focus on recent developments and trends in underwater imaging,” Mar. Technol. Soc. J. **42**, 52–67 (2008). [CrossRef]

**26. **A. Jordt-Sedlazeck and R. Koch, “Refractive structure-from-motion on underwater images,” in *IEEE International Conference on Computer Vision* (2013), pp. 57–64.

**27. **C. Ancuti, C. O. Ancuti, C. De Vleeschouwer, R. Garcia, and A. C. Bovik, “Multi-scale underwater descattering,” in *IEEE International Conference on Pattern Recognition* (IEEE, 2016), pp. 4202–4207.

**28. **M. Sheinin and Y. Y. Schechner, “The next best underwater view,” in *IEEE Conference on Computer Vision and Pattern Recognition (CVPR)* (2016), pp. 3764–3773.

**29. **H. Singh, J. G. Bellingham, F. Hover, S. Lemer, B. A. Moran, K. Von der Heydt, and D. Yoerger, “Docking for an autonomous ocean sampling network,” IEEE J. Ocean. Eng. **26**, 498–514 (2001). [CrossRef]

**30. **Y. Wang, S. Negahdaripour, and M. D. Aykin, “Calibration and 3D reconstruction of underwater objects with non-single-view projection model by structured light stereo imaging,” Appl. Opt. **55**, 6564–6575 (2016). [CrossRef]

**31. **T. Nakatani, S. Li, T. Ura, A. Bodenmann, and T. Sakamaki, “3D visual modeling of hydrothermal chimneys using a rotary laser scanning system,” in *IEEE Symposium on Underwater Technology and Workshop on Scientific Use of Submarine Cables and Related Technologies* (IEEE, 2011), pp. 1–5.

**32. **M. Leonardi, A. Stahl, M. Gazzea, M. Ludvigsen, I. Rist-Christensen, and S. M. Nornes, “Vision based obstacle avoidance and motion tracking for autonomous behaviors in underwater vehicles,” in *IEEE OCEANS* (IEEE, 2017), pp. 1–10.

**33. **E. Kelasidi, K. Y. Pettersen, J. T. Gravdahl, and P. Liljebäck, “Modeling of underwater snake robots,” in *IEEE International Conference on Robotics and Automation* (IEEE, 2014), pp. 4540–4547.

**34. **J. Gregson, M. Krimerman, M. B. Hullin, and W. Heidrich, “Stochastic tomography and its applications in 3D imaging of mixing fluids,” ACM Trans. Graph. **31**, 52 (2012). [CrossRef]

**35. **A. Vainiger, Y. Y. Schechner, T. Treibitz, A. Avni, and D. S. Timor, “Optical wide-field tomography of sediment resuspension,” Opt. Express **27**, A766–A778 (2019). [CrossRef]

**36. **A. Aides, A. Levis, V. Holodovsky, Y. Y. Schechner, D. Althausen, and A. Vainiger, “Distributed sky imaging radiometry and tomography,” in *IEEE International Conference on Computational Photography (ICCP)* (IEEE, 2020), pp. 1–12.

**37. **A. Levis, Y. Y. Schechner, A. B. Davis, and J. Loveridge, “Multi-view polarimetric scattering cloud tomography and retrieval of droplet size,” Remote Sens. **12**, 2831 (2020). [CrossRef]

**38. **A. Aides, Y. Y. Schechner, V. Holodovsky, M. J. Garay, and A. B. Davis, “Multi sky-view 3D aerosol distribution recovery,” Opt. Express **21**, 25820–25833 (2013). [CrossRef]

**39. **J. Frank, “Single-particle reconstruction of biological macromolecules in electron microscopy–30 years,” Q. Rev. Biophys. **42**, 139–158 (2009). [CrossRef]

**40. **R. R. Lederman and A. Singer, “A representation theory perspective on simultaneous alignment and classification,” Appl. Comput. Harmon. Anal. (2019). [CrossRef]

**41. **L. Wang, A. Singer, and Z. Wen, “Orientation determination of Cryo-EM images using least unsquared deviations,” SIAM J. Imaging Sci. **6**, 2450–2483 (2013). [CrossRef]

**42. **A. Levis, Y. Y. Schechner, and R. Talmon, “Statistical tomography of microscopic life,” in *IEEE Conference on Computer Vision and Pattern Recognition (CVPR)* (2018), pp. 6411–6420.

**43. **I. Greenberg and Y. Shkolnisky, “Common lines modeling for reference free Ab-initio reconstruction in cryo-EM,” J. Struct. Biol. **200**, 106–117 (2017). [CrossRef]

**44. **E. C. Orenstein, D. Ratelle, C. Briseño-Avena, M. L. Carter, P. J. Franks, J. S. Jaffe, and P. L. Roberts, “The Scripps Plankton Camera system: a framework and platform for in situ microscopy,” Limnol. Oceanogr. Methods **18**, 681–695 (2020). [CrossRef]

**45. **D. A. Agard, Y. Hiraoka, P. Shaw, and J. W. Sedat, “Fluorescence microscopy in three dimensions,” Methods Cell Biol. **30**, 353–377 (1989). [CrossRef]

**46. **D. A. Agard, “Optical sectioning microscopy: cellular architecture in three dimensions,” Ann. Rev. Biophys. Bioeng. **13**, 191–219 (1984). [CrossRef]

**47. **M. Alterman, Y. Y. Schechner, M. Vo, and S. G. Narasimhan, “Passive tomography of turbulence strength,” in *European Conference on Computer Vision (ECCV)* (2014), pp. 47–60.

**48. **V. Holodovsky, Y. Y. Schechner, A. Levin, A. Levis, and A. Aides, “In-situ multi-view multi-scattering stochastic tomography,” in *IEEE International Conference on Computational Photography (ICCP)* (2016), pp. 1–12.

**49. **D. Gürsoy, Y. P. Hong, K. He, K. Hujsak, S. Yoo, S. Chen, Y. Li, M. Ge, L. M. Miller, Y. S. Chu, and V. De Andrade, “Rapid alignment of nanotomography data using joint iterative reconstruction and reprojection,” Sci. Rep. **7**, 11818 (2017). [CrossRef]

**50. **I. Ihrke and M. Magnor, “Image-based tomographic reconstruction of flames,” in *ACM SIGGRAPH/Eurographics Symposium on Computer Animation* (Eurographics Association, 2004), pp. 365–373.

**51. **A. Levis, Y. Y. Schechner, and A. B. Davis, “Multiple-scattering microphysics tomography,” in *IEEE Conference on Computer Vision and Pattern Recognition (CVPR)* (2017), pp. 6740–6749.

**52. **D. Ren, C. Ophus, M. Chen, and L. Waller, “A multiple scattering algorithm for three dimensional phase contrast atomic electron tomography,” Ultramicroscopy **208**, 112860 (2020). [CrossRef]

**53. **B. D. Trifonov, “Tomographic reconstruction of transparent objects,” Ph.D. thesis (University of British Columbia, 2006).

**54. **A. Veeraraghavan, A. V. Genkin, S. Vitaladevuni, L. Scheffer, S. Xu, H. Hess, R. Fetter, M. Cantoni, G. Knott, and D. Chklovskii, “Increasing depth resolution of electron microscopy of neural circuits using sparse tomographic reconstruction,” in *IEEE Conference on Computer Vision and Pattern Recognition* (2010), pp. 1767–1774.

**55. **Y. Shkolnisky and A. Singer, “Viewing direction estimation in cryo-EM using synchronization,” SIAM J. Imaging Sci. **5**, 1088–1110 (2012). [CrossRef]

**56. **R. Gordon, R. Bender, and G. T. Herman, “Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography,” J. Theor. Biol. **29**, 471–481 (1970). [CrossRef]

**57. **W. van Aarle, W. J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. De Beenhouwer, K. J. Batenburg, and J. Sijbers, “Fast and flexible X-ray tomography using the ASTRA toolbox,” Opt. Express **24**, 25129–25147 (2016). [CrossRef]

**58. **Z. V. Finkel, J. Beardall, K. J. Flynn, A. Quigg, T. A. V. Rees, and J. A. Raven, “Phytoplankton in a changing world: cell size and elemental stoichiometry,” J. Plankton Res. **32**, 119–137 (2009). [CrossRef]

**59. **C. L. Lawson, M. L. Baker, C. Best, C. Bi, M. Dougherty, P. Feng, G. Van Ginkel, B. Devkota, I. Lagerstedt, S. J. Ludtke, and R. Newman, “EMDataBank.org: unified data resource for CryoEM,” Nucleic Acids Res. **39**, D456–D464 (2010). [CrossRef]

**60. **H. Sosik, E. E. Peacock, and E. Brownlee, “Annotated plankton images data set for developing and evaluating classification methods,” WHOAS2015, http://darchive.mblwhoilibrary.org/handle/1912/7341.