Framework for denoising Monte Carlo photon transport simulations using deep learning

Abstract. Significance The Monte Carlo (MC) method is widely used as the gold-standard for modeling light propagation inside turbid media, such as human tissues, but combating its inherent stochastic noise requires one to simulate a large number photons, resulting in high computational burdens. Aim We aim to develop an effective image denoising technique using deep learning (DL) to dramatically improve the low-photon MC simulation result quality, equivalently bringing further acceleration to the MC method. Approach We developed a cascade-network combining DnCNN with UNet, while extending a range of established image denoising neural-network architectures, including DnCNN, UNet, DRUNet, and deep residual-learning for denoising MC renderings (ResMCNet), in handling three-dimensional MC data and compared their performances against model-based denoising algorithms. We also developed a simple yet effective approach to creating synthetic datasets that can be used to train DL-based MC denoisers. Results Overall, DL-based image denoising algorithms exhibit significantly higher image quality improvements over traditional model-based denoising algorithms. Among the tested DL denoisers, our cascade network yields a 14 to 19 dB improvement in signal-to-noise ratio, which is equivalent to simulating 25× to 78× more photons. Other DL-based methods yielded similar results, with our method performing noticeably better with low-photon inputs and ResMCNet along with DRUNet performing better with high-photon inputs. Our cascade network achieved the highest quality when denoising complex domains, including brain and mouse atlases. Conclusions Incorporating state-of-the-art DL denoising techniques can equivalently reduce the computation time of MC simulations by one to two orders of magnitude. Our open-source MC denoising codes and data can be freely accessed at http://mcx.space/.

assessment in the bedside or natural environments. 1 However, the main challenge of using lowenergy NIR photons is the high degree of complex interactions with human tissues due to the presence of high scattering, which is much greater than that of x-rays. As a result, the success of many emerging NIR-based imaging or intervention techniques, such as diffuse optical tomography, 2 functional near-infrared spectroscopy, 3 photobiomodulation, 4 etc., requires a quantitative understanding of such complex photon-tissue interactions via computation-based models.
The Monte Carlo (MC) method is widely regarded as the gold-standard for modeling photon propagation in turbid media, 5 including human tissues, due to its accuracy and flexibility. 6 It stochastically solves the general light propagation model-the radiative transfer equation (RTE)-without needing to build large simultaneously linear equations. 7 As an approximation of RTE, the diffusion equation (DE) can be computed more efficiently using finite element-based numerical solvers, 8 and DE is known to yield problematic solutions in regions that contain low-scattering media. 9 In addition to the accuracy and generality, simplicity in implementing MC algorithms compared with other methods has made MC a top choice not only for teaching tissue-optics but also for developing open-source modeling tools.
MC methods have attracted even greater attention in recent years as simulation speed has increased dramatically due to the broad adoptions of massively parallel computing and graphics processing unit (GPU) architectures. The task parallel nature of MC algorithms allows it to be efficiently mapped to the GPU hardware. 10 Current massively parallel MC photon propagation algorithms are capable of handling arbitrary 3D heterogeneous domains and have achieved hundreds-fold speedups compared with traditional serial simulations. [11][12][13][14][15] This breakthrough in the MC algorithm has allowed biophotonics researchers to increasingly use it in routine data analyses, image reconstructions, and hardware parameter optimizations, in addition to its traditional role of providing reference solutions in many biophotonics domains.
A remaining challenge in MC algorithm development is the presence of stochastic noise, which is inherent in the method itself. Because an MC solution is produced by computing the mean behaviors from a large number of photon packets, each consisting of a series of random samplings of the photon scattering/absorption behaviors, creating high-quality MC solutions typically requires simulations of tens to hundreds of millions of photons. This number depends heavily on the domain size, discretization resolution, and tissue optical properties. This translates to longer simulation times because the MC runtime is typically linearly related to the number of simulated photons. From our recent work, 16 a 10-fold increase of photon number typically results in a 10 decibel (dB) signal-to-noise ratio (SNR) improvement in MC solutions, suggesting that MC stochastic noise is largely shot-noise bound. From this prior work, we have also observed that the MC stochastic noise is spatially varying and, in highly scattering/absorbing tissues, exhibits a high dynamic range throughout the simulation domain.
To obtain high-quality simulation results without increasing the number of simulated photons, signal processing techniques have been investigated to remove the stochastic noise introduced by the MC process. This procedure is commonly referred to as denoising. 16,17 In the past, model-based noise-adaptive filters have been proposed to address the spatially varying noise in the radiation dosage estimation context and computer graphics rendering. [18][19][20] However, improvements provided by applying these filtering-based techniques have been small to moderate, creating an equivalent speedup of only three-to fourfold. 16 Recent works on denoising raytraced computer graphics and spatially variant noisy images in the field of computer vision focus mainly on machine-learning (ML)-based denoising methods, more specifically convolutional neural networks (CNNs). 17 Despite their promising performance compared with traditional filters, no attempt has been made, to the best of our knowledge, to adapt denoisers designed for the two-dimensional (2D) low bit-depth image domain to high dynamic range MC fluence maps. 16,21 Our motivation is therefore to develop effective CNN-based denoising techniques, compare them with state-of-the-art denoisers in the context of MC photon simulations, and identify their strengths compared with traditional model-based filtering techniques.
In recent years, the emergence of CNNs has revolutionized many image-processing-centered applications, including pattern recognition, image segmentation, and super-resolution. CNNs have also been explored in image denoising applications, many targeted at removing additive white Gaussian noise (AWGN) from natural images 22 and, more recently, real camera noise. 23,24 Compared with classical approaches, CNNs have also demonstrated impressive adaptiveness to handle spatially varying noise. 25,26 In a supervised setting, given a dataset representative of media encountered in real-life simulations, CNNs have shown to better preserve sharp edges of objects without introducing significant bias, compared with model-based methods. 22,27,28 Finally, due to extensive efforts over the past decade to accelerate CNNs on GPUs, modern implementations of CNN libraries can readily take advantages of GPU hardware to achieve high computational speed compared with traditional methods. Nonetheless, there has not been a systematic study to quantify CNN image denoiser performance over MC photon transport simulation images, either in 2D or 3D domains.
The contributions of this work are the following. First, we develop a simple generative model that uses the Monte Carlo eXtreme (MCX) 12 software to create a synthetic dataset suited for supervised training of an image denoiser, providing ample opportunities for learning its underlying noise structure. Second, we develop and characterize a novel spatial-domain CNN model that cascades DnCNN 26 (an effective global denoiser) and UNet 29 (an effective local denoiser). Third, we adapt and quantitatively compare a range of state-of-the-art image denoising networks, including DnCNN, 26 UNet, 29 DRUNet, 28 deep residual-learning for denoising MC renderings 30 (referred to as ResMCNet hereinafter), and our cascaded denoiser, in the context of denoising 3D MC simulations. We assess these methods using a number of evaluation metrics, including mean-squared error (MSE) and structural similarity index measure (SSIM). For simplicity, other deep-learning (DL)-based denoising methods that do not operate in the spatial domain 31,32 or require specialized knowledge from their target domain 33 are not investigated here and are left for future work. Finally, a range of challenges encountered during the development of our approach are also discussed, providing guidance to future work in this area.

