1. Introduction

Molecular multiplexed and superplexed optical imaging techniques are emerging with potential to enhance exploration of basic biology and impact biomedical applications such as optical biopsy [1,2] and image-guided surgery [3–5]. Examples of recent advances in optical molecular multiplexed imaging include: multicolor single-molecule tracking super-resolution microscopy [6]; excitation-emission [7] and spatial-spectral multiplexing multiphoton microscopy [8,9]; simultaneous label-free autofluorescence-multiharmonic (SLAM) microscopy for contrast-rich optical biopsy [1] using only a single laser; deep tissue hyperspectral multiphoton microscopy using angle-tuned bandpass filters and three excitation wavelengths to image up to 10 different fluorophores simultaneously [10]; fast phasor-based hyperspectral microscopy that facilitates spectral Fourier phasor analysis in hardware [11]; whole-body hyperspectral fluorescence cryo-imaging [12]; multiplexed Raman endoscopy [13]; varied length Stokes shift fluorophores [14]; targeted, activatable probes [15,16]; 20-color stimulated Raman scattering microscopy using a carbon rainbow, "Carbow", with potential for supermultiplexing using barcoding; and, multicolor carbon nanotubes [17]. Hyperspectral detection helps to support many of these technologies by capturing multiple overlapping spectral signatures simultaneously using a contiguous array detector. This contrasts with traditional light microscopy for which dichroic mirrors and bandpass filters are often used to capture discrete signals with control experiments to verify that spectral cross-talk among the channels is insignificant. Hyperspectral approaches enrich contrast among intrinsic and extrinsic probes and isolate background signals for removal.

The rapid growth of this field towards the applications listed above is likely to push the limits of imaging depth and speed such that photon-starved regimes will become more prevalent than taking snapshots with long integration times. Regardless, a general challenge for hyperspectral imaging is to maximize the accuracy of spectral decomposition (spectral unmixing to quantify and display the individual spectral signatures) for any given imaging conditions. Here, we introduce a broadly applicable denoising concept and algorithm specific to hyperspectral microscopy in which a convolution filter is applied in the spectral domain rather than the spatial domain. This approach, termed spectral vector denoising, features denoising power that scales with the size of the data set and the number of observations of a sparse set of spectral signatures. Of course, traditional spatial convolution filters may be combined with this approach when the pixel spatial sampling exceeds the optical resolution for combined spectral-spatial denoising that preserves spatial information but with greater denoising power than either approach alone.

Spectral vector denoising maps each pixel’s discrete spectrum to a point in spectral space and utilizes spectral space coordinates to denoise the raw hyperspectral data cube. Here, we apply a discrete Fourier transform hyperspectral image analysis, developed by others [11,18–20], and we introduce discrete Chebyshev polynomial transform hyperspectral image analysis that may be used to translate pixel spectra to spectral space. Each pixel’s spectrum is recalculated by performing a Gaussian weighted average in which the Gaussian weight is determined by the pixel’s neighbors and their separation in spectral space. Thus, pixels with similar spectral shapes are averaged together, reducing noise in individual spectra. Spectral vector denoising is a simple hyperspectral denoising method that relies solely on spectral information, does not risk spatial oversmoothing, and preserves edges in hyperspectral images, in contrast to traditional methods for denoising microscopy images using spatial filters. Commonly, denoising methods for microscopy images are applied to each 2D channel of the hyperspectral image individually, i.e., the spectral channels are treated separately, and the spectral information is not utilized to denoise [21]. Beyond common methods in microscopy, state-of-the-art denoising methods developed for hyperspectral imaging in other fields, such as satellite imaging, include sparse representation and low-rank constraint utilizing local redundancy and correlation [22], singular value decomposition [23], principal component analysis [24], a spectral-spatial adaptive total variation model [21], and noise-adjusted iterative low-rank matrix approximation [25]. Most of these methods combine spatial and spectral information to reduce noise, and typically some spatial information is sacrificed. Spectral vector denoising does not sacrifice spatial information and is particularly useful for hyperspectral fluorescence microscopy, where denoising is applied to raw hyperspectral data prior to fluoresence spectral decomposition analyses (spectral unmixing) such as non-negative least squares fitting. Since there is no reliance on relative spatial information, multiple images of the same biological sample may be denoised as a group, which increases denoising power. Here, spectral vector denoising was first developed illustratively using a simulated hyperspectral image composed of five fluorescent signals using measured Alexa Fluor (AF) dye emission spectra with simulated noise. We then tested and characterized the denoising performance to extend the accuracy of spectral unmixing with imaging depth using hyperspectral multiphoton microscopy of AF labeled epithelial cancer cells embedded in a tissue phantom.

