1. Introduction

Multiphoton microscopy (MPM) is a powerful and widely utilized technique for non-invasive and label-free imaging of biological tissues [1–8]. Multiple contrasts can be acquired by MPM using different nonlinear optical effects. For example, in two-photon microscopy (2PM), two-photon-excitation-fluorescence (2PEF) can be acquired from intrinsic fluorophores such as nicotinamide adenine dinucleotide hydrogen (NADH) in cells, and second harmonic generation (SHG) can be acquired from non-central symmetric molecules such as collagen fibers; in three-photon microscopy (3PM), three-photon-excitation-fluorescence (3PEF) can also be acquired from intrinsic fluorophores, and third harmonic generation (THG) can be acquired from interfaces such as the lipid aqueous medium interface [9–15]. Being able to acquire multiple contrasts is important in label-free imaging as it can provide a comprehensive understanding of tissue morphology and function. With the advancement in lasers as well as signal detection techniques, combining 2PM and 3PM in a single system has gained a lot of interest as it can provide complementary information about tissues including cells, collagen fibers, lipids, etc., and form information-rich images that are essential for label-free tissue imaging [16–21].

Several groups have demonstrated combined 2PM and 3PM imaging systems. F. Akhoundi et al. demonstrated a 2PM (SHG) and 3PM (3PEF and THG) probe using a single excitation wavelength at 1700nm [22]. However, their THG and 3PEF signals were acquired sequentially by a single photomultiplier tube (PMT), and samples were hematoxylin and eosin (H&E) stained for 3PEF imaging. At 1700-nm excitation wavelength, it is not possible to excite 2PEF from tissues without staining, and 3PEF is a higher-order nonlinear effect with much lower efficiency than 2PEF. Although excitation by a single wavelength is simple, it has many limitations in combined 2PM and 3PM imaging. On one hand, lasers with long wavelengths (>1 µm) cannot efficiently excite 2PEF from intrinsic fluorophores in tissues [23,24]. On the other hand, a short excitation wavelength (<1 µm) is not optimal for THG imaging as the signal wavelength will fall inside the ultraviolet regime [23].

To achieve both 2PM and 3PM from tissues label-free, several other groups including our group have used two excitation wavelengths, where a shorter wavelength (<1 µm) is for 2PM and a longer wavelength (>1 µm) for 3PM, respectively [16,17,24–27]. D. Septier et al. demonstrated 2PM (SHG) and 3PM (3PEF and THG) in mouse skin excited by 920 and 1300 nm, respectively [26]. By using a single PMT, all three images were acquired one after another. F. Aptel et al. performed 2PM (2PEF and SHG) and 3PM (THG) imaging of human cornea using dual excitation wavelengths [27], where a 730-nm wavelength excited 2PEF and a 1200-nm wavelength excited SHG and THG. The two wavelengths were applied on tissue in sequential mode. Previously, our group has demonstrated combined 2PM (2PEF and SHG) and 3PM (THG) using dual excitation wavelengths at 790 and 1580 nm obtained from a frequency-doubled erbium-doped fiber laser [16,17,28]. 3D volumetric 2PM and 3PM image stacks were acquired in sequential mode, where the 790- and 1580-nm excitation wavelengths were applied to tissue one after another. There are several drawbacks to the sequential acquisition of 2PM and 3PM images. First, it doubles the image acquisition time because the 790- and 1580-nm excitation wavelengths are applied one after another. Second, it is vulnerable to motion artifacts and mechanical drifts which can cause spatial mismatch when merging the 2PM and 3PM images. The acquisition of corresponding frames in 2PM and 3PM stacks are separated by the duration of imaging a whole stack which can last several minutes depending on the number of frames in a stack [16]. Meanwhile, motion artifacts are especially critical during in-vivo imaging.

To reduce motion artifacts and mechanical drifts in the sequential acquisition of 2PM (2PEF and SHG) and 3PM (THG) imaging, A. Filippi et al. used computer-controlled mechanical shutters to alternate wavelengths for 2PM and 3PM between consecutive frame acquisition [29]. The 2PEF and THG images were excited by wavelengths of 800 and 1200 nm, respectively, and the SHG images were excited by both wavelengths [29]. The frame-level interleaving mitigated inconsistency caused by sample motion or mechanical drift between 2PM and 3PM images. However, for each frame (depth) in a stack, the 2PM and 3PM images were still acquired sequentially, and the total imaging time would still be doubled compared with simultaneous acquisition.