Training Dataset Overview
To train and evaluate CNN denoisers in a supervised fashion, a series of datasets that provided one-to-one mappings between "noisy" and "clean" simulations were generated. The training dataset was created using our MCX software package, 12 in which simulations of a range of configurations with different photon levels were included. The 3D fluence maps generated from the highest number of photons were treated as clean data, and the rest were regarded as noisy. For this work, all configurations were simulated with photon numbers between 10 5 and 10 9 with a 10-fold increment. Simulations with 10 9 photons were selected as the "ground truth," because they provide the closest estimate to the noise-free solutions. Therefore, the CNN denoisers are tasked to learn a mapping between simulations with photon numbers lower than 10 9 to results simulated with 10 9 photons.

Generation of training and validation datasets
To efficiently generate a large and comprehensive corpus of representative MC training data, first a volume generation scheme was designed. In such a scheme, arbitrarily-shaped and -sized polyhedrons and random 3D American standard code for information interchange (ASCII) characters with arbitrary sizes are randomly placed inside a homogeneous background domain with random optical properties. Using combinations of ASCII characters and polyhedrons produces a wide variety of complex shapes, while keeping the data generation process efficient. A similar letterbased random domain generation approach has been previously reported for training networks for fluorescence lifetime imaging. 34 A diagram showing the detailed steps for creating a random simulation domain for generating training data is shown in Fig. 1.
Specifically, a random number (M ¼ 0 to 4) of randomly generated shapes, either in the form of 3D polyhedrons or 3D ASCII letters, are first created as binary masks, with the same size as the target volume. Then, the binary mask is multiplied by a label-a unique identification number assigned to each object-and subsequently accumulated in a final volume, in which voxels marked with the same label belong to the same shape. In the process of accumulation and generation of binary masks for each shape, if two or more objects intersect, this process creates new inclusions for the overlapping regions. We generate all training datasets on 64 × 64 × 64 (in 1-mm 3 isotropic voxels) domains, while the datasets for validation are 128 × 128 × 128 voxels. This allows us to observe the scalability of the networks to volume sizes that are different than the training dataset. A total of 1500 random domains are generated for training and 500 random domains for the "validation." During training, the average global metrics (explained in Sec. 2.4.1) of the model computed over the validation dataset are saved over single epoch intervals. At the end of the training, the model with the best overall metrics is selected as the final result.
To create random 3D polyhedrons, a number of points (N ¼ 4 to 10) are determined on a sphere of random location and radius using the algorithm provided by Deserno. 35 The convexhull of the point set is computed and randomly rotated and translated in 3D. This convex-hull is subsequently rasterized into a binary mask.
For ASCII character inclusions, first, a random character in either lower or upper cases of English alphabet is selected. A random font size is chosen from a specified range, and the letter is rendered/rasterized in a 2D image with a random rotation angle and position. This binary 2D mask is further stacked with a random thickness to form a 3D inclusion. Finally, a 3D random rotation/translation is applied to the 3D ASCII character inclusion.
After generating a random volume, a random simulation configuration is generated to enable simulations with MCX. This includes determining the optical properties, including the absorption (μ a ), scattering (μ s ) coefficients, anisotropy (g), and refractive index (n), for each of the labels inside the generated volume, as well as the light source position and launch direction for the simulation. For the training and validation datasets, only isotropic sources are used for simplicity. The source is randomly positioned inside the domain.
The random optical properties are determined in ranges relevant to those of biological tissues, including (1) μ a ¼ jNð0.01; 0.05Þj mm −1 , where Nðμ; σÞ is a normal distribution with mean μ and standard deviation σ; (2) g is a uniform random variable between 0.9 and 1, where the reduced scattering coefficient μ 0 s ¼ jNð1; 1Þj mm −1 ; and 4) n is a uniformly distributed random variable between 1 and 10. For all data, we simulate the continuous-wave fluence for a time-gate length randomly selected between 0.1 and 1 ns. Each simulation uses a random seed. In Fig. 2, we show a number of image slices (log-10 scale) from 3D simulation samples ranging from homogeneous domains to heterogeneous domains containing multiple polyhedral or letter-shaped inclusions.

