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).

 

Fig. 1. Specimens viewed in nature have random sizes and poses. Light passes through, yielding a single 2D image per specimen. The goal is to tomographically recover the 3D volumetric distribution of the plankton (on average) and estimate the population size distribution.

Download Full Size | PDF

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.

 

Fig. 2. (a) Illustration of the Scripps Plankton Camera (SPC). (b) An image of the SPC in situ. The dark-field microscope’s sample volume is between the ports.

Download Full Size | PDF

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,

 

Fig. 3. 3D Fourier domain and two central slices ${{{\tilde {\textbf y}}}_n},{{{\tilde {\textbf y}}}_m}$. Each central slice is perpendicular to the direction-vector of the respective projection in the microscope. The two slices intersect at a common tilt axis (common line) in the 3D Fourier domain.

Download Full Size | PDF

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

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

Using this formulation, we propose an algorithm for empirically estimating the rotations matrices $\{{{{\hat {\textbf R}}}_n}\} _{n = 1}^N$.

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

 

Fig. 4. Central slice in 3D Fourier domain (seen in Fig. 3) is equivalent to a 2D Fourier transform of an acquired image. A common line (also seen in Fig. 3) is at respective corresponding angles ${\psi _{n \to m}},{\psi _{m \to n}}$ in the 2D Fourier domains. Noisy data (e.g., due to non-rigid deformation of specimens) leads to erroneous estimation of these angles, ${\hat \psi _{n \to m}},{\hat \psi _{m \to n}}$.

Download Full Size | PDF

In Cryo-EM literature [43], a parametric model had been suggested to express the distribution of $\Delta {\psi _{nm}}$. The model dictates the following:

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

 

Fig. 5. Histogram of the angular error $\Delta {\psi _{nm}}$ (in degrees) with ${V_{{\max}}} = 9\%$, when data is affected by simulated random non-rigid deformations. The red curve is a fit of the histogram to Eq. (18). The fit has a coefficient of correlation 0.98.

Download Full Size | PDF

 

Fig. 6. Influence of the deformation coefficient ${V_{{\max}}}$ on the parameters of the statistical model [Eq. (18)]. The stronger the deformations, the lower is the probability that correlation-based angle pairs are indicative. Pairs that are indicative have a larger uncertainty as ${V_{{\max}}}$ increases.

Download Full Size | PDF

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:

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).

 

Fig. 7. Top left, the pre-computed PDF of a triplet score $s$, if the triplet of estimated relative rotations is arbitrary. Top right, the pre-computed PDF of a triplet score $s$, if the triplet of estimated relative rotations is indicative. Here estimated relative rotations have errors whose standard deviation is $\sigma = {6^\circ}$. Bottom, empirical histogram ${\cal H}(s)$ of triplet scores. It is modeled as a convex sum of the above PDFs, controlled by $P$. Fitting this parametric model to ${\cal H}(s)$ yields an estimate of ${{\hat P}},\hat \sigma$.

Download Full Size | PDF

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:

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}$:

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:

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,

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.

 

Fig. 8. (a) 3D volumetric object is randomly scaled, deformed shifted, and projected onto 100 random orientations. (b) Here we show several image samples from this simulation.

Download Full Size | PDF

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.

 

Fig. 9. Recovery using data described in Fig. 8. Top, visualization of the 3D object density estimated by (a) our method versus (b) that of Ref. [42]. (c) and (d) The respective scatter plots of $\hat \beta$ versus true $\beta$ (values drawn from 1% random voxels). The red line represents a perfect fit. Our approach clearly yields a better fit to ground-truth than [42]. (e) The statistical model [Eq. (18)] fitted to the empirical histogram of angular errors, similarly to Fig. 6. In the process described in Section 3.D, estimated angle pairs are better constrained. After a single iteration, our error model (red plot) has a narrower standard deviation and a higher indicative probability, compared to errors (blue plot) obtained in the method of Ref. [42].

Download Full Size | PDF

Table 1. Results Corresponding to Data Described in Fig. 8^a

View Table

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$.

 

Fig. 10. Experiment on a population of Daphnia sp., imaged using our system. Top, a few raw dark-field images. Middle, corresponding cross-validation re-projections, as described in Section 3.A, using our tomographically estimated object. Bottom, corresponding cross-validation re-projections, using the estimation resulting from the method of Ref. [42].

Download Full Size | PDF

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.

 

Fig. 11. Experiment on a population of Oithona sp., imaged using our system. Top, a few raw dark-field images. Middle, corresponding cross-validation re-projections of our tomographically estimated object. Bottom, corresponding cross-validation re-projections of the object resulting from the method of Ref. [42].

Download Full Size | PDF

 

Fig. 12. Volumetric density recovery: Oithona sp.

Download Full Size | PDF

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.

 

Fig. 13. Resulting size distribution histograms (blue) and fitted [Eq. (31)] density (red). (a) Oithona sp. The fit 95% confidence intervals are $\xi \in [6.48\;6.5]$, $\nu \in [0.138\;0.151]$. (b) Daphnia sp. The fit 95% confidence intervals are $\xi \in [1.12\;1.21]$, $\nu \in [0.244\;0.308]$.

Download Full Size | PDF

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.

 

Fig. 14. Volumetric density recovery: Pyramimonas longicauda.

Download Full Size | PDF

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.

 

Fig. 15. Experiment on a population of Pyramimonas longicauda, from dataset [4]. Top, a logarithm operation on raw bright-field images yields a linear model as in Eq. (1). A few samples are shown. Middle, corresponding cross-validation re-projections of our tomographically estimated object. Bottom, corresponding cross-validation re-projections of the object resulting from the method of Ref. [42].

Download Full Size | PDF

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.