2.1.

Multichannel Mosaicking Algorithm

The presented algorithm applies pairwise normalized cross-correlation²¹ (NCC) as a basis for frame registration. Following Bedard et al.,⁵ cross-correlation maps are calculated between consecutive image frames [Figs. 1(a) and 1(b)] and the maps are normalized to an in-place autocorrelation (i.e., the theoretical maximum correlation). Correlation space then maps the probability of alignment across the spatial domain spanned by the intersection of two images being aligned, from which the location of maximum likelihood is selected as the “true” alignment location for the subsequent stitching step. We calculate the spatial cross-correlation map separately for each additional channel and average the maps in correlation space, producing a single map [Figs. 1(c) and 1(d)].

Fig. 1

Exemplary cross-correlation maps demonstrate improved noise tolerance for multichannel versus single-channel micromosaicking. (a) Single-channel NCC map shows typical landscape with global maximum representing location of the correct image alignment. (b) Gaussian noise (${\sigma }^2=0.02$) added to the images correlated in (a) results in fluctuations in the NCC map that obscure the location of the correct alignment. (c) NCC landscape averaged over 10 spectral channels reduces the probability of an incorrect alignment. (d) The location of the true alignment is not obscured after Gaussian noise [${\sigma }^2=0.14$, a 7-fold increase from (b)] is added to the 10-channel NCC configuration in (c), demonstrating robust noise tolerance compared to the single-channel analog. Image insets provide a visual juxtaposition of noise strength.

The code enabling this algorithm was developed in MATLAB R2019a (The MathWorks) in a Linux environment (Ubuntu 18.04LTS) with an Intel i7 6800K 3.40 GHz CPU and 32 GB DDR4 2400 MT/s RAM. The substantial computational cost of calculating multiple cross-correlations per frame-pair was alleviated by using MATLAB’s Parallel Computing Toolbox v7.0 for GPU interfacing. Speed benchmarks were carried out using a GeForce GTX 1080 Ti graphics card (NVIDIA). The code implementing this algorithm is available at the GitHub repository https://github.com/springlabnu/Multi-Channel-Mosaicking.

2.2.

Simulated Data Generation and Mosaicking Performance

Quantitative analysis and comparison to equivalent single-channel methods was carried out using a simulated video-rate data set derived from images obtained by empirical microscopy. We acquired a hyperspectral image cube (26 spectral bands, images shown in Fig S1 in the Supplemental Materials) of fixed Madin–Darby canine kidney cells stained with three dyes for visualization of the nucleus (Hoechst 33342), cytoskeleton (Alexa Fluor 488-Phalloidin), and mitochondria (MitoTracker Red CMXRos) using a laser scanning confocal microscope (Olympus FV3000). An artificial hyperspectral video stream was then generated from this data by cropping the field-of-view and traversing the image in a continuous raster scan with 94% area overlap between pairs of images [raster scan pattern illustrated in Fig. S2(a) in the Supplemental Materials]. This created a sequence of 65 hyperspectral image cubes ($512\times 512\text{pixels}$) with a known ground truth alignment for quantitative algorithm characterization. The image cubes were also linearly unmixed into their three molecular constituents, creating three unmixed channels that could also be used for comparative analysis.

Having the ground truth location for each frame means that, for each trial, the alignments may be checked against the lookup table of correct locations without the need to visualize the full mosaic. From here, we quantified mosaicking accuracy $\eta$ as the fraction of the 64 frame pairs (one less than the number of frames) that were correctly aligned compared to the ground truth (i.e., the original image). The alignment was counted as correct only for an exact match. Zero-mean additive Gaussian noise was added to each image in every frame, in order to degrade image quality and to mimic challenging experiments that could potentially cause alignment mistakes. In order to quantify noise power of the simulated images, we use the peak SNR, defined as

The mosaicking accuracy $\eta$ is approximately equal to the probability that a single pair of frames are correctly aligned [$P(\mathrm{correct})$]; that is, correct alignment is a Bernoulli variable with probability $\eta$ for which the amount of noise is a predictor. The probability of Bernoulli variables is given by a logistic function;²² this allows us to connect the fraction of correct alignments ($\eta$) to the noise strength causing the misalignments (PSNR) by