Data augmentation
To increase the diversity of the generated dataset and avoid overfitting, data augmentation 36 was used. Our data augmentation consisted of 90-deg rotation and flipping. Each transformation was applied independently over a randomly selected axis. Transforms were identically applied to both inputs and labels of the training data. Both transforms were randomly selected and applied, with a probability of 0.7. This on-the-fly strategy multiplied the data encountered by the models during training by 256 without performing any time-consuming MC simulation.

Test datasets
Three previously used standard benchmarks, 16 (B1) a 100 × 100 × 100 mm 3 homogeneous cube with a 1-mm voxel size, (B2) the same cubic domain with a 40 × 40 × 40 mm 3 cubic absorber, and (B3) the same cubic domain with a refractive inclusion, were employed to characterize and compare the performance of various denoising methods. The optical properties for the background medium, the absorbing and refractive inclusions, can be found in Sec. 3 of our previous work. 16 Each of the benchmarks was simulated with 100 repetitions using different random seeds. Additionally, the Colin27 12,37 atlas (B4), Digimouse 38 atlas (B5), and University of Southern California (USC) 19.5 39 atlas (B6) from the Neurodevelopmental MRI database 40 were selected as examples of complex simulation domains to test our trained MC denoisers. In addition, we also included a benchmark (B7) containing a ball-lens to test the performance of our denoisers in low-albedo media. In this benchmark, a cubic domain of 100 × 100 × 100 grid of 0.1 mm 3 isotropic voxels is filled with medium of μ a ¼ 0.01∕mm, μ s ¼ 1∕mm, g ¼ 0.95, and n ¼ 1. A spherical-lens of 2-mm radius is placed in the center of the domain and filled with a medium of μ a ¼ 0.01∕mm, μ s ¼ 1∕mm, g ¼ 0.9, and n ¼ 1.4. A pencil beam pointing toward the þz axis is located at (4, 5, 0) mm. The domain volume is pre-processed to utilize the splitvoxel MC algorithm, 41 which can accurately handle curved media boundaries rasterized using a voxelated domain.