Temporal multiplexing can further improve the simultaneous acquisition of 2PM and 3PM signals at the pixel level by interleaving 2PM and 3PM excitation pulses [30,31]. D. G. Ouzounov et al. interleaved 920-nm pulses from a Ti:sapphire laser modulated with an electro-optical modulator with 1300-nm pulses from a noncollinear optical parametric amplifier (400-kHz repetition rate) in 2.5-µs periods, respectively [30]. They performed simultaneous recording of neuron activity with 2PEF and 3PEF on GCaMP6-labeled neurons at the same site and time in order to validate functional imaging of 3PM by comparing it with 2PM [30]. The temporal multiplexing approach allowed the simultaneous acquisition of 2PM and 3PM signals within µs time delay. In order to increase the imaging volume in neuroscience study, S. Weisenburger et al. utilized spatiotemporal multiplexing of 2PM and 3PM to image multiple sub-volumes at the same time [31]. A custom laser system outputted 960 nm for 2PM and 1300 nm for 3PM excitation, at repetition rates of 4.7 and 1.0 MHz, respectively. They interleaved four 960-nm pulses and one 1300-nm pulse, each delayed by 8 ns through optical delay lines, forming five beamlets, each targeting a different sub-volume. Specifically, the 1300-nm beamlet was used to target a deeper axial location due to its higher penetration depth. Temporal multiplexing provides an alternative approach for signals that cannot be separated in the spectral domain. In both Refs. [30,31], 2PEF and 3PEF signals were acquired from GCaMP6-labeled neurons where the fluorescence emission signals could not be separately spectrally. However, temporal multiplexing generally requires high system complexity due to the requirements on modulators, optical delay lines, precise timing, and synchronization with the detection unit.

In this paper, we will address the challenges of applying two excitation wavelengths on tissue simultaneously for multimodal 2PM and 3PM imaging. We are aiming for a compact fiber-based system design and all label-free imaging. When two excitation wavelengths are shone on the sample simultaneously, the separation of 2PM and 3PM signals becomes difficult due to a partial wavelength overlapping between 2PEF and THG. Other challenges include the risk of thermal damage due to the additional excitation beam, and the difficulty to align two beams accurately for imaging the same field of view if two separate laser sources are utilized.

In addition to the challenges in the hardware design for data acquisition, the data processing of combined 2PM and 3PM imaging also has its unique challenges. While the different contrasts (e.g. 2PEF, SHG, and THG) visualize complementary tissue information, the signal level for each contrast can vary significantly, giving different detection channels highly differentiated signal-to-noise ratios (SNRs). This renders it difficult to visualize weak-signal channels, especially in merged multimodal images. For example, 3PM signals can be orders of magnitude weaker than 2PM signals [32]. Meanwhile, THG images from 3PM can reveal important information on tissue structures that lack 2PM contrasts, such as aqueous-lipidic interfaces in animal tissues and spores in plant tissues [16,20,33]. As another example, certain samples may only generate very weak (but still meaningful) signals on certain contrasts due to their structural and/or biochemical properties [16]. To reduce the amount of noise that would be introduced to the merged multimodal MPM images, images from the low-signal channels require denoising and SNR improvement before being merged with other channels.

Traditional denoising methods such as outlier removal and median filters are usually not effective in removing noise from those low-signal images while preserving tissue feature patterns [34]. Machine learning methods have also been utilized to denoise MPM images. Niu et al. demonstrated a neural-network-based method that effectively reduced noise from THG and 3PEF images of Hoechst blue labeled RAW 264.7 cells [35]. However, the feasibility of such methods typically relies on the availability of abundant training data [36]. In addition, the effectiveness of the trained neural network model may typically be limited to images with signal patterns similar to the training data. In recent years, self-supervised machine-learning methods have been developed to enable training using datasets without ground truth [37]. Noise2Void (N2 V) is a widely used self-supervised denoising algorithm that allows training using a single image without ground truth [38,39]. Nevertheless, machine-learning methods have the well-known black box problem, meaning that the trained model is typically not interpretable, which renders it difficult to explain how the image is processed [40].

In this paper, we will demonstrate a compact fiber-based multimodal MPM system with the capability of simultaneously acquiring 2PM (2PEF and SHG) and 3PM (THG) images. First, we will address the challenges related to data acquisition of simultaneous excitation using two wavelengths and the effective separation of 2PM and 3PM signals. Second, we will address the data processing challenges related to multimodal imaging. Specifically, a kernel-based nonlinear scaling (KNS) denoise and contrast enhancement method is introduced and demonstrated for improving images from low-signal channels and forming clean multimodal images. Comparisons are conducted between our denoising method and several other approaches including common traditional methods and N2 V. Finally, combined 2PM and 3PM imaging which simultaneously acquires 2PEF, SHG, and THG contrasts is demonstrated on biological samples, and the multimodal images are processed by our proposed denoising method to show the potential of this compact multimodal MPM system.

2. System setup