Testing of the present algorithm was carried out using $n=[\mathrm{2,26}]$ raw spectral channels. The unmixed channels were only analyzed by calculating their $D$ score to compare their predicted performance. For each $n$ value, $D$ was calculated for every permutation of $n$ channels, from which the permutation with the maximal $D$ value was selected as the optimal subset. As a control, a permutation was selected subjectively by an operator based on prioritizing channels with the brightest signal corresponding to each molecular constituent, intended to mimic selection of channels without knowledge of $D$. Histograms of all possible $D$ values (all possible permutations) are shown for these cases in Fig. S3 in the Supplemental Materials to demonstrate the optimization versus the user-selected permutations. Additive Gaussian noise with ${\sigma }^2=[0.003,1.000]$ was added to each channel in each image within the configuration, after which the algorithm was run five times, with $\eta$ being saved after each trial. The time for each trial to complete the full mosaic was also saved. The applied variance was translated to PSNR using Eq. (1), fit to Eq. (2), from which the ${\sigma }_{50}^2$ was extracted. Independent two-sample $t$-tests (Welch’s $t$-test) were calculated to find the significance of the improvements in ${\sigma }_{50}^2$ through optimization of subset permutation, as well as other comparisons.

These experiments were repeated for multiplicative Gaussian and Poisson noise. For these noise types, the amount of noise power delivered depends on signal, so the ${\sigma }_{50}^2$ value extracted is signal dependent as well. The characteristic performance is thus amended to ${\sigma }_{50}^2/I^2$ (i.e., normalized by the intensity) accordingly for these noise types.

2.3.

Empirical Data Generation

While the simulated data analysis thoroughly examines and characterizes the present methods, the ultimate test of the multichannel algorithm and the predictive power of $D$ is to verify that they produce useful images from an empirical data set subject to real-data constraints (Fig. 2). To create an empirical data set, ovarian cancer cells (POWDER cell line, Cellaria Biosciences) were plated and stained with two fluorescent dye stains (Hoechst 33342 and DiI) and two fluorescent dye–antibody conjugate markers (PerCP-CD44 and AF647-EGFR). Using a confocal fluorescent microscope with a 32-channel hyperspectral detector (Zeiss LSM880, channel bands shown in Table S1 in the Supplemental Materials) 30 frames were generated by taking separate images after manual movements of the stage, generating a set of frames with a wide range of overlap amounts (estimated at 77% to 96% based on alignments reported by the algorithm). Mosaicking was carried out for a 10-channel configuration, whose channel constituents were chosen as the set that maximizes $D$ (i.e., the optimal 10-channel subset). This was compared to a single-channel control representing the analogous single-channel method. The same area was imaged using a tiling protocol available in the microscope’s software (ZEN Blue 2.6) to create a gold standard image for qualitative comparison of the mosaics.

Fig. 2

For real data applications, ground truth knowledge is unavailable, requiring users to apply the principles shown here to improve mosaicking accuracy under noise loads without quantitative feedback. Dimensionality scores ($D$) must be calculated from sample frames, which is shown to have some effect on subset optimization (see Fig. S4 in the Supplemental Materials). In addition, the cutoff value $ϵ$ must be estimated using a histogram rather than determined quantitatively. If the channels calculated to maximize $D$ are not sufficient to prevent unacceptable mosaicking errors, users may include additional channels or attempt to fine-tune the estimation of $ϵ$ until results are satisfactory. Note that the $\oplus$ symbol indicates a logical “or.”

3.1.

Dimensionality Score

In an effort to predict the optimal subset of channels for multichannel mosaicking, we developed a new metric termed “dimensionality score” that quantifies the degree of independence of a given grouping of data channels, and by extension, predicts the characteristic noise tolerance. The optimal channel set, then, is predicted to be that which maximizes this score (thereby maximizing noise tolerance). This metric must capture the primary means by which the present multichannel algorithm suppresses the effects of noise to retain mosaicking accuracy.