Pre-Processing of Monte Carlo Data
Many of the reported DL denoising techniques were developed to process natural images of limited bit-depth that usually do not present the high dynamic range as in MC fluence maps. To allow CNNs to better recognize and process unique MC image features and avoid difficulties due to limited precision, we applied the following transformation to the fluence images before training or inference: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 2 0 7 where x is the MC fluence map, c is a user-defined constant, and the output y serves as the input to the CNN. This transformation serves two purposes. First, it compresses the floating-point fluence values to a limited range while equalizing image features across the domain. Second, it compensates for the exponential decay of light in lossy media and reveals image contracts that are relevant to the shapes/locations of the inclusions, assisting the CNN to learn the features and mappings. The addition of 1 in Eq. (1) ensures that tðxÞ does not contain negative values. An inverse transform t −1 ðy 0 Þ ¼ ðe y 0 − 1Þ∕c is applied to the output of the CNN (y 0 ) to undo the effect of this transform. Moreover, when training a CNN on 8-bit natural image data, a common practice is to divide the pixel values by the maximum value possible (i.e., 255) to normalize the data. From our tests, applying such an operation on floating-point fluence maps resulted in unstable training; therefore our training data were not normalized.
Additionally, due to limited data precision, we noticed that all tested CNN denoisers exhibit reduced denoising image quality when processing voxel values (before log-transformation) that are smaller than an empirical threshold of 0.03. To address this issue and permit a wider input dynamic range, two separate copies of the fluence maps were denoised during inference-the first copy was denoised with c set to 1 and the second one with c set to 10 7 . The final image is obtained by merging both denoised outputs: voxels that originally had fluence values larger than 0.03 retrieve the denoised values from the first output and the rest are obtained from the second output. This variable-gain approach allowed us to process MC fluence images containing both high and low floating-point values.

Cascaded MC Denoising Network that Combines DnCNN and UNet Networks
In this work, we designed a cascaded CNN denoiser, as shown in Fig. 3, specifically optimized for denoising our 3D MC fluence maps by combining two existing CNN denoisers: a DnCNN denoiser is known to be effective for removing global or spatially invariant noise, especially AWGN, without any prior information, 26 whereas a UNet denoiser is known to remove local noise that is spatially variant. 28,29 Therefore, in our cascaded DnCNN/UNet architecture, referred to as "cascade" hereinafter, the CNN first learns the global noise of an MC fluence image and attempts to remove it. The remaining spatially variant noise can then be captured and removed using a UNet. In both stages, the noise is learned in the residual space, meaning that, instead of mapping a noisy input to a clean output directly, the network maps the noisy input to a noise map and then subtracts it from the input to extract the clean image. We want to mention that cascaded denoisers similar to the above design have been proposed for processing real-world images. 25,42 In these works, a model-based filter (BM3D) serves as the global denoiser to provide an improved prior to a CNN-based local denoiser. In comparison, our method utilizes CNN denoisers for both global and local denoising stages, making it possible to train and accelerate fully on GPUs while automatically adapting to varying levels of noise.