To perform simultaneous 2PM (2PEF and SHG) and 3PM (THG) imaging, a dual-wavelength excitation approach is used, where excitation wavelengths of ∼790 nm and ∼1580 nm are obtained from a frequency-doubled erbium-doped fiber laser with a pulse width of ∼80 fs and a repetition rate at 50 MHz. Our laser system enables fiber-delivery of femtosecond pulses, which brings about advantages in compactness and portability for a potentially handheld imaging head. With the frequency-doubling scheme, both 2PM and 3PM can be achieved with the single laser source and that further enhances the compactness.

Figure 1 shows the optical layout of the combined 2PM and 3PM imaging system. The 1580-nm laser is partially frequency-doubled to 790 nm by a periodically poled MgO:LiNbO₃ (PPLN) crystal (MSHG1550-0.5-0.3, Covesion). The maximum conversion efficiency for the PPLN is 35%. However, as THG can be 10 times weaker than 2PM signals [32], the PPLN is placed slightly out of focus to preserve a larger portion of 1580 power for THG excitation. Thus, the conversion efficiency is controlled at ∼18%. The laser power applied on samples is ∼9 and ∼40 mW at wavelengths of 790 and 1580 nm, respectively. The 790- and 1580-nm wavelengths are used to excite 2PM and 3PM signals, respectively. Compared with using two separate laser sources, the frequency-doubling solution avoids the difficulty of accurately aligning two beams since the fundamental and frequency-doubled beams are naturally overlapped. A miniaturized micro-electromechanical system (MEMS) mirror (13Z2.1-2400, Mirrorcle Technologies) performs XY scanning of the excitation beams. The pixel dwell time is set at 10 µs and that gives a frame rate of ∼0.4 FPS. The miniaturized objective is custom-designed and consists of an aspheric lens (352080, LightPath) and a plano-convex lens (45469, Edmund Optics) [16,41]. The objective has a focal length of f_OBJ = 3.6 mm and a numeric aperture (NA) of ∼0.5. The aspheric lens is mounted inside a shape memory alloy (SMA) actuator which enables axial scanning. Depth-resolved image stacks are acquired by scanning the aspheric lens in the depth direction (z). The lateral and depth resolutions for 2PM imaging are ∼0.7 and ∼9 µm, respectively; and those for 3PM imaging are ∼1 and ∼15 µm, respectively.

 

Fig. 1. Optical layout of the compact multimodal MPM imaging system with the illustration of the spectral relationship of the excitation and emission signals. L – lens; PPLN – periodically poled MgO:LiNbO₃; MEMS – micro-electromechanical system; OBJ – objective lens pack; MMF – multimode fiber; DBS – dichroic beam splitter; PMT – photomultiplier tube; λ_e – excitation wavelength.

Download Full Size | PDF

The emitted 2PM and 3PM signals from the sample are collected backward through the same objective and are separated from the excitation beams by a long-pass dichroic beam splitter DBS1 (FF665-Di02, Semrock) with an edge wavelength of 665 nm. Signals are then coupled into a multimode fiber and delivered to the detection unit. Before being detected by PMTs, the multimodal signals are separated by two additional dichroic beam splitters DBS2 and DBS3. DBS2 is a long-pass dichroic beam splitter (FF414-Di01, Semrock) with an edge wavelength of 414 nm. SHG signal is reflected by DBS2 and then detected by PMT1 (H6780-20, Hamamatsu). The longer-wavelength lights containing 2PEF and THG signals are further separated by a notch dichroic beam splitter DBS3 (ZT532dcrb, Chroma). DBS3 has a central reflection wavelength of 532 nm and a full-width-half-maximum (FWHM) notch bandwidth of 25 nm. THG signal is reflected by DBS3 and fed into PMT3 (H7422, Hamamatsu). As shown in the spectral relationship in Fig. 1, due to a partial spectral overlapping between 2PEF and THG, a small portion of the 2PEF signal within the 532/25 nm band is also fed into PMT3. Thus, PMT3 contains THG + 2PEF signals. Most part of the 2PEF signal passes through DBS3 and then reaches PMT2 (H9305-03, Hamamatsu). The spectral relationship of the two excitation wavelengths and the three nonlinear signals of 2PEF, SHG, and THG are illustrated in Fig. 1. A spectrograph (SP-2150i, Princeton Instruments) and an EMCCD camera (ProEM 512, Princeton Instruments) are used to measure and verify the spectrum of the emitted signals.

The multimodal MPM system in Fig. 1 contains multiple considerations for a compact design, including the fiber-delivered excitation laser and fiber-based signal collection, the miniaturized MEMS scanner and SMA actuator for XY and Z scanning, respectively, and the miniaturized objective with custom two-lens design. In the following sections, ex-vivo samples are imaged to demonstrate the proof of principle. Nevertheless, the compact design of this system will have the potential for in-vivo imaging.