The multichannel integrated NCC map (the mean NCC map; Fig. 1) suggests that the present approach offers two means of noise suppression. Adding noise to the image channels generates similar noise levels but with uncorrelated patterns in the individual correlation maps, and comparison of Figs. 1(b) and 1(d) shows that the variance in NCC space is reduced by taking the mean of multiple NCC maps, since random noise fluctuations area averaged out and reduced over several spectral channels. However, comparison of Figs. 1(a) and 1(c) reveals that this process also populates (sharpens) the central true peak feature while reducing the magnitudes of false peaks. Hence, multichannel mosaicking is more immune to noise compared to the single-channel case not only through averaging out random noise in NCC space but also because this spectral averaging apparently suppresses the false NCC peaks due to separate biological objects that share similar structural features, and therefore enhances the strength of the true peak versus false peaks. The strength of this suppression of false peaks is dependent on averaging over diverse NCC landscapes. Averaging multiple NCC maps from weakly correlated or uncorrelated channels causes true peak enhancement by utilizing a set of independent NCC maps with the same true peak location while stochastic false peaks generally are not amplified and tend to average to zero.

By this reasoning, the noise tolerance of a given set of channels depends primarily on the number of channels and the degree of dependence between those channels. The coefficient of correlation is already a well-known descriptor of the degree of dependence between two samples,²³ so it forms a natural basis for constructing such a descriptor. Since the samples here are images, we can calculate the two-dimensional (2-D) coefficient of correlation $R_{2\mathrm{D}}$, which is defined for 2-D matrices (images) $A$ and $B$ as

The selection of $ϵ$ has a significant effect on the performance of $D$ as a predictor for ${\sigma }_{50}^2$. An appropriate value of $ϵ$ realizes the separation of $R_{ij}$ values characteristic of channels that image nearly identical features. For a matrix of $R_{ij}$ values, this effect results in a characteristic division of the $R_{ij}$ matrix into distinct neighborhoods [Fig. 3(a)]. We demonstrate this effect for $ϵ=0.89$, leading to three such areas that likely correspond to the three molecular constituents imaged in this dataset. Without the knowledge of a ground truth for grading performance, an appropriate value of $ϵ$ may be estimated (Fig. 2) through analysis of a histogram of $R_{ij}$ [Fig. 3(b)] with the intention of separating the population of strongly correlated channel pairings from the rest. Here, with knowledge of the ground truth, the $ϵ$ value that provides the greatest predictive strength may be found by maximizing the Pearson correlation between $D$ and ${\sigma }_{50}^2$ (the characteristic noise tolerance), with a peak value at $ϵ=0.89$ giving a squared Pearson correlation of $R_p^2=0.97$ [Fig. 3(c)]. However, our results indicate that a considerable range [shaded area in Figs. 3(b) and 3(c)] of $ϵ$ values correspond to a central peak of Fig. 3(c), and therefore retain near-optimal predictive performance. Therefore, we expect that an estimation of $ϵ$ using the aforementioned histogram method will still ensure that $D$ remains a reliable and useful predictor of mosaicking performance. Note that the limits of $ϵ\to 0$ and $ϵ\to 1$ both reduce $D$ to $1+(1-\overline{R})*n$, which is an intuitive descriptor of dimensionality; simply the mean of one minus cross-channel correlation, scaled by $n$. However, in this limit, there is only one operative term ($\overline{R}$). Therefore, it is apparent that both noise suppression modes are not captured in this limit, resulting in a weakened (lower Pearson’s correlation) prediction of noise tolerance.

Fig. 3

Selection of the cross-channel correlation cutoff, $ϵ$, affects the predictive accuracy of $D$ by changing the strictness of “independence” definition. (a) A heatmap of the $R_{ij}$ matrix shows $ϵ=0.89$ (dotted line) segregating strongly correlated channel pairings into three major areas corresponding to the three molecular constituents imaged. (b) A histogram of $R_{ij}$ values shows that $ϵ$ values within the shaded area separate a population of strongly correlated channel pairings ($R_{ij}$ values) from the rest (dotted line indicates $ϵ=0.89$). (c) For data with known ground truth, predictive power may be optimized directly by maximizing Pearson’s correlation ($R_p^2$) of $D$ and ${\sigma }_{50}^2$ versus dimension score cutoff ($ϵ$). This reveals an optimal $ϵ=0.89$ corresponding to $R_p^2=0.97$; the shaded area from (b) now outlines a central peak, indicating that subprime $ϵ$ values estimated from (b) will still deliver near-optimal predictive strength.

