2.1.

Sample Acquisition

All activities pertaining to cervical tissue biopsy handling were done in accordance with approved Tufts Health Sciences IRB protocol #10283. Patients over the age of 18 with a recent low-grade squamous intraepithelial lesion (LSIL) or high-grade squamous intraepithelial lesion (HSIL) pap smear diagnosis undergoing a colposcopy or loop electrosurgical excision procedure were recruited to the study. Informed consent was acquired from all study subjects before participation. During the routine procedure, a second biopsy from a colposcopically abnormal region of the cervix was taken and placed in a custom-built tissue carrier containing keratinocyte serum-free media (Lonza). Biopsies were transported via personal vehicle to the Tufts Advanced Microscopy Imaging Center for imaging. All imaging was conducted within 4 h post-biopsy. Immediately after imaging, biopsies were fixed in 10% neutral buffered formalin. Biopsies were returned within 5 business days to the Tufts Medical Center Department of Pathology for standard histopathological diagnosis.

Patients over the age of 18 undergoing hysterectomies for benign gynecological disease were also recruited to the study as healthy controls. The only difference between healthy and precancerous biopsy acquisition was in the actual biopsy excision. Healthy biopsies were sampled from the resected cervix by a pathologist after macroscopic inspection to rule out abnormalities.

2.2.

Deep Learning Dataset Details

A total of 151 regions of interest (ROIs) (image stacks) were collected from 54 patients. The training and validation sets were comprised of 100 ROIs featuring 5 to 50 optical sections (OSs) per ROI. 75% of the ROIs were randomly selected for training, and the remaining 25% were set aside as the validation set (1657 training OSs and 554 validation OSs). To prevent data leakage, training and validation OSs were separated on an ROI basis. The test set featured 51 ROIs (with 10 to 50 OSs per ROI) and was excluded from all training (1018 OSs). For $k$-fold validation, the 100 training and validation ROIs were shuffled and split again using the same 75:25 ratio for up to five times to ensure robustness of denoising on a constant test set (see Fig. S5 in the Supplementary Material). The dataset features images from tissues with three diagnoses: benign, LSIL, and HSIL. The test set was composed of 25 benign ROIs (49.02%), 14 LSIL ROIs (27.45%), and 12 HSIL ROIs (23.53%). The training and validation sets were composed of 55 benign ROIs (54.45%), 25 LSIL ROIs (24.75%), and 21 HSIL ROIs (20.79%). Based on training/validation splitting seed, these values could range from 52% to 57.3% benign, 25.3 to 26.7% LSIL, and 18.7% to 22.7% HSIL in the training set and 48% to 64% benign, 20% to 24% LSIL, and 16% to 28% HSIL in the validation set. An alternative training scheme was initially attempted. In this scheme, only benign ROIs were used in training with 112 ROIs of mixed diagnosis being used in the test set and 39 benign ROIs being used for training. The training set was later modified as it became evident that training on a mixture of diagnoses resulted in superior restoration of downstream metrics (Fig. 5).

2.3.

Optical Instrumentation and Image Acquisition

Images were collected using a commercially available Leica SP8 inverted microscope system equipped with an Insight fs laser. Tissue biopsies were placed epithelial side down onto a glass bottom dish and light was delivered using an epi-illumination scheme. 10 to 60 mW of power were delivered at the tissue surface with power being varied linearly through the depth of the tissue during acquisition. The rate of power increase was determined based on achieving minimal pixel saturation in the most superficial and most basal OSs. The 60 mW threshold was motivated by the finding that this power equates to a 0.6 minimal erythma dose. The threshold of sunburn development is a 1.0 minimal erythma dose.⁴² The maintenance of safe excitation power levels supports the translational relevance of the present dataset. Images were acquired using bidirectional scanning with a 600 Hz line scan frequency and a 400 ns pixel dwell time. Tissue biopsies were excited with 755 and 860 nm light. Two hybrid photodetectors were set up to collect the two-photon autofluorescence signal from NAD(P)H ($460\pm 25\mathrm{nm}$) and FAD ($525\pm 25\mathrm{nm}$). Hybrid photodetectors were used in place of photomultiplier tubes due to their increased sensitivity to low intensity fluorescence signal compared with photomultiplier tubes.^43,44 Two photomultiplier tubes were set up to collect the second harmonic generation signal from collagen fibers ($430\pm 12\mathrm{nm}$) and the red autofluorescence signal ($624\pm 20\mathrm{nm}$). Light was delivered and collected using a 40X/1.1 NA water-immersion objective lens ($290\times 290\text{-}\mu {\mathrm{m}}^2$ field-of-view). Images were collected through the full thickness of the epithelium using a depth-sampling rate of $4\mu \mathrm{m}$. Six individual frames were collected at each depth. On average, 3 to 5 ROIs were sampled from each biopsy.