3. Simultaneous acquisition of 2PM and 3PM signals

One challenge in the simultaneous dual-wavelength excitation comes from chromatic aberration of the objective lens between the two excitation wavelengths. For a compact design, the objective consists of a miniature aspheric lens and a plano-convex lens. While this lens design minimizes spherical aberration, there is still obvious chromatic aberration between the wavelengths of 790 and 1580 nm. Specifically, the 790- and 1580-nm wavelengths are focused to two different focal depths with a focal shift of ∼80 µm as shown in Fig. 2(a). The amount of focal shift is lens-specific, and it is calculated by Zemax simulation and verified by experimental calibration using a silicon photonic chip that contains structures with a thickness of hundreds of nanometers. Therefore, the acquired 2PM and 3PM image stacks would have a ∼80 µm depth offset, which can be compensated in post-processing by shifting the THG image stack downward by 80 µm. Figure 2(b) shows the verification of 2PM and 3PM signal spectra using the photonic chip sample. When the 790 nm is focused on the photonic chip, the SHG peak and the broadband 2PEF emission are detected in the spectrum; when the 1580 nm is focused on the photonic chip, the THG peak is observed in the spectrum.

 

Fig. 2. Acquisition and verification of the multimodal MPM with simultaneous 2PM and 3PM imaging. (a) Focal shift between the two excitation wavelengths; (b) 2PM and 3PM signal spectra; (c) process for crosstalk removal and adjusting the depths between 2PM and 3PM stacks. The merged image is the superposition of single-channel images with the blue, red, and green marks, respectively. THG image is denoised before merging. Red – 2PEF; blue – SHG; green – THG. Scale bar is 40 µm.

Download Full Size | PDF

Another challenge in simultaneous 2PM and 3PM imaging is the separation and simultaneous detection of multimodal signals. Based on the spectral relationship of the two excitation wavelengths and the 2PEF, SHG, and THG signals, dichroic beam splitters are carefully selected to spectrally separate the 2PEF, SHG, and THG signals. However, due to the partial spectral overlap between 2PEF and THG signals, PMT1, PMT2, and PMT3 detect SHG, 2PEF, and THG + 2PEF, respectively. Samples with both 2PEF and THG signals presenting in a considerable depth range may cause crosstalk from the 2PEF channel to the THG channel. Thin samples like the photonic chip are not significantly affected by this crosstalk issue. However, as shown in Fig. 2(c), the problem starts to appear when the sample thickness and imaging penetration reaches the focal shift of ∼80 µm. The relatively low depth resolution of the miniature objective may cause the excitation of out-of-focus signals which can add difficulty to the separation of the THG signal. To obtain a clean THG image stack, an image post-processing method is proposed to remove the leaked 2PEF signal from THG images by performing a subtraction of I_THG = I_PMT3 – r×I_PMT2 [42]. Here I_THG is the signal intensity of the THG image, and I_PMT3 and I_PMT2 stand for the signal intensity received by PMT3 and PMT2, respectively. In the equation of subtraction, the intensity of PMT2 is multiplied by a correction factor r because the two channels use different PMTs with different gains and the leaked portion of 2PEF signals can vary with sample.

The simultaneous acquisition of 2PM and 3PM signals doubles the imaging speed compared to sequential acquisition. Therefore, it can effectively reduce motion artifacts and mechanical drift between the 2PM and 3PM stacks. Figures 3(a) and 3(b) show the superposed SHG and THG images of a mouse femur bone acquired sequentially and simultaneously, respectively. The images are from stacks with 90 frames and acquired at a frame rate of 0.4 FPS. In the sequential acquisition, a relatively long time delay of ∼5 minutes exists between the 2PM and 3PM images from the same depth. The delay includes time for stack acquisition, manual wavelength switch, and saving the files. This long separation in acquisition time makes the imaging vulnerable to mechanical drift. In the current design of our compact MPM system, the sample is mounted vertically during imaging, and immersion water is added between the objective and the sample. Gravity pulling on the sample and immersion water can cause mechanical drift. Mechanical drift can also originate from tension and elastic deformation in tissues. In Fig. 3, the SHG contrast comes from collagen fibers while the locations of osteocytes lack SHG signal and appear dark. The oval-shaped structures showing THG signals correspond to osteocytes which should fit into the dark holes shown in the SHG image. However, as can be seen in Fig. 3(a), the mechanical drift of the sample during the time between imaging the corresponding frames in 2PM and 3PM stacks (∼5 min) caused motion artifacts on the superposed image. The location for the same osteocyte shifted ∼6 µm between the SHG and THG images. This motion artifact is minimized by the simultaneous acquisition of 2PM and 3PM stacks and the misalignment between SHG and THG images is significantly reduced as shown in Fig. 3(b).

 