2. Methods

2.1 Spectral vector space coordinates

Inspired by Fourier phasor analyses for fluorescence lifetime imaging microscopy [19,20,26,27] and for fluorescence spectral imaging [11,18], here we additionally explore n-dimensional Chebyshev space for non-periodic spectral signatures in a similar fashion. The 2D Fourier phasor $z_f(n)$ (Eq. (1)) has been termed a "phasor" because it is derived from discrete Fourier transform calculations $G(n)$ (Eq. (2)) and $S(n)$ (Eq. (3)) [18,27,28] which yield an amplitude and a phase for the phasor $z_f(n)$. To extend a similar coordinate calculation to Chebyshev polynomials which do not naturally yield phasors, we suggest "spectral vectors" as an appropriate generalized terminology to include higher dimensional spaces resulting from discrete Chebyshev or Fourier transforms. The Chebyshev spectral vector space provides an analogous but distinct spread of the spectral signatures with complete separation of a set of Alexa Fluor dyes in 3 to 5-dimensional space (Figs. 1 & 3). Here, Chebyshev vector space outperforms Fourier space for separation of this dye set and denoising, however Fourier may outperform Chebyshev for alternative dye sets. We suggest Chebyshev vector analysis as a complementary approach. We note that other methods for calculating a basis for clustering pixels in a spectral space, like principle component analysis (PCA) [23,24], could be utilized to produce optimal separation of dye sets and methods like PCA would minimize the dimensionality of the vector space. Note that PCA requires discovery of the optimal basis for each set of chromophores or each image. Here, we chose to use Fourier and Chebyshev vector spaces because they provide a sufficient basis for fluorophore separation and because these transforms are independent of the fluorophore set used, allowing for broad applicability.

 

Fig. 1. (a) Delta functions at 10 nm intervals across the visible light spectrum are plotted in exemplary 2D Fourier phasor space vs. 2D and 3D Chebyshev vector spaces. Red corresponds to a delta function at 700 nm and violet to a delta function at 400 nm. (b) An intensity-dependent multiplicative noise model was used to add noise to groups of 16 channel hyperspectral pixels set to one of five AF emission spectra or one of five mixtures of various AF spectra. The ten noisy spectra groups were plotted in 2D and 3D Chebyshev vector spaces.

Download Full Size | PDF

The spectral vector denoising algorithm is performed on the raw hyperspectral image cube. The user may select to calculate spectral space coordinates using either Fourier or Chebyshev vectors. We first briefly review calculation of spectral phasor coordinates into the commonly used 2D Fourier phasor space and then introduce the analogous approach using Chebyshev vectors and higher dimensional spectral spaces. In 2D Fourier phasor space, Fourier phasors are calculated using Eq. (1). The 2D Fourier phasor $z_f(n)$ (Eq. (1)) is taken directly from Cutrale et. al.’s methods of Fourier phasor space analysis of fluorescence spectra [18]. $G(n)$ (Eq. (2)) and $S(n)$ (Eq. (3)) are normalized variations of the two components of the discrete Fourier transform of the spectrum, where $n$ corresponds to the $n^{th}$ harmonic of the Fourier series expansion, $x$ is the index of a hyperspectral channel, $I(x)$ is the intensity of the $x^{th}$ channel in the hyperspectral pixel, and $m$ is the number of channels in the hyperspectral image. Wavelength range of the input spectra is scaled from 0 to 2$\pi$ to span the phasor plot. Phasor coordinates are also scaled by the intensity of the pixel, so the phasor location reflects the shape of the spectrum and not its amplitude.