Global performance metrics
The global resemblance between the denoised volume and the ground truth (in this case, simulations with 10 9 photons) can be used to measure the performance of a denoiser. A number of metrics measuring such a similarity have been used by others to evaluate image restoration networks or measure convergence. 21,26,43,44 Typically, these metrics are defined for 2D images; in this work, we extended the definitions to apply to 3D fluence maps. Fig. 3 Overview of the cascaded DnCNN + UNet architecture. Each block in the dashed squares represents a group of CNN layers that are applied sequentially. The number on the square block indicates the number of channels for the respective output tensor. Conv3D, TransConv3D, and BN stand for 3D convolution, 3D transposed convolution, and batch normalization layers, respectively. PyTorch function log 1pðcxÞ is a stable implementation of function lnðcx þ 1Þ.
The most commonly used objective functions for denoising networks are the mean least squared error (L 2 ) and mean absolute error (L 1 ): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 7 1 1 where K is the number of noisy-clean fluence map pairs sampled from the dataset, referred to as the "mini-batch" size; θ contains all parameters of the network; F denotes the network itself; n is either 1 or 2; and ðx i ; y i Þ denotes the i'th noisy-clean pair of data in the mini-batch. These error metrics are widely used in supervised denoising networks, including the DnCNN, DRUNet, and ResMCNet models, as well as several other studies. 25,26,28,30,42 L 1 and L 2 may have different convergence properties. 43 The L 1 loss has gained more popularity in the DL community due to its good performance and low computational costs. 30,43 For this work, however, to penalize large errors more, the L 2 loss was used instead to train the networks.
In contrast to L n distances, SSIM 45 provides a perceptually motivated measure that emulates human visual perception for images. The SSIM for a pixel in an image is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 5 4 4 SSIMðpÞ where μ x and σ x are the mean and standard deviation of the image x, respectively, and σ xy is the co-variance of images x and y. The statistics are calculated locally by convolving both volumes with a 2D Gaussian filter with σ G ¼ 5. Small constants C 1 and C 2 are used to avoid division by zero. The SSIM value of two images is the average SSIM computed across all pixels, with a value of 1 suggesting that the two images are identical, and a value of 0 suggesting that the two images are not correlated. This definition can also be applied to 3D fluence maps using a 3D Gaussian kernel to calculate neighborhood statistics. Another metric, peak signal-to-noise ratio (PSNR), measures the ratio between the maximum power of a signal and the power of the noise. 46 The PSNR for two volumes x and y is expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 3 7 7 PSNRðx; yÞ ¼ 20 log 10 I max ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kx − yk 2 p : (4) Larger PSNR values indicate smaller L 2 distances between volumes. The I max value is the maximum value that a voxel can have in a fluence map after the transformation in Eq. (1). Therefore, in this work, we set I max to 40.

Local performance metrics
A number of locally (voxel-bound) defined performance metrics have been used in our previous MC denoising work. 16 The SNR of the denoised volumes for each voxel measures the efficacy of the denoiser of spatially adaptive noise. For a simulation running k photons, we first run multiple (N ¼ 100) independently seeded MC simulations and compute SNR in dB with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 2 1 3 SNR k ðrÞ ¼ 20 log 10 μ k ðrÞ σ k ðrÞ ; where μ k and σ k are the mean and standard deviation of voxel values at location r across all repetitions, respectively. The average SNR difference before and after applying the denoising filter, ΔSNR, is subsequently calculated along selected regions of interest.
Our previous work 16 suggests that the noise in MC images largely follows the shot-noise model; therefore, increasing the simulated photon number by a factor of 10 results in ∼10 dB improvement in SNR on average. We have previously proposed a photon number multiplier 16 M F to measure equivalent acceleration using the average SNR improvement ΔSNR: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 8 6 M F ¼ 10 ΔSNR 10 : For example, a ΔSNR ¼ 20 dB gives M F ¼ 100, suggesting that the denoised result is equivalent to a simulation with 100 times the originally simulated photon number, which is equivalent to accelerating the simulation by a factor of 100 if the denoising run-time is ignored.

BM4D and ANLM
Block-matching four-dimensional collaborative filtering (BM4D) and our GPU-accelerated adaptive non-local means (ANLM) 16 are used as representative state-of-the-art model-based denoisers and used to compare against CNN-based denoisers. For BM4D, a Python interface developed based on the filter described by Mäkinen et al. 47 was used, whereas for the ANLM filter, a MATLAB function developed previously by our group 16 was used.

CNN training details
All CNN denoising networks were re-implemented for handling 3D data using the open-source DL framework, PyTorch. 48 For most of the studied CNN denoisers, our implementations largely follow their originally published specifications but replacing the 2D layers with their 3D variants. Small adjustments were made. For UNet, for example, 3D batch normalization (BN) layers were introduced in between the 3D convolution, the convolution transpose, and the pooling layers to address the covariance shift problem. 49 Additionally, we simplified ResMCNet by removing the auxiliary features needed for computer graphics renderings purposes, making the kernel size of the first layer 3 instead of 7.
All networks in this study were trained for 1500 epochs on a single NVIDIA DGX node equipped with eight NVIDIA A100 GPUs, each with 40 GB of memory and NVLink 2.0 connection. Leveraging the PyTorch scaling wrapper, PyTorch Lightning 50 was used to simplify the implementation process. We need high-performance hardware because a forward propagation of the CNN for a 64 × 64 × 64 voxelated volume requires around 6 GB of GPU memory; to use a batch size of 4 per GPU (i.e., processing four data pairs in parallel), at least 24 GB of memory is necessary. Furthermore, using all 8 GPUs in parallel combined with the high-speed NVLink connection reduces the average training time from 10 days (on a single A100 GPU) to 24 h for each network tested-the cascade and DRUNet usually require longer training times compared with those of DnCNN and UNet.
The networks were all trained using the Adam with weight decay regularization optimizer, 51 with a weight decay of 0.0001 for the parameters in all layers, except for the BN parameters and bias parameters. The learning rate was scheduled with a cosine annealing learning rate, 52 using 1000 linear warm-up mini-batch iterations to added learning stability. 53 A batch-size of four per GPU was selected to maximize the effective use of GPU memory resources. The base learning rate was set to 0.0001. The gradient clipping value was set to 2 for BN layers and 1 for other layers to avoid exploding gradients and faster training. 54 The optimization, data augmentation, and configuration sections of the codebase for this work were inspired by the open-source PyTorch Connectomics package 55 for easier prototyping of the trained models.