Fig. 3. Motion artifacts comparison between (a) sequential and (b) simultaneous acquisition of 2PM and 3PM images of a mouse femur.

Download Full Size | PDF

While the chromatic aberration of the miniaturized objective causes the scanning focal shift of ∼80 µm and makes adjustment of stack images necessary when creating superposed images, the simultaneous scanning of two planes can potentially provide advantages in reducing thermal damage during multimodal imaging. To illustrate this, we have performed simulations of laser heating on tissue with and without the focal shift. The simulation is carried out using Ansys Workbench (Ansys Inc.) and the results are shown in Fig. 4. A volume of 200 × 100 × 100 µm³ of water is used to simulate the biological tissue and the environmental temperature is set to be 22 °C. Figure 4(a) shows the configuration of the simulated focal volume of single-plane imaging with a heat power of 4 mW applied. Figure 4(b) shows the configuration of dual-plane imaging where the two focal points are separated by 80 µm and a heat power of 2 mW is applied to each focal volume. The simulated temperature distribution across the section of the simulated volume under the settings of single- and dual-plane imaging is shown in Figs. 4(c) and 4(d), respectively. With the same total heat power, the highest temperature for single- and dual-plane imaging is ∼57 °C and ∼40 °C, respectively. Therefore, by imaging from two planes separated by 80 µm, the focal point temperature rise has been reduced by ∼50%, which can dramatically reduce tissue damage. This is critical in MPM where a large amount of laser power is needed in order to excite 2PM and 3PM signals.

 

Fig. 4. Configuration of thermal analysis and obtained temperature distribution of single vs. dual focal-plane imaging under the same total heat power. (a) and (b) are for single focal-plane imaging; (c) and (d) are for dual focal-plane imaging. The total heat power is 4 mW and the environment temperature is set to be 22 °C.

Download Full Size | PDF

4. Kernel-based nonlinear scaling denoising algorithm

To address the challenge of the large variation in the SNR among the multimodal images and to visualize weak-signal channels in merged multimodal images, we propose a kernel-based nonlinear scaling (KNS) denoising algorithm. Some preliminary results of the KNS algorithm were reported in our previous publication Ref. [16]. Here we will explain the KNS algorithm in detail, and provide improved results and more comprehensive comparisons with various denoising approaches.

The principle of this denoising algorithm is to adjust pixel brightness (intensity) based on a comparison between the global average brightness of the image and the local average brightness of the targeted pixel's surrounding kernel. Based on the main tissue features (e.g., epidermis cells, osteocytes, Volkmann’s canals) shown in the image, a square-shaped kernel size of n_K × n_K pixels is empirically chosen to be ∼50% of the size of the tissue feature. For each pixel (x, y) on an image with a dimension of X × Y pixels, a kernel centered on the pixel (x, y) with a dimension of n_K × n_K is then generated. The relationship between the scaled brightness I_s and the original brightness I_o of pixel (x, y) is defined as

Here I_GM is the global mean intensity of the entire image, and I_KM (x, y) is the local mean intensity within the kernel centered on the pixel (x, y). I_GM and I_KM (x, y) are respectively calculated as

This brightness scaling is performed on every pixel of the image. In case the targeted pixel is close to the edge of the image that causes the kernel (centered on that pixel) with the pre-selected dimension to exceed the image boundary, the actual effective kernel will be cropped by the image boundary to a reduced dimension.

Figure 5 shows representative images of low-SNR denoised by KNS and various commonly used denoising methods. The raw images include two THG images of mouse ear pinnae, an SHG image of the lamina of an Acer rubrum (red maple) leaf, and a THG image of a mouse femur bone. In our previous work [16], we demonstrated a comparison between an early-version KNS (the scaling ratio was based on a lookup table, and KNS has no iteration) and outlier removal using the mouse pinna and maple lamina images. Here in Fig. 5, additional denoising methods such as the Gaussian filter and bounded-variable least squares (BVLS) deconvolution algorithm have been added for comparison [43]. Furthermore, KNS is improved to be based on the empirical analytical scaling function Eq. (2) and be able to operate with multiple iterations. To show a more realistic SNR in generally dark images with few scattered signals, we also updated the definition of signal in SNR calculation from the global average brightness in Ref. [16] to the highest local average brightness. In this study, SNR is calculated by taking the ratio between the highest local average and the lowest local average, where the local averages are calculated for each image block with 20 × 20 pixels. The THG contrast of mouse pinna shows the cartilage layer, where the large (diameter of ∼30 µm) honeycomb structures could originate from the interfaces created by chondrocyte lacunae [44]. The SHG contrast in the maple lamina could be generated from starch in chloroplasts and cellulose in cell walls [45,46]. The THG image of mouse femur bone shows oval-shaped structures that indicate the location of osteocytes. The THG contrast could originate from the interfaces created by the osteocyte lacunae which are part of the lacuna-canalicular network. All the images contain 328 × 512 pixels. For KNS denoising, the kernel size is selected to be 20 × 20 pixels. KNS is operated on Matlab and the other methods are operated on ImageJ.

 