2.4.

Morphological and Functional Metrics

Images were calibrated and processed as described in detail previously to extract images that represented NAD(P)H and FAD TPEF intensity fluctuations.^23–25,45 At each optical depth, NAD(P)H and FAD images were used to define a corresponding redox ratio for each pixel of the field, given as

From the RR distributions for each OS, we calculated the mean RR and the interquartile range (IQR) as metrics of the overall oxidation-reduction tissue state and the corresponding heterogeneity, respectively. The mean and sample variance (variability) of the mean OS RR and the OS RR IQR for all images in an epithelial stack were calculated to assess the depth-dependence of these metrics.

NAD(P)H images were analyzed as described previously^17,18,21,22 using a Fourier based approach to extract a value for the parameter $\beta$ as a metric of the level of mitochondrial fragmentation and networking, which also depends highly on the metabolic activity of the tissue. Briefly, an inverse power law was fit to the power spectral density (PSD) of the two-dimensional Fourier transform of the cytoplasmic NAD(P)H intensity fluctuation images, given as

2.5.

Deep Learning Model Description

The basic structure of the CARE network has been described extensively [Fig. S3(a) in the Supplementary Material].³⁴ The network was implemented through Keras and TensorFlow.^46,47 A copy of the CSBDeep repository (available in a Github repository: https://github.com/CSBDeep/CSBDeep) was locally imported into an anaconda environment.⁴⁸ The network was configured to take a $256\times 256\times 1$ input image and generate a $256\times 256\times 1$ denoised image. A 40-gigabyte Nvidia Tesla A100 GPU card was used for all training and evaluation. Typically, a $1\times 512\times 512\times z$-depth image stack was split into $4\times 256\times 256\times z$-depth image patches before training using a $2\times 2$ grid as previously described.⁴⁹ A starting learning rate of $1\times {10}^{-5}$ was used with an Adam optimizer.⁵⁰ Training was allowed to continue for 300 epochs with a scheduler reducing the learning rate when the network performance stagnated for more than 20 epochs. Early stopping was not implemented to allow the model to improve with lower learning rates. Model improvement was generally observed to stagnate after 75 to 100 epochs with the best model weights being saved. The loss functions were varied to find the optimal function to improve the downstream analysis performance. The loss functions used include SSIM loss, R2 loss, FFL, MAE (L1) loss, MSE (L2) loss, Redox ratio loss, and combined losses, such as a combined SSIM + L2, SSIM + FFL, SSIM + R2 loss.⁵¹ Combined loss functions were weighted using an $\alpha$ term to control the contribution of SSIM loss and L2, FFL, R2 loss. An $\alpha =0.84$ was used in this study, in line with other studies, to balance the contribution of both loss functions.⁵² Six down-sampling and up-sampling layers were generated with the first layer expanding the single-channel images to 32 channels. Residual connections were used to preserve encoded information from each down sampled layer and pass it forward to the decoder layers (see Fig. S3 in the Supplementary Material).

For the wavelet U-Net (WU-net) architecture, four CARE networks, one per sub-band, were built as described above. A discrete wavelet transform (DWT) was used to decompose a $1\times 256\times 256$ optical section patch (OSP) into $1\times 128\times 128\times 4$ frequency band images. The four frequency bands would then be individually input to each CARE network for denoising. After denoising, an inverse DWT (iDWT) was used to reconstruct the $1\times 256\times 256$ OSP (for greater detail see Fig. S4 in the Supplementary Material).

The training time typically varied from 1 to 2 h, with an evaluation time of $\sim 24\mathrm{s}$ per image stack. For all trained CARE networks, 3D SSIM, PSNR, mean $\beta$, $\beta$ variability, mean RR, RR variability, RR IQR, and RR IQR variability were analyzed. All final metrics were assessed using a single frame input, denoised, and ground truth (six frame averaged) images with built-in and custom MATLAB (MathWorks; Natick, Massachusetts) functions.

2.6.

Statistics

For Figs. 3Fig. 4–5(c) and Table 7, the Fisher $r$-to-$z$ transformation was used to convert Pearson’s correlation coefficients ($r$) to $z_r$ values.⁵³ This transform was calculated as

The ${\mathrm{z}}_r$ value, unlike $r$, belongs to a normal distribution, allowing for the calculation of a $Z$-statistic to determine confidence intervals. The test $Z$-statistic for comparison of ${\mathrm{z}}_r$ values to determine significance was calculated as

3.1.

Identification of the Optimal Deep-Learning Model Architecture for Denoising Label-Free, Optical TPEF Images to Enable Recovery of Metabolic Function Metrics

Human cervical tissue biopsies were collected from 54 patients and imaged immediately upon excision, as described in Sec. 2.3 (Fig. 1). Several ROIs were imaged from each biopsy. Multiple OSs were imaged from each ROI at distinct depths. At each OS, we acquired TPEF images at a combination of two excitation wavelengths (755 and 860 nm) and three or four emission bands. Images collected at 755 nm excitation and 435 to 485 nm emission were attributed primarily to NAD(P)H, whereas images at 860 nm excitation and 500 to 550 nm emission were considered to contain signal primarily from FAD and FAD bound to lipoamide dehydrogenase. Six frames were acquired at each wavelength setting. To reduce the contribution of noise, these six frames were averaged together. The decision to use six frames was made based on data collected from freshly excised rat oral squamous epithelial tissue. Analysis of this dataset highlighted that six-frame images conveyed quite accurately the depth dependence of the metabolic function metrics when compared with the analysis of 32 frame-averages (see Fig. S6 in the Supplementary Material). Metrics extracted from these averaged images were previously observed to enable highly sensitive and specific detection of cervical pre-cancer.²⁵ The averaged image was therefore considered the ground truth used for training and testing the denoising success of single frames. Single frames, the corresponding denoised images, and ground truth images were analyzed using established procedures to extract the RR and mitochondrial clustering ($\beta$) (Fig. 1). All models (Fig. 1) were trained and evaluated with identically generated image stacks. Various combinations of model architectures, loss functions, data transformations, and training data combinations, as outlined in Table 1, were evaluated on 3229 total OSs (1657 training OSs + 554 validation OSs + 1018 testing OSs) representing healthy/benign cervical tissues as well as precancerous (low-grade and high grade) squamous intraepithelial lesions (LSIL and HSIL, respectively).

Fig. 1

Summary of deep learning pipeline. Human cervical tissue biopsies are collected and subsequently imaged within 4 h post-excision. Collected biopsies are plated on glass bottom dishes and imaged using a Leica SP8 commercial microscope. At a minimum, three ROIs are imaged per sample. At each ROI, multiple OSs are imaged at distinct depths through the epithelium. Depth-resolved, two-photon OSs are collected using two excitation wavelengths and several bandpass-filtered detectors. Six images are captured for each excitation/emission wavelength and every OS at a given depth, $z$. These six images are averaged together to generate the ground truth image set. A random image from the six per depth $z$ is selected as the input (RAW) image. The paired image stacks are provided to the neural network for training and denoising. Four-leading denoising networks are used in this study to denoise input images: a previously described CARE model, an RCAN model, an SRGAN model, and a WU-Net.^34–36,38 Denoised images and input images are compared against 6× averaged images to determine 3D PSNR and SSIM along with metabolic metrics. $\text{Scale bar}=25\mu \mathrm{m}$.

Table 1

Summary of all parameters explored during training and optimization of the final model (highlighted in bold). Results shown below are focused on the optimized model, but all combinations were trained and evaluated.

PSNR and SSIM improvements are standard metrics of image visual quality and have been used in other studies focused on denoising biomedical images as indicators of model success.^34,35,37,38 We aim to assess whether images restored by models that yield optimized PSNR and SSIM values result in accurate recovery of metabolic metrics (Fig. 2). For evaluation of the model architecture, loss function, and signal type, only results from models trained on NAD(P)H data from tissues of known benign status were included.

Fig. 2

$\beta$ and RR metrics are extracted from each OS (see Sec. 2.4 for greater detail). Depth-dependent trends across the multiple cell layers of the cervical squamous epithelium are assessed for input images (RAW 1X), denoised images (restored) and six-frame averaged, ground truth images (GT 6X). Measurements of mean values and corresponding variability across all depths are shown as a dashed line and shaded region in the mitochondrial clustering, $\beta$, and RR (FAD/(NAD(P)H+FAD)) panels for the GT 6X image. The distribution of RR values for each OS is used to extract the IQR, representing the range of values within the 25% and 75% of the RR distributions and providing an assessment of intra-field RR heterogeneity. IQR variability is a metric of inter-field (depth-dependent) RR heterogeneity.

Leading denoising model architectures were selected for evaluation based on a comprehensive literature search. CARE, RCAN, and SRGAN (Fig. 1) models were trained as described in Sec. 2.5 and Supplementary Methods: Deep Learning Performance Benchmark in the Supplementary Material. A representative OS from an LSIL biopsy is shown in Fig. 3(a). The results shown were generated by models trained using an SSIM + MSE (or L2) loss function. A summary of all parameters used to generate the figures and tables is listed in Table 2.

Fig. 3

(a) A $290\times 290\mu {\mathrm{m}}^2$ field of view from an LSIL cervical tissue biopsy. NAD(P)H and FAD images for the same region are shown along with the corresponding denoised image from each of the three trained models (CARE, RCAN, and SRGAN). $\text{Scale bar}=50\mu \mathrm{m}$. (b) A $44.2\times 44.2\mu {\mathrm{m}}^2$ field of view (white square in a) of three cells. NAD(P)H and FAD images are shown for all models and the input and ground truth images. $\text{Scale bar}=10\mu \mathrm{m}$. (c) Bar plots of the coefficient of determination of all downstream metrics for images denoised by all models and RAW 1X versus the GT 6X image. Fisher $r$ to $z$ transformation was used to measure significance. *$p<0.05$ and **$p<0.01$.

Table 2

Summary of parameters used to generate Figs. 3–6 and Tables 3–6. Parameters are bolded when all combinations from Table 1 are used.

Prior to denoising, standard image quality metrics were calculated for input (RAW 1X) images by comparing the RAW 1X images to ground truth (GT 6X) images. PSNR and SSIM values were calculated using the GT 6X image as a reference and RAW 1X or denoised images as the distorted image.⁵⁵ Across all images, FAD image PSNR was greater than NAD(P)H image PSNR (Table 3), even though FAD images featured lower cytoplasmic signal compared with NAD(P)H images [Fig. 3(a)]. During PSNR calculation, the reduced signal intensity led to smaller differences between RAW 1X and GT 6X images and yielded a greater observed PSNR value. This observation was also consistent with the results from other studies.³⁷ SSIM values were consistent between NAD(P)H and FAD images (Table 3). Corruption of the GT 6X images for both channels by noise was expected to have similar effects on structural similarity and calculated SSIM values.

Table 3

Summary of standard metrics of image quality for RAW 1X images and denoised images generated from various model architectures. Values are reported for mean performance (± standard deviation) across all test set ROIs.

We used 777 and 109 RAW 1X NAD(P)H OSs for training and validation of the models, respectively. Each $512\times 512$ OS was patched into four-$256\times 256$ image patches (OSP) prior to training and validation (3108 and 436 OSPs, respectively). All three models were trained before being evaluated on an independent set of 2343 OSs (9372 OSPs). Metrics of image quality and metabolic function were calculated as described in Supplementary Methods: Deep Learning Metrics in the Supplementary Material and Sec. 2.4.

CARE-generated image stacks demonstrated higher PSNR for FAD images and higher SSIM for both NAD(P)H and FAD images compared with restored-image stacks generated by RCAN and SRGAN. Across all test set images, standard metrics of image quality (Table 3) and visual inspection [Fig. 3(b)] suggested the RCAN- and CARE-denoised images had similar image quality. Across the entire test set, we observed that SRGAN failed to restore cellular features in the RAW 1X that are found within the GT 6X images [Fig. 3(b)] and underperformed even relative to RAW 1X images in standard image quality metrics (Table 3). Perceptual loss was believed to impact content restoration in the SRGAN architecture.³⁵ Inputs for perceptual loss calculations have been shown to impact significantly the SRGAN performance and were likely the cause of SRGAN’s poor recovery of the image quality.³⁵

To assess the restoration of metabolic activity, depth-dependent optical RR and mitochondrial clustering ($\beta$) values were calculated for the restored images, input (RAW 1X) images, and ground truth (GT 6X) images (Fig. 2). Perturbations in metric values are reflective of changes in the metabolic functional state of the different cell layers in tissue.^{21,22,24,25,39,40} Pearson correlation coefficient values were calculated between the metabolic function metrics from the GT 6X and either the RAW 1X or restored images. Statistical significance was derived from Fisher-$r$-to-$z$ transformation for all metrics of interest. Interestingly, analysis of the RAW 1X images led to very high correlations with metrics of RR intra- and inter-field variability compared with GT 6X images. We hypothesized that similar sources of noise in both FAD and NAD(P)H images led to this outcome because RR metrics were calculated using a ratio of FAD and NAD(P)H intensity measurements. It was for this reason that, in this initial comparison, we trained models on NAD(P)H images and applied the weights to NAD(P)H and FAD images. RCAN-generated images demonstrated a statistically significant recovery of $\beta$ variability (${\sigma }^2$ ($\beta$)) [Fig. 3(c)]. However, recovery of mean RR variability by this model was poor [Fig. 3(c)]. CARE-denoised images, overall, featured higher (albeit not statistically significant) correlations with RR metrics compared with all other models [Fig. 3(c)]. Thus, the U-Net architecture of CARE was utilized for all further optimization steps.

3.2.

Multiscale Image Transformation Enhances Quantification of Mitochondrial Clustering

Although denoising improved the restoration of ${\sigma }^2$ ($\beta$), the mean $\beta$ ($\overline{\beta }$) values of the denoised images were not well correlated with the values from the GT 6X images. We considered DWT to enhance high spatial frequency restoration necessary for $\beta$ metric calculations. A single level DWT transformed each image into four sub-band images: a coarser scale approximation and three detail images, one horizontal, one vertical, and one diagonal.⁵⁶ To generate the three subband images, a basis function, called a mother wavelet, was convolved along both dimensions of the original image, and an associated scaling function was used to generate the coarser approximation. During standard wavelet-based denoising, thresholds are used to remove noise from wavelet-transformed detail images, before implementing an inverse-transform to recover the restored image.^57,58 The DWT has been shown to be advantageous compared with traditional low-pass filtering as the pixel-by-pixel convolution with the mother wavelet preserves correlations of high frequency features. In this study, we used deep learning models trained on each of the transformed images to adaptively learn the best threshold for denoising of low frequencies (approximation) and high frequencies (details) rather than relying on arbitrary thresholding for denoising (see Discussion S2 in the Supplementary Material).⁵⁹ As with any DWT-denoising model, the selection of the correct mother wavelet played a significant role in model performance. For all models, mother wavelets from the biorthogonal, coiflets, and Daubechies families were evaluated. These mother wavelet families were selected due to their frequent use in denoising tasks.⁶⁰ Multiple models were trained and evaluated, with biorthogonal 1.1 yielding the highest recovery of metabolic metrics (data not shown). As such, biorthogonal 1.1 was used for all subsequent model optimization.

The application of DWT before training four CARE models and iDWT after evaluation yielded images with improved FAD and NAD(P)H PSNR with slight decreases in SSIM [Fig. 4(a)]. Across the entire test set, NAD(P)H PSNR improved using WU-Net, whereas FAD PSNR and SSIM both decreased compared with CARE (Table 4). All loss functions were evaluated for WU-Net with SSIM + R2 loss (results shown in Fig. 4) and SSIM + FFL loss (see Table S4 in the Supplementary Material) yielding the best overall performance. WU-Net denoised NAD(P)H images extracted similar cellular structures as the CARE derived images but featured a lower background signal and small fluctuations in cytoplasmic signal, leading to the observed higher PSNR values [Fig. 4(b)]. WU-Net led to statistically significant improvements in the correlation of extracted ${\sigma }^2$ ($\beta$) with GT 6X images relative to the analysis of the RAW 1X images. Extracted $\overline{\beta }$ values were also better correlated to GT 6X images, albeit improvements were not significant.

Fig. 4

(a) A $290\times 290\mu {\mathrm{m}}^2$ field of view from an LSIL cervical tissue biopsy. NAD(P)H and FAD images for the same region are shown along with the corresponding denoised image from each signal pre-processing method utilized (single frame image and wavelet transformation). $\text{Scale bar}=50\mu \mathrm{m}$. (b) A $44.2\times 44.2\mu {\mathrm{m}}^2$ field of view (white square in a) of three cells. NAD(P)H and FAD images are shown for all models based on the corresponding signal pre-processing method used during training and the input and ground truth images. $\text{Scale bar}=10\mu \mathrm{m}$. (c) Bar plots of the coefficient of determination of all downstream metrics for images denoised by all models trained based on the corresponding signal pre-processing method used during training and RAW 1X versus the GT 6X image. Fisher $r$ to $z$ transformation was used to measure significance. *$p<0.05$ and **$p<0.01$.

Table 4

Summary of standard metrics of image quality for RAW 1X images and denoised images generated after signal pre-processing. Values are reported for mean performance (± standard deviation) across all test set ROIs.

Fig. 5

(a) A $290\times 290\mu {\mathrm{m}}^2$ field of view from an LSIL cervical tissue biopsy. NAD(P)H and FAD images for the same region are shown along with the corresponding denoised image from each data type used as training data for the WU-Net model (NAD(P)H healthy only-, NAD(P)H mixed diagnosis-, and FAD mixed diagnosis-wavelet transformed images). All models were equally constructed with only the data type and diagnosis type varied. $\text{Scale bar}=50\mu \mathrm{m}$. (b) A $44.2\times 44.2\mu {\mathrm{m}}^2$ field of view (white square in a) of three cells. NAD(P)H and FAD images are shown for all data types used during training and the input and ground truth images. $\text{Scale bar}=10\mu \mathrm{m}$. (c) Bar plots of the coefficient of determination of all downstream metrics for images denoised by models trained using varying data types and diagnosis types and RAW 1X versus the GT 6X image. Fisher $r$ to $z$ transformation was used to measure significance. *$p<0.05$ and **$p<0.01$.

Comparing WU-Net with an identical CARE model, we observed that WU-net achieved an improved performance on $\beta$ metrics while maintaining the recovery of RR metrics [Fig. 4(c)]. The overall improved $\beta$ restoration suggested that WU-Net was better able to capture true signals from noise in the high spatial frequencies found in NAD(P)H images. WU-Net further preserved the relationship between NAD(P)H and FAD channel images, enabling high correlations for RR metrics. Due to the observed performance of WU-Net for $\beta$ metric recovery, we explored further optimization of WU-Net, which could be achieved by varying the training datasets.

3.3.

Selection of Training Data

Initial model development focused on a limited training set of cervical tissues of known benign status (healthy). Benign tissue samples comprise cell layers with consistent changes in differentiation as a function of depth among image stacks. Training on such images was expected to enable the model to learn characteristics of noise without having to account for feature heterogeneity found in pre-cancerous cervical tissue samples. We further sought to assess whether training on a data set that was expanded to include image stacks from tissues with both benign and pre-cancerous lesions (Mix) impacted performance. In this new training set, 1657 and 554 RAW 1X NAD(P)H OSs (6628 and 2216 OSPs) were used for training and validation of the models, respectively. An independent test set of 1018 OSs (4072 OSPs) was used to evaluate the model performance after training.

An additional consideration that we explored was the impact of the source of image contrast, i.e., NAD(P)H or FAD, used for training. NAD(P)H images featured greater structural information compared with FAD images, and they were utilized in our analysis for the extraction of mitochondrial clustering-focused metabolic function metrics (Figs. 3Fig. 4–5). Thus, training was focused on NAD(P)H images, and optimized model weights from NAD(P)H image training were used to denoise FAD images for extraction of RR-based metrics. However, because similar noise characteristics were assumed to be present in both RAW 1X NAD(P)H and FAD images, we sought to confirm that training on NAD(P)H images was optimal. Thus, we used FAD images to train WU-Net models using the same hyperparameters and settings as the ones used when NAD(P)H images were used. Post-training, NAD(P)H images were denoised using the weights of the FAD image trained model to extract RR and mitochondrial clustering-based metrics.

The use of training sets with mixed diagnosis images resulted in minimal differences in the denoised images when compared with training just on healthy sample images [Fig. 5(a)]. PSNR and SSIM values for images were observed to be nearly identical because of these insignificant differences (Table 5). Both models led to denoised images with consistent cell boundary and intracellular structures given the same RAW 1X images [Fig. 5(b)] and had similar levels of restoration of downstream metrics, with the mixed diagnosis dataset leading to slightly improved correlations in most cases [Fig. 5(c)]. The increase in correlation could be attributed to the large training set available for a mixture of diagnoses compared with only training on healthy data.

Table 5

Summary of standard metrics of image quality for RAW 1X images and denoised images generated after training models on various data types. Values are reported for mean performance (± standard deviation) across all test set ROIs.

An identical model was trained using the FAD image data from the mixed diagnosis dataset. The denoised images from the FAD-trained model looked like those from the corresponding NAD(P)H-trained model [Figs. 5(a) and 5(b)]; however, the standard metrics of image quality were slightly lower. Images denoised by the FAD-trained model demonstrated a higher background signal compared with images denoised by NAD(P)H-trained models [Fig. 5(b)]. However, despite FAD images lacking much of the structural and morphological information of their NAD(P)H counterparts, their use in training led to further improvements in $\beta$ metric recovery and mean RR restoration from the RAW 1X images [Fig. 5(c)]. We hypothesize that high frequency information in the FAD images originated primarily from noise in comparison with NAD(P)H images. As a result of the high frequency information containing primarily noise, the model improved in its learning of noise characteristics in the images, enabling improved denoising and recovery of metrics of metabolic activity [Fig. 5(c)].

3.4.

Summary of Final Model Performance

Across all models, image quality improved after denoising based on PSNR and SSIM (Table 6). Based on standard image quality metrics of all models discussed in this study, it could be assumed that models trained using NAD(P)H images and the CARE architecture with standard loss functions of MAE and MSE would perform best at the restoration of downstream metrics [Fig. 6(a)]. CARE models trained with MAE and MSE loss functions both demonstrated statistically significant improvement in denoised FAD and NAD(P)H image PSNR and SSIM ($p<0.001$). Comparatively, wavelet-transformed-FAD images denoised using WU-Net with SSIM + R2 loss had poorer standard metric performance (Table 6). Images restored with this model did not achieve statistically significant improvement of FAD image PSNR and SSIM [Fig. 6(a)]. As PSNR and SSIM are commonly used as indicators of model performance, it was expected that improvements in these metrics would correspond to better recovery of downstream metabolic metrics. However, the WU-net model trained on mixed diagnosis, FAD images led to denoised images with extracted metabolic metrics that were consistently correlated with the metrics extracted from GT 6X images [Fig. 6(b)]. The final correlations of the models shown in Fig. 6 are reported in Table 7.

Table 6

Summary of standard metrics of image quality (PSNR and SSIM) for RAW 1X images, standard implementation of CARE, and the best performing model from this study. Values are reported for mean performance (± standard deviation).

Fig. 6

(a) Box and whisker plots of PSNR and SSIM of 40 test set ROIs. Denoised images demonstrated an improvement in standard metrics of image quality. (b) Bar plots of the coefficient of determination of all downstream metrics for images denoised by models trained using various data types, loss functions, and diagnosis types versus the GT 6X image. A one-way ANOVA with Tukey Kramer post-hoc test was used to measure significance of PSNR and SSIM. Fisher $r$ to $z$ transformation was used to measure significance of improvement in metabolic metric correlations *$p<0.05$, **$p<0.01$, and ***$p<0.001$.

Table 7

Correlation values of models in Fig. 6. Fisher r to z transformation was used to measure significance.