Denoising Performance
In Fig. 4, we visually compare the fluence maps before and after denoising for each tested denoiser and photon number (10 5 to 10 8 ) for three standard benchmarks 16 (B1, B2, and B3). Table 1 summarizes the global metrics derived from the outputs of each denoiser; computed local metrics including mean ΔSNR and M F are given in Table 3. Each entry in both tables is averaged from 100 independently seeded repeated simulations. In both tables, the best-performing metrics are highlighted in bold. Similarly, a visual comparison between those from more complex domains, including Colin27, Digimouse, USC-19.5 atlases, and the ball-lens benchmark, are shown in  Table 2. Due to limited space, in Fig. 5, we only show representative images with 10 5 and 10 7 photons, and we removed DnCNN and BM4D due to their relatively poor performance. For the same reason, BM4D global metric results were removed from Table 2.
From the denoised images shown in Fig. 4, we can first confirm that all CNN-based denoisers show noise-adaptive capability similar to ANLM and BM4D-they apply a higher level of smoothing in noisy areas within low-photon regions and apply little smoothing in areas with

Table 1
Average global metrics derived from three basic benchmarks: (B1) a homogeneous cube, the same cube with (B2) as absorption, and (B3) refractive index inclusion; each data point was averaged over 100 repetitions. The best performing models are highlighted in bold.   Fig. 4(c), we can also observe that all CNN denoisers show edge-preservation capability, again similar to ANLM and BM4D. Both noise-adaptiveness and edge-preservation are considered desirable for an MC denoiser. 16 Because all CNN networks were trained on images of 64 × 64 × 64 voxels whereas all three benchmarks shown in Fig. 4 are 100 × 100 × 100 voxel domains, these results clearly suggest that our trained networks can be directly applied to image domain sizes that are different from the training domain size. By visually inspecting and comparing the denoised images in Figs. 4 and 5, we observed that all CNN-based methods appear to achieve significantly better results compared with modelbased denoising methods (BM4D and GPU ANLM); such difference is even more pronounced in low-photon simulations (10 5 and 10 6 photons). Although the CNN denoisers were trained on shapes with less complexity, the images in Fig. 5 indicate that they are clearly capable of denoising novel structures that are significantly complex, yielding results that are close to the respective ground truth images. However, we also observe that the denoiser's ability to recover fluence maps varies depending on the photon level in the input data-in areas where photons are sparse, the denoisers understandably create distortions that deviate from the ground truth. This can be seen clearly in the results of the complex benchmarks B4 to B7 at 10 5 . Nevertheless, these distorted recovered areas are still significantly better estimates than the input in the same area without denoising.
To confirm that CNN denoisers can produce unbiased images, the means and SNRs from benchmarks B1, B2, and B3 along the line x ¼ 50 and y ¼ 50 were calculated and plotted in