Fig. 5. Comparison of MPM images before and after denoising by outlier removal, Gaussian filter, BVLS deconvolution and KNS. SNR is calculated and shown in each image. σ – radius of the Gaussian filter; N – numbers of iteration. Cyan squares on KNS denoised images indicate the size of the kernel. Scale bar is 40 µm.

Download Full Size | PDF

As shown in Fig. 5, KNS effectively reduces noise on those low-signal images and can result in an SNR increase of more than 10 dB without significant loss of feature details. In comparison, commonly used traditional methods such as outlier removal, Gaussian filter, and BVLS fail to provide a noticeable SNR improvement to those low-signal MPM images. Thus, KNS outperforms those traditional methods. KNS is able to preserve the brightness of tissue features while reducing the background noise. Most tissue features such as the honeycomb structures in the mouse pinna are well-preserved, and the contrast of tissue features such as circular structures on the leaf lamina is increased. KNS can also be applied iteratively. In Fig. 5, KNS results with iteration numbers N = 1 and 2 are shown. Iteratively applying KNS may exponentially increase the SNR, as shown in Fig. 5. However, increasing the iteration number can result in a loss of image details.

Figure 6 compares the denoising results obtained by KNS (with iteration number N = 1 or 2) and N2 V. N2 V is a self-supervised deep-learning method that can be trained by a single image without ground truth [38,39]. N2 V training and prediction are operated using the CSBDeep-N2 V plugin of ImageJ [38,47]. The N2 V model is trained by one MPM image and then applied to three other low-signal MPM images as well as the training image itself. In N2 V training, the number of epochs, number of steps per epoch, batch size per step, patch shape, and neighborhood radius are set to be 300, 200, 56, 64, and 5, respectively. The KNS method is applied to the same images for comparison. The N2 V model in Fig. 6(a) is trained by a human bone 2PEF image with a relatively high signal level. The image shows osteocytes and their processes with fluorescence from endogenous fluorophores like NADH in the cytoplasm. The denoising result shows that except for the training image, KNS outperforms N2 V in all three low-signal testing images. On the human bone image (also the training image), N2 V shows better results as certain details including some canaliculi are dimmed or lost in the KNS results.

 

Fig. 6. Comparison between KNS and N2 V in processing low-signal noisy images. The N2 V models used in (a) and (b) are trained by (a1) a human bone 2PEF image with a relatively high SNR and (b1) a mouse pinna THG image with an ultra-low SNR, respectively. (a2) and (b2) show the training and validation loss through the training of these two models, respectively.

Download Full Size | PDF

In Fig. 6(b), the N2 V model is trained by a mouse pinna THG image which has an ultra-low signal level. The training does not show a typical successful process as it can be seen that while the training loss decreases through the training process, the validation loss increases. As a result, models from the last epoch and from the lowest validation loss show significantly different performances, and the last-epoch model fails to denoise images. The lowest-loss model is able to preserve some features but also introduces bright background artifacts. For all four images, KNS outperforms the N2 V model in removing the noise while preserving tissue features. On the computation side, both methods ran on a computer equipped with an Intel i9 11900F CPU and an Nvidia RTX 3080 GPU. The N2 V model costed ∼21.5 hours to train and took ∼1 s for each prediction; KNS took ∼34 ms per iteration for each image. While KNS cannot be generally expected to outperform deep-learning methods such as N2 V, in this comparison, KNS shows its specialty in processing images with ultra-low signal levels. Moreover, its simplicity enables fast computation and makes video-rate real-time processing possible.