Even with our parallelized implementation, calculating the cross-correlation remains the dominant computational bottleneck, and the mosaicking computation time thus linearly scales with the number of channels used. $D$ is employed to optimize the selection of a high-priority subset of channels in an effort to lessen this computational load. The optimal configuration is found by maximizing $D$ over channel selection choices that are sufficiently speedy (i.e., limiting the number of channels to reach a target mosaicking rate). This enables $D$ to be a practical tool for optimizing the selection of subsets of available channels. A brute-force calculation of $D$ scores for all possible permutations of the desired configuration (of which there may be hundreds or thousands) can be computed in a reasonable time frame, providing comparative values that predict each permutation’s noise tolerance. The problem of selecting the optimally performing subset of channels for mosaicking is, to the best of our knowledge, otherwise intractable, except by empirical trial-and-error by the operator combined with a heuristic quality criterion.

3.2.

Quantitative Mosaicking Performance

Increasing the number of spectral channels improves the mosaicking accuracy under noise load, as indicated by plotting $\eta$ versus PSNR [Fig. 4(a); see also Eq. (2)]. For a given PSNR, more channels always perform better (higher $\eta$), and at a given $\eta$, having more channels allows us to tolerate stronger noise levels (see Table S3 in the Supplemental Materials). These results are replicated for multiplicative Gaussian noise and Poisson noise, as shown in Fig. S5 in the Supplemental Materials. Therefore, the multichannel algorithm is expected to achieve greater accuracy than equivalent single-channel methods for noise models relevant in microscopy. Two permutations were tested for each configuration (number of channels): the optimal permutation predicted by $D$ and an operator-selected control (these permutations are shown in Table S2 in the Supplemental Materials). The control was assembled by selecting subsets of channels that evenly sampled the three molecular constituents, intended to represent an operator’s educated guess at the optimal configuration of channels (the alternative to using the dimensionality score for optimization). In all cases, the optimized permutation was superior (lower ${\mathrm{PSNR}}_{50}$) to the control [Fig. 4(b)].

Fig. 4

Quantification and characterization of mosaicking accuracy, dimensionality score, and computational speed under noisy imaging conditions. (a) Mosaicking accuracy under noise load is modeled using a logistic curve [Eq. (2)] where the logistic midpoint (dashed line) characterizes the noise tolerance. Performance is improved both by adding more channels and optimizing (opt.) the subset of channels utilized versus user-selected channels (ctrl.) based on even sampling of molecular constituents. (b) ${\sigma }_{50}^2$ values extracted from midpoints of (a) using Eq. (1) display improved noise tolerance through selection of a channel subset that maximizes $D$; asterisks represent levels of significance (* signifies $p<0.05$, **** signifies $p<0.0001$). (c) $D$ correlates with ${\sigma }_{50}^2$ for optimized permutations ($R^2=0.95$); data labels denote number of channels used. (d) Mosaicking speed decreases with increasing dimensionality score as a power law. Predicted performance of optimized channel sets is better than that of the user-selected control sets, and unmixed channels provide the best predicted performance. Dotted line denotes video-rate (15 fps) threshold.

The correlation and thus predictive strength remains high between the ${\sigma }_{50}^2$ values and $D$ [Fig. 4(c), $R_p^2=0.95$] for the optimized subsets. This is very similar to the previous correlation value of 0.97 between ${\sigma }_{50}^2$ values and $D$ for the operator control channel sets. Therefore, $D$ retains its prediction strength among the various subsets of optimized permutations. In all cases tested, the optimized permutations provide superior performance over the operator-selected ones [Fig. 4(d); at least $p<0.05$]. While the improvement is most striking for small numbers of channels, the effect is present for all cases, confirming an increase in noise tolerance by use of the dimensionality score to predict the optimally performing subset of channels.

The dimensionality score $D$ plotted against the mosaicking speed [Fig. 4(d)] demonstrates the trade-off between computational speed (increasing $y$) and dimensionality score (increasing $x$) for both control and optimized permutations of the raw channels, as well as the optimized permutations of unmixed channels. These plots were fit to power law curves, and two-tailed $t$-tests were evaluated for the coefficients of these fits to show that: (1) unmixed channels outperform optimized raw channels ($p<0.001$) and (2) optimized raw channels outperform operator-selected channels ($p<0.001$). Here, higher performance is denoted by a more elevated trendline, delivering faster performance and forecasting better noise tolerance.