The 2D Chebyshev vector is plotted with coordinates $(U_n,U_{n+1})$ (Eq. (4)), where $m$ is the number of channels in the hyperspectral image and $T_n(x)$ are Chebyshev polynomials of the first kind. The $U_n$ (Eq. (4)) are calculated using a discrete Chebyshev transform of the spectrum using Chebyshev polynomials of the second kind, hence the $U_n$ notation. Wavelength range of the input spectra is scaled from −1 to 1 to span the spectral vector space. Just as for the familiar phasor coordinates, spectral vector coordinates are scaled by the intensity of the pixel, so the vector location reflects the shape of the spectrum and not its amplitude.

A 2D Fourier phasor plot for $z_f(2)$ and a Chebyshev vector plot for $(U_2,U_{3})$ vectors are shown in Fig. 1(a), where delta functions spanning the visible light spectrum have been translated to and plotted in the two spectral vector spaces. This figure illustrates that spectral peaks are spread differently throughout spectral space depending on the type of spectral vectors used. Extending spectral vector spaces to higher dimensions often has the advantage of further separating vector points corresponding to different spectral peaks and shapes. This concept is illustrated in Fig. 1(a), where there is increased separation between spectral delta functions in 3D Chebyshev vector space as compared to 2D Chebyshev vector space. In Fig. 1(b), a signal-dependent multiplicative noise model (see Supplemental Method A) was used to add noise to a 16-channel hyperspectral image composed of 2500 pixels set to one of 10 emission spectra: five were single AF dye emission spectra (AF633, AF647, AF660, AF680, and AF700) and five were various combinations of these five AF dyes in equal proportions, as identified in 1(b). The addition of noise to the AF spectra results in a cluster of points in spectral vector space with a centroid at the position of the original noise-free spectra. We can see in Fig. 1(b) that the AF633 and AF647 noisy clusters are highly overlapped in 2D Chebyshev vector space but well separated in 3D Chebyshev space. In Fourier phasor space, when two distinct spectra are summed and then translated to phasor space, the corresponding phasor will be the sum of the vectors corresonding to each of the summed spectra [18,20]. Each vector coordinate in both Fourier and Chebyshev spaces is calculated as a sum over the number of hyperspectral channels of a normalized intensity value for that channel that is scaled by a factor calculated from the channel number and order of the vector component. Thus similarly for Chebyshev vector space, when two spectra are summed, the corresponding vector in spectral space will be a sum of the two original spectral vectors. This math falls directly from the calculation of the vector coordinates. As expected, in Fig. 1(b) the spectral space centroids of the AF mixtures are positioned linearly between their component spectra and correspond to the relative amplitudes of each dye in the spectrum, just as for Fourier phasor space [18,20]. Fourier vectors were extended to $n$-dimensional space for even $n$ by using the $n$-dimensional coordinates $(G(2),S(2), G(3),S(3),\ldots,G(\frac {n}{2}+1),S(\frac {n}{2}+1))$. Chebyshev vectors were extended to $n$-dimensional space for any positive integer $n$ by using the $n$-dimensional coordinates ($U_2,U_3,\ldots,U_{n+1}$). Note that just like the Chebyshev vector spaces created here, one can use odd dimensionality to create Fourier vector spaces using mixtures of Fourier harmonic G and S components.

2.2 Spectral vector denoising algorithm

Our spectral vector denoising algorithm recalculates the spectra of each pixel in a hyperspectral image by performing a weighted average with the spectra of pixels that are its nearest neighbors in spectral vector space. First, all pixels in the image are translated into spectral vector space. For each pixel $p_{i}$, where $i$ is an index used to identify a single pixel, with corresponding spectral vector $\vec {v_{i}}$, a weighted average is calculated using the spectra of all pixels within 3$\sigma$ of $\vec {v_{i}}$ in spectral vector space, where $\sigma$ is a user-input Euclidean distance in spectral vector space that defines the width of the spectral space Gaussian filter. The Gaussian weight $w$ (Eq. (5)) used in the weighted average is a function of $x_{ij}$, the Euclidean distance in spectral space between spectral vectors $\vec {v_{i}}$ and $\vec {v_{j}}$.