Table 3
Overall average SNR improvements (ΔSNR all in dB) and those (ΔSNR eff ) in the effective region (where ΔSNR > 0.5 dB) as well as the photon number multipliers (M F ) in the three basic benchmarks (B1-B3).  Fig. 6. For brevity, we only report the results from the cascade network as representative of all CNN methods in this plot. These plots confirm that the cascade method does not alter the mean fluence of the simulations over the plotted cross section, while providing a consistent SNR improvement across a wide range of photon numbers. It also demonstrates the adaptiveness of CNN denoisers in that SNR improvement starts to decline in areas with high fluence value (thus lower noise due to shot-noise). The ∼12-dB SNR improvement shown by denoising simulations with 10 9 photons (purple dotted lines over purple solid lines in the SNR plots) indicate that the cascade denoiser is capable of further enhancing image quality even it was not trained using simulations with more than 10 9 photons. Such an SNR improvement is not as high as that reported from low-photon simulations, yet it is still significantly higher than the best SNR improvement produced using the GPU ANLM denoiser (dashed lines) of all tested photon numbers.
Our earlier observation that most CNN-based denoisers outperform model-based denoisers (GPU ANLM and BM4D) is also strongly evident by both the global metrics reported in Table 1 (a) (b) (c) Fig. 6 Plots of the means (left) and SNRs (right) before (solid) and after denoising using cascade network (dotted) and GPU-ANLM (dashed) in 3 benchmarks (a) B1, (b) B2, and (c) B3 along a cross section. and the local metrics reported in Table 3. Among all tested CNN filters, the cascade network offers the highest performance in all tests with 10 5 photons and comes close to the best performer, ResMCNet, among the 10 6 test sets. Among the 10 7 photon levels, DRUNet is a strong performer, with ResMCNet and cascade coming close to or surpassing it in some cases. Among the real-world complex domain benchmarks shown in Table 2, cascade reports the best performance in almost all cases, with UNet performing slightly better on USC-195 with 10 5 photons and ResMCNet giving better SSIM results.
From Table 3, we can observe that all CNN-based denoisers appear to offer a five-to eightfold improvement in SNR enhancement compared with our previously reported model-based GPU ANLM filter 16 ; our cascade network reports an overall SNR improvement between 14 to 19 dB across different benchmarks and photon numbers. This is equivalent to running 25× to 35× more photons in heterogeneous domains, and nearly 80× more photons for the homogeneous benchmark (B1). In other words, applying our cascade network for an MC solution with 10 5 photons can obtain a result that is equivalent to running ∼2.5 × 10 6 photons. In fact, except for DnCNN, the majority of our tested CNN-based denoisers can achieve a similar level of performance.

Assessing Equivalent Speedup Enabled by Image Denoising
In Table 4, we report the average runtimes (in seconds) of MC simulation and denoising (i.e., inference for CNN denoisers). Each test case runs on a single NVIDIA A100 with 40 GBs of memory with over 100 trials, and the time needed to transfer data between the host and the GPU is included. As we mentioned in Sec. 2.2, to obtain every denoised image, we apply CNN inference twice to handle the high dynamic range in the input data. Table 3 suggests that, on average, about a 20 to 30 photon multiplier (M F ) is to be expected for most CNN denoisers, meaning the denoised simulations will have 20 to 30 times more photons than its input. Therefore, our goal is to identify cases in which the sum of the runtime of the baseline MC simulation running on N photons, T MC ðNÞ, and that of the denoiser (T f ) is shorter than an MC simulation running M F × N photons, i.e., T MC ðNÞ þ T f < T MC ðM F × NÞ. Due to space limitations, we are unable to list all combinations of simulations that satisfy the above condition. However, our general observations include the following: (1) the CNN inference runtime is independent of the number of simulated photons; (2) DnCNN is typically faster than other CNN denoisers, but it also has the poorest performance among them from Table 3; (3) the larger the domain size is, the longer it takes for CNN denoisers to run; and (4) generally speaking, applying CNN denoisers to simulations with 10 7 photons or above can result in a significant reduction of total runtime. In our previous work, 16 we had also concluded that 10 7 photon is a general threshold for GPU-ANLM to be effective; however, from the runtime data reported here using NVIDIA A100 GPUs, GPU-ANLM appears to also benefit simulations with 10 6 photons, likely due to the high computing speed of the GPU. Nonetheless, comparing to most tested CNN denoisers, the GPU-ANLM denoiser offers dramatically less equivalent acceleration despite its fast speed.