To further evaluate the denoising performance, we also compare the peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM) using simulated data. In this simulation, a real image with high SNR is used as the ground truth (reference image); additive Gaussian noise is added to the ground truth to synthesize a noisy input image; KNS and N2 V methods are then applied to the noisy input image; finally, PSNR and SSIM are computed between the denoised image and the ground truth. Figure 7 shows the results of this comparison. A relatively clean THG image of Amanita muscaria (fly agaric mushroom) spores shown in Fig. 7(a) is used as the ground truth. In Fig. 7(b), Gaussian noise with a mean of zero and a standard deviation of 35 is added to the reference image to synthesize a noisy image which is then inputted to the denoising algorithms. Figure 7(c) shows the KNS denoising result with iteration number N = 2. Figures 7(d) and 7(e) are the denoising results using the N2 V model shown in Fig. 6(a) and another N2 V model trained by the input image Fig. 7(b), respectively. Each image is calculated for its PSNR and SSIM with reference to the ground truth. The PSNR and SSIM calculations are performed by the SNR plugin [48] and SSIM index plugin [49] of ImageJ, respectively. The input noisy image has a PSNR and SSIM of 17.51 dB and 0.251, respectively, compared to the ground truth. After denoising by KNS, the PSNR and SSIM improved to 23.54 dB and 0.459, respectively. In contrast, the N2 V processed two results reached only slightly improved PSNRs (18.06 and 18.11 dB) and even deteriorated SSIMs (0.220 and 0.190). In addition, the two N2 V models resulted in observably different results. Compared with the model trained by the bone sample with a good SNR (shown in Fig. 6(a)), the model trained by the input image (Fig. 7(b)) gives a much-blurred prediction (Fig. 7(e) compared with Fig. 7(d)). This shows that because of the black-box effect of machine-learning-based methods, they may generate inconsistent and unpredictable results due to different training and/or input data.

 

Fig. 7. Comparison of PSNR and SSIM obtained by KNS and N2 V. (a) Ground truth (reference image). (b) Noisy input image created by adding Gaussian noise with a standard deviation of 35 onto (a). (c) KNS (iteration N = 2) denoised result. (d) Denoised result by an N2 V model trained by a relatively high SNR image (shown in Fig. 6(a)). (e) Denoised result by an N2 V model trained by the noisy input image (b) with relatively low SNR.

Download Full Size | PDF

Denoising is particularly important in multimodal imaging. The multiple contrast images may have very different SNRs which can make the merged image very noisy or the weak channel(s) hard to distinguish. Figure 8 shows the contribution of the KNS denoising algorithm to the merging of multimodal images. A mouse pinna and a mouse femur are imaged by 2PEF, SHG, and THG. The THG signal is much weaker than the 2PEF and SHG signals. Without denoising, the noise level is close to the signal level on the THG images, and a large amount of noise is introduced to the merged images, affecting the readability of the other channels. In comparison, the KNS denoising algorithm effectively reduced the noise in the THG images. The merged images show significant improvement with tissue features well preserved and presented.

 

Fig. 8. Comparison of merging multimodal MPM images with and without KNS denoising on the low-signal THG channel. The gray and blue arrows indicate examples of noise and tissue features, respectively.

Download Full Size | PDF

5. Combined 2PM and 3PM imaging of biological samples

Multimodal MPM with simultaneously acquired 2PM and 3PM images are demonstrated in Fig. 9 on biological samples. The 2PM (2PEF and SHG) and 3PM (THG) signals are acquired simultaneously and label-free, enabling a comprehensive assessment of intricate microstructures while minimizing motion artifacts. Due to the focal shift, the THG image stack is shifted by 80 µm when merging the three channels.

 

Fig. 9. Simultaneous 2PM and 3PM imaging of biological samples. Red – 2PEF; blue – SHG; green – THG. Scale bar is 40 µm.

Download Full Size | PDF

Figure 9(a) shows multimodal MPM images of a fresh yew leaf. The SHG and THG channels show low-level signals and are denoised by KNS (iteration N = 1). Cell walls can be clearly seen as bright stripes in 2PEF on the leaf surface and dark stripes beneath the leaf surface. Sub-cellular structures can be identified as dark circular features on the images from a deeper depth (>20 µm). THG contrast shows chloroplasts inside leaf cells underneath the surface.

Figure 9(b) shows multimodal MPM images of a healthy human trabecular bone tissue fixed by formalin and mounted into a wax mold. The THG images are denoised by KNS (iteration N = 1). Osteocytes and canaliculi (lacuno-canicular networks) are visualized mainly by THG but also 2PEF contrasts from near the surface up to ∼100 µm deep into the tissue. Osteocytes play a vital role in bone homeostasis and mechano-sensation and are shown as oval-shaped structures in THG and 2PEF [50,51]. As THG requires interfaces [20,33], osteocytes may generate THG signals due to their complex morphology and refractive index variations created by the gap of 0.5–1.0 µm between the osteocyte cell body and the lacuna [52,53]. The 2PEF from osteocytes could be generated from endogenous fluorophores like NADH in the cytoplasm. The tiny fibrous structures in the image shown from depths 43, 55, and 74 µm are most likely canaliculi (canals where processes of osteocytes are embedded), which will maintain some extent of their micro-structures even with no metabolism. Given the sensitivity of THG to interfaces [20,33] and the role of canaliculi as essential communication channels for osteocytes, it is plausible that the observed THG signals originate from interfaces in the intricate networks of those cellular processes. The THG channel does not have sufficient SNR to image deeper beyond ∼100 µm while the 2PEF channel can image up to the depth of ∼170 µm. This is because three-photon nonlinear effects are generally weaker than two-photon nonlinear effects and the focal spot at a longer excitation wavelength (1580 nm, for 3PM) is not as tight as that at a shorter excitation wavelength (790 nm, for 2PM). The 2PEF contrast also visualizes endogenous fluorophores in other biomolecules and structures such as cartilage within the bone matrix. The SHG contrast delineates the collagen bundles in different depths, showing relatively regular distribution and organization, which confirms the healthy status of the sample.