The weighted average calculation is given in Eq. (6), where $p_{i}$’s discrete spectra is represented in vector form by $\vec {c_{i}}$, where the elements are the intensity values for each of the $m$ hyperspectral channels. $\vec {c_{i}'}$ yields the denoised spectra. The intensity of each pixel considered in the weighted average is normalized by the sum of its discrete spectra, $N_{i}$. $p_{i}$’s denoised spectra is renormalized to $p_{i}$’s original intensity.

The recalculated pixel spectra are placed into a new, denoised hyperspectral image cube. This process is illustrated in Fig. 2.

 

Fig. 2. Illustration of the spectral vector denoising algorithm. The discrete noisy spectra from each pixel in the raw hyperspectral image cube are translated to spectral vector space. Spectra are recalculated using a weighted average based on distance from the other pixels’ vector space coordinates. Denoised spectra are placed into the denoised hyperspectral image cube.

Download Full Size | PDF

Component spectra in hyperspectral fluorescence microscopy images were unmixed using the hyperspectral fluorescence unmixing software, HyperViewer, previously developed by us [29]. Each pixel in the input hyperspectral image has a discrete raw fluorescence spectrum from measurements, forming an image cube with two spatial dimensions and one spectral. A known "basis" set of fluorophores is input to HyperViewer. This basis set is the known set of fluorophores used in the sample being imaged. The fluorescence spectra of the basis set are then linearly combined to fit the measured spectra in each pixel using a GPU-accelerated fast non-negative least squares fit, yielding the proportions of each fluorophore that are present in each pixel. The spectrum is thus decomposed, or "unmixed", to quantify the amplitudes of each fluorophore’s emission spectra in each pixel. The original hyperspectral image cube is remapped to a data cube using the basis set of fluorophores as a basis instead of the wavelength-intensity data in the original input cube.

3. Results and discussion

3.1 Simulated image test

Spectral vector denoising was first tested on a simulated image with artificial noise added in order to quantify denoising effects with a known noise-free, ground truth version of the image (Fig. 3). Measured emission spectra of five Alexa Fluor (AF) dyes (AF633, AF647, AF660, AF680, and AF700) were used to create the simulated hyperspectral image with 16 wavelength channels. The image was gridded into 10 rectangles with various mixtures of AF dyes in each rectangle. The simulated image with no noise added was decomposed into basis fluorescence spectra using HyperViewer [29] and is shown in Fig. 3(a) with a breakdown of the AF dye mixtures in each rectangle of the image grid. The five AF dye spectra and a background signal spectra were used as the basis spectra to unmix the hyperspectral image. Each rectangle in the top row of the simulated image is composed of one AF dye. The colors assigned to each AF dye to create the composite "unmixed" image are indicated in the top row. Basis intensity line profiles were taken across the middle of each row of rectangles and are shown in the far right of Fig. 3(a). For the noise-free simulated image, the rectangles have been unmixed uniformly and with very high accuracy. In the center column of Fig. 3(a), the raw (pre-unmixing) noise-free image has been translated to 3D Chebyshev vector space. Each rectangle appears as a dot in spectral vector space, as the spectrum of every pixel in each rectangle is identical and thus translated to the same point in vector space. The five AF basis spectra are also identified in vector space. In the right-most column of Fig. 3(a), a basis intensity line profile is plotted vs. x pixel coordinate for both the top and bottom row. The basis intensity values for each pixel along a horizontal line drawn through the center of either the top or bottom row are plotted on the y-axis. The basis intensity line profile indicates close to 100% accuracy in unmixing results for the noise-free image in Fig. 3(a).

 

Fig. 3. (a) Unmixed image, 3D Chebyshev vector plot, and basis intensity line profiles shown for simulated image grid of various AF dyes. (b) Noisy image (super-Poissonian for > 49 photons/px). (c) Traditional spatial Gaussian filter (st. dev. 1.1 px) applied to image in (b). (d) 3D Chebyshev spectral vector denoising ($\sigma$ = 0.4) combined with a conservative spatial Gaussian filter (st. dev. 0.6 px) applied to the noisy image in (b).

Download Full Size | PDF

A signal intensity-dependent multiplicative noise model that accounts for shot noise and detector multiplicative noise [20] was used to test spectral vector denoising on the simulated image (see Supplemental Method A). Detector amplification common in photomultiplier tubes used in microscopy, among other sources of noise, results in super-Poissonian noise, meaning the variance of the noise is greater than mean number of photon counts per pixel. The noise model used in Fig. 3 is super-Poissonian for images with 49 photons per pixel or greater. The noisy simulated image is illustrated in Fig. 3(b) where its unmixed composite image, translation to 3D Chebyshev vector space, and basis intensity line profiles for each row in the image are shown. The added noise is evident in the unmixing results shown in the composite unmixed image as well as the basis intensity plots. The noise is also reflected in spectral space, where the single points from the noise-free spectral plot (Fig. 3(a)) have spread out in spectral space. A photon-starved example with simple shot noise is shown in Supplemental Fig. S1 for a mean of 10.98 photons per pixel.

Multiple denoising parameters were tested to maximize denoising without over-blurring in spectral space. User-chosen inputs to the denoising algorithm are dimensionality of spectral space and denoising $\sigma$ value that defines the width of the spectral space Gaussian filter. In Fig. 3, 3D Chebyshev vector denoising ($\sigma$ = 0.4) combined with a small spatial Gaussian filter (standard deviation of 0.6 pixels) (Fig. 3(d)) is compared to traditional denoising with a simple spatial Gaussian filter with a larger width (standard devation of 1.1 pixels) (Fig. 3(c)). For Chebyshev vector denoising, we chose the dimensionality and $\sigma$ value that maximized the peak signal to noise ratio of the denoised image as compared to the no-noise simulated image. The larger spatial Gaussian filter applied to the noisy image from Fig. 3(b) is shown in Fig. 3(c). The resulting unmixed images and basis intensity line profiles indicate that the noisy image has become more uniform and less splotchy, but the edges between rectangles have been blurred. The 3D vector denoising combined with a conservative spatial filter is shown in Fig. 3(d). The combined spectral-spatial method results in enhanced denoising compared to the spatial Gaussian filter alone, as well as a greater degree of edge preservation. In Supplemental Fig. S2, we demonstrate the utility of spectral vector denoising using an additional simulated image featuring a continuous distribution of dyes across the image, in contrast to the discrete mixtures used in Fig. 3. Note that the denoising power of the spectral vector denoising algorithm increases with discrete clustering of fluorescent signal groups in spectral space. Therefore, the algorithm may be most useful for experimental conditions in which clustering of fluorescent signals is expected.

3.2 Tissue phantom test

Multi-dimensional spectral vector denoising was also tested on real images of a tissue phantom [30]. Ovarian cancer Ovcar5 cells were tagged with various AF dyes and suspended in a material with similar diffractive properties to tissue. Hyperspectral images were taken of two tissue phantoms at a range of depths, allowing us to assemble 3D hyperspectral images. As depth increases, signal is reduced due to the diffractive properties of the suspension material, mimicking imaging at depth in tissue. One tissue phantom was composed of 6 types of fluorescently labelled Ovcar5 cells, each tagged with a different AF dye, and imaged using a 790 nm laser (Fig. 4(a)). For images of this phantom unmixed using HyperViewer, each cell should be unmixed to a single dye, as each cell was tagged with only one dye. Cells incorrectly unmixed to a mixture of dyes are a result of noise. Thus denoising is expected to increase the uniformity of pixel unmixing to a single dye constituent per pixel, and to increase the percentage of signal correctly unmixed to one specific dye within each cell. The other tissue phantom was composed of 2 types of fluorescently labelled Ovcar5 cells, both types tagged with a mixture of AF546 and AF568 but each with a different ratio of the two dyes, and imaged using a 1090 nm laser (Figs. 4(d), 4(e)). For images of this phantom, each cell should be unmixed to one of the two ratios of AF546 to AF568 dyes. Cells unmixed to AF dyes not present in the phantom or cells non-uniformly unmixed are again due to noise, and denoising is also expected to increase the uniformity of pixel unmixing within these cells to the correct dye ratio.

 

Fig. 4. (a) XZ-projection of denoised and unmixed images of a tissue phantom at various depths for which Ovcar5 cells are tagged with one of 6 AF dyes. Scale bar 50 $\mu$m. (b) Individual cells from the unmixed tissue phantom in (a) are shown before and after 5D Chebyshev vector denoising. Scale bar 10 $\mu$m. (c) Percent reduction in unmixing error after denoising vs. Z-depth ($\mu$m) is plotted for the images in (a). (d) XZ-projection of denoised and unmixed images of a tissue phantom at various depths for which Ovcar5 cells are tagged with two AF dyes, AF546 and AF568, in one of two ratios shown in (d). Scale bar 50 $\mu$m. (e) Individual cells from the unmixed tissue phantom in (d) are shown before and after denoising. Scale bar 10 $\mu$m. (e) Percent reduction in unmixing error after denoising vs. Z-depth ($\mu$m) is plotted for the image in (d).

Download Full Size | PDF

In Fig. 4, 5D Chebyshev vector denoising was used to denoise the tissue phantom images. We chose the dimensionality and $\sigma$ value that maximized unmixing accuracy, calculated based on the known tissue phantom being used. Projections along the x plane of two of the tissue phantom samples are shown in Fig. 4(a) and 4(d). Individual cells from the z-stack slices are picked out in Figs. 4(b) and 4(c) and each of the basis images are shown to illustrate the unmixing results and improvement in unmixing accuracy post-denoising. For example in the right cell highlighted in Fig. 4(b), the correct dye present in the cell is AF594, but for the raw unmixed image there is signal present in the basis channels for AF568, AF610, and AF633, indicated some pixels have been incorrectly unmixed. After denoising, signal in these incorrect channels is reduced and signal in the correct channel is increased. Additionally, edge preservation is observed. In the raw image, a clear edge of the cell is observed on the pixel level. After denoising, the exact same pixel edge is maintained with no blurring at the edges of the cell. Spectral vector denoising improves unmixing accuracy with no spatial blurring. This concept is similarly illustrated in the right cell highlighted in Fig. 4(b), where the correct unmixing channel for the two cells is AF610. In Fig. 4(e), one cell is highlighted for each of the two dye mix ratios. "Mix 1" contains a higher proportion of AF546 to AF568 and an example cell is shown on the left of Fig. 4(e). Here, there is signal present in each of the 4 incorrect channels, and the AF568 channel looks particularly noisy and heterogeneous, though the cell should have relatively uniform dye concentration. After denoising, signal in the incorrect unmixing channels is reduced and signal in the correct channels is increased and appears more uniform. The edge of the cells is the same in the raw and denoised image, illustrating edge preservation. "Mix 2" contains a more even ratio of AF546 to AF568 and an example cell is shown on the right of Fig. 4(e). Again, incorrect signal is observed in the basis channels for the raw image and unmixing error is reduced in the denoised image while maintaining the edges of the cell. For each of these cells, spectral vector denoising reduces the amount of signal in incorrect unmixing channels and increases signal and uniformity in the correct unmixing channels, thus reducing unmixing error. All spatial information is preserved. Percent reduction in unmixing error is plotted vs. z-depth in Figs. 4(c) and 4(f) for the images in Figs. 4(a) and 4(d) respectively (see Supplemental Fig. S3 for additional data). The image in Fig. 4(a) demonstrates up to 40% reduction in unmixing error with denoising (Fig. 4(c)). Reduction in error increases with depth reflecting the increase in unmixing error accompanying less signal at greater imaging depth. The image in Fig. 4(d) demonstrates up to 70% reduction in unmixing error with denoising (Fig. 4(f)).

4. Conclusion

We have developed a spectral vector denoising program for hyperspectral images that reduces spectral noise without sacrificing spatial information. This new denoising method builds on phasor spectral space analyses that enable facile visualization of complex hyperspectral imaging data through dimension reduction and clustering of pixels based on spectral content. Like phasor analysis, spectral vector denoising translates the discrete spectrum of each pixel in a hyperspectral image to multidimensional spectral space using a discrete Fourier or Chebyshev transform. Denoising is then performed using the spectral space coordinates to calculate a weighted average of pixel spectra and reduce spectral noise. The effectiveness of spectral vector denoising was demonstrated by testing the algorithm on both a simulated image with artificial noise and a real image of a tissue phantom. Both tests indicated that spectral vector denoising is an effective method for reducing noise in hyperspectral images obtained using fluorescence microscopy (Figs. 3 & 4). We anticipate that this algorithm may be useful for other types of hyperspectral optical images, although we have only evaluated the algorithm for use in hyperspectral fluorescence microscopy. A key advantage of this method is that the spectral dimension is utilized in order to minimize the loss of spatial resolution when denoising hyperspectral microscopy image cubes.