3.3.

Empirical Mosaicking

Based on these results, the multichannel algorithm and optimization of channel selection should enable greater accuracy for real-data mosaicking. Section 2.3 outlines the generation of an empirical data set (without the privilege of ground truth) for this purpose. This empirical set consists of a series of images with varying interframe overlap of a primary tumor cell culture stained using four spectrally distinct fluorescent probes and captured by 32-channel microscopy using a precision stage. We followed the previously discussed steps to optimize the choice of channels for video-rate mosaicking (i.e., limited to up to 10 channels to retain video-rate computation). We estimated $ϵ$ using a histogram [Fig. 5(a)] of correlation values ($R_{ij}$) from the first video frame, finding that channels 1, 2, 12, 13, 21, 22, 26, 27, 28, and 32 maximized $D$. The images were analyzed by examining the unmixed composite images created by unmixing the raw data into the four endmembers (dyes) and creating a color map to represent all four dyes in one image. The Zeiss tiling mentioned earlier [Fig. 5(b)] provides a gold standard for qualitative comparisons of mosaicking accuracy and for evaluating the presence of artifacts.

Fig. 5

Micromosaicking of an empirical hyperspectral data set (described in Sec. 2.3). (a) Histogram of correlation values for 32 hyperspectral channels listed in Table S1 in the Supplemental Materials reveals that a choice of $ϵ=0.75$ (dotted line) separates a population of strongly correlated channels, informing calculation of $D$ [Eq. (4)]. (b) Tiling of four unmixed dye channels (annotated on image) serves as a qualitative gold standard to compare mosaics. (c) One-channel mosaicking produces severe artifacts such as repeated features (red circle), sharp edges (red rectangle), and a notably different aspect ratio compared to (b) despite being displayed at the same scale. These are indicative of significant misalignments and this configuration ultimately fails to provide an accurate mosaic. (d) Micromosaicking using the optimal 10 channels results in an image free of noticeable artifacts; qualitative comparison to the gold standard (b) indicates that the images have only minute differences (Video 1, 4.93 MB, MP4 [URL: https://doi.org/10.1117/1.JBO.24.12.126002.1]).

The single channel method produced an image [Fig. 5(c)] with many severe artifacts, including sharp edge artifacts and clear examples of feature duplication. These errors are presumed to result from alignment mistakes caused by noise distorting similar areas between frame pairs, ultimately hindering the NCC alignment step. The 10-channel optimized mosaicking [Fig. 5(d)] qualitatively outperformed the one-channel control, producing an image largely free of artifacts, substantially mimicking the gold standard. Additionally, step-by-step comparisons of the video feed and the associated intermediate mosaic can be visualized in Video 1.

This empirical test demonstrates that significant improvement was achieved in a real-data scenario compared to the single-channel analog. Cross-correlation-based algorithms have been implemented successfully for handheld microendoscopes in the literature,⁵ so the methods here should achieve similar results for multichannel versions of these probes. However, the primary limitation of the cross-correlation-based method implemented here is that alignments are limited to translational corrections, and therefore nontranslational motion must be negligible. To overcome this limitation, micromosaicking algorithms have been developed far beyond single-channel correlations and now typically utilize feature-based alignment strategies that allow for more complex frame-pair alignments such as affine transformations⁷ and local warping corrections.⁶ The success of the multichannel NCC method here motivates integration of these advanced methods with the present multichannel mosaicking algorithm in a manner that similarly takes advantage of the unique information available in multispectral and hyperspectral data. A limitation of this method is that smooth image layering such as graph-cut-based stitching^7,24 or alpha blending of overlapping areas⁴ have not yet been integrated. The present algorithm also does not compensate for compounding errors (or drift error), where small misalignments may build up due to the pairwise image alignment method and cause a significant discrepancy after many alignments, as it cannot correct itself after misalignments, as demonstrated in Fig. S2(b) in the Supplemental Materials. This problem has been addressed in earlier publications⁶ by referencing a global alignment system and refining pairwise alignments to fix these issues. Here, the results motivate the integration of such methods with the present algorithm in future work.