6. Discussion and conclusions

In this study, a compact multimodal MPM system has been demonstrated to perform 2PM and 3PM imaging simultaneously. The dual-wavelength excitation enables the acquisition of multiple contrasts including 2PEF, SHG, and THG all label-free, which can obtain complementary information about tissues. To achieve simultaneous acquisition of 2PM and 3PM signals, several challenges are addressed. The 2PEF, SHG, and THG channels are separated spectrally by properly selected dichroic beam splitters. The challenge of partial wavelength overlapping between 2PEF and THG has been discussed and a signal-subtraction method is proposed to remove the leaked signal from images. The benefit of reducing motion artifacts and mechanical drift by simultaneous acquisition of multiple channels has been demonstrated. While the 2PM and 3PM image stacks are acquired simultaneously in our study, there is a focal shift of ∼80 µm between the 2PM and 3PM stacks. Thus, there still exists a time delay between the corresponding 2PM and 3PM frames from the same tissue depth. This time delay depends on the focal shift, depth scanning step size, and frame rate. To further reduce this time delay and its impact on motion artifacts during in-vivo imaging, several approaches can be carried out in future work. First, a faster scanning speed of up to ∼4 fps can be applied to reduce the time delay by 10 times. Second, the focal shift can be engineered and reduced by an achromatic design of the objective. Moreover, the focal shift can also be controlled if separate beam paths for the 790 and 1580 nm wavelengths are implemented. Nevertheless, this will increase the complexity and dimension of the system and thus is not ideal for the compact system. On the other side, the focal shift between 790 and 1580 nm can potentially be utilized to acquire images from two planes simultaneously. In addition, since the two beams are focused at two different depths, it can bring about a higher endurance of the total laser power compared to when the two wavelengths are focused on the same plane. When the focal shift is reduced, the laser-induced temperature rise in tissue can increase. Based on our thermal simulation, a focal shift of ∼20–40 µm could be a good balance.

To improve the image quality of low-signal channels and help form clean merged images, the KNS denoising method has been developed. With kernel size chosen by tissue feature dimensions, a kernel is generated centered on each pixel on an image, and the brightness of the central pixel is scaled by a nonlinear equation with variables based on the global and local (kernel) averages. The denoising method is generally universal to different tissue features and has been demonstrated on low-signal THG and SHG images acquired from animal and plant tissue samples. Comparisons between various denoising methods including outlier removal, Gaussian filter, deconvolution, and N2 V have been conducted to verify the denoising performance. KNS outperforms the above-mentioned denoising methods on low-signal MPM images. The reason for KNS to work is that the method judges a pixel’s possibility to be noise by its surrounding brightness rather than the pixel itself. This idea provides KNS the capability to distinguish signal and noise from images with signal intensities that are so low that even the part showing signals consists of highly sparkled pixels. As a result, KNS is particularly effective in processing ultra-low signal MPM images. However, for input images already with a relatively high SNR, KNS may cause artifacts that dim some signals. As a scaling method, KNS can alter the relative brightness between different tissue features shown on an image and may not be suitable for quantitative analyses. This effect can be more obvious on bright images with most parts showing signals rather than background noise. With the increased number of iterations, KNS may also cause a loss of image details, saturations, and gridded artifacts. Overall, in the merging of MPM images consisting of multiple channels with highly inconsistent SNRs, KNS improves the image quality by significantly reducing the excessive amount of noise introduced to the superposed image from low-signal noisy channels. KNS provides an effective and fast denoising approach for low-signal images whenever ground truth images are not available.

To demonstrate its capability, depth-resolved multimodal MPM images from biological samples are demonstrated with 2PEF, SHG, and THG signals acquired simultaneously by our compact system. The multimodal imaging approach offers a comprehensive view of microstructures in tissues, facilitating a deeper understanding of tissue composition and potential cellular interactions. Our technique of simultaneous acquisition of 2PM and 3PM signals shows great potential in improving MPM imaging by accelerating the imaging process and reducing possible image artifacts caused by tissue motion and mechanical drift during sequential acquisition. The proposed KNS denoising method further improves the quality of multimodal images and has the potential to be integrated into real-time clinical imaging systems due to its simplicity in computation.