Conclusion
In summary, we have developed a framework for applying state-of-the-art DL methods for denoising 3D images of MC photon simulations in turbid media. A list of supervised CNN denoisers, including DnCNN, UNet, ResMCNet, and DRUNet, were implemented, extended for processing 3D data, and tested for denoising MC outputs. In addition, we have developed a customized cascaded DnCNN/UNet denoiser combining the global-noise removal capability of DnCNN and local-noise removal capability of UNet. All developed MC denoising networks were trained using GPU accelerated MCX simulations of random domains to learn the underlying noise from MC outputs at a range of photon numbers. A simple yet effective synthetic training data generation approach was developed to produce complex simulation domains with random inclusions made of 3D polyhedral and ASCII characters with random optical properties and simulation parameters. In addition to following current best practices of contemporary CNN and DL development, we have also specifically fine-tuned and customized our MC denoisers to better handle the unique challenges arising in denoising 3D MC data. For example, to handle the high dynamic range in MC fluene maps using CNNs, a reversible log-mapping scheme was applied to each volume before being fed to the models. In addition, we have also applied inference twice and combined the results to further enhance the dynamic range of the input data. All reported CNN MC denoisers have been implemented in the Python programming language using the PyTorch framework, with both source codes and training data freely available to the community as open-source software.
To evaluate the efficacy of these proposed CNN denoisers, we have constructed seven standard benchmarks-three simple domains and four complex ones-from which we have derived and reported both global performance metrics (such as SSIM and PSNR) and local performance metrics (such as ΔSNR and M F ). From our results, all tested CNN-based denoisers offered significantly improved image quality compared with model-based image denoisers such as GPU-ANLM and BM4D in this particular application. Overall, most CNN denoisers provide a 10-to 20-dB SNR improvement on average, equivalent to running 10-to 100-fold more photons. Among these CNN denoisers, our proposed cascade network outperformed most of the state-of-the-art spatial domain denoising architectures and yielded the best image quality for low-photon simulations with 10 5 and 10 6 photons. Its performance is on-par with or only slightly inferior to DRUNet in high-photon simulations (10 7 photon) in simple domain tests. For all benchmarks involving real-world complex domains, the cascade network yielded the highest global metrics in nearly all tests. In comparison, some of the most effective model-based image denoisers such as the GPU-ANLM filter that we proposed previously 16 only yielded 3 to 4 dB improvement, despite being relatively fast to compute. It is worth noting that the cascade network yielded an impressive 80-fold equivalent speedup when processing low-image-feature simulations such as a homogeneous domain.
From our tests, CNN denoisers demonstrate superior scalability to input data sizes and input image qualities. Although our training data were produced on a 64 × 64 × 64 voxelated space with relatively simple shapes, all tested CNN denoisers show no difficulty in handling images of larger sizes or significantly more complex inclusions. Our cascade network also reported a 12-dB average SNR improvement when being applied to denoise baseline simulations with 10 9 photons-the level of photon number that was used as the ground truth for training. This suggest that these CNN denoising architectures may not be strictly limited by the quality of the data that they are trained on.
From our results on runtimes, most CNN denoiser inference (including two passes) time ranges between less than a second to a dozen seconds, regardless of the input data quality. We concluded that, to yield an overall shorter total runtime, applying CNN denoisers to processing MC images generated from 10 7 photons or more can generally lead to significantly improved computational efficiency.
One of the limitations of the current work is the relatively long training time. To train each denoising network using our synthetic dataset of 1500 random domains (each with five photon number levels with multiple rotated views) requires on average a full day (24 h) if running on a high-end 8-GPU server with large-memory NVIDIA A100 GPUs (40 GB memory allows for using a batch-size of 4 for acceleration). If running on a single GPU node, we anticipate that the required training time is around 10 to 12 days on a single A100 GPU and even longer for lowmemory GPUs. Experimenting with the number of layers in each model to reduce the number of intermediate tensors while retaining the performance benefits reported in this work, as well as the development of new and significantly more compact DL-based denoisers, will be the focus of our future work. Moreover, some of the training parameters were determined empirically and deserve further optimization. For example, we trained the networks over 64 × 64 × 64 domains. It will be significantly faster if we can reduce the training data size while still retaining the scalability to arbitrarily sized domains. Additionally, the landscapes of CNN architecture and denoising networks are constantly being updated and improved over the past few years. We cannot exhaust all emerging CNN denoisers and would be happy to extend this work with newer and more effective CNN denoising architectures in the future.
To conclude, we strongly believe that investigating high-performance image denoising techniques offers a new direction for researchers seeking for the next major breakthrough in speed acceleration for MC simulations. DL-and CNN-based image denoising techniques have demonstrated impressive capabilities compared with the more traditional model-based denoising methods and have yielded notable image quality enhancement that is equivalent to running 10 to 50 times more photons, which can be directly translated to a 10-to 50-fold speedup, in most of our tested benchmarks. Our cascade denoising network even reported a nearly 80-fold equivalent speedup when denoising homogeneous domain results-a level of acceleration that we were only able to witness when migrating MC from single-threaded computing to massively parallel hardware over a decade ago. [11][12][13] With the combination of advanced image processing methods and new simulation techniques, we anticipate that MC will play an increasingly important role in today's biomedical optics data analysis and instrument development.

Disclosures
No conflicts of interest, financial or otherwise, are declared by the authors.