\subsubsection{Synthetic Setting}
\label{sec:denoising_synth_discussion}

\begin{figure}[t]
\centering
\includegraphics[width=\linewidth]{figures/denoising_imgs/synthetic/synth_comparison.png}
\caption{Denoising visual and quantitative comparison of synthetic data: (a) The results of data setting (100, 2). (b) The results of data setting (2000, 50). Note that in these two settings our method (Ours) does not apply data augmentation.}
\vspace{-1.5em}
\label{denoised_results_synthetic_setting}
\end{figure}

Using the synthetic data with ground truth in the setting (100, 2) shown in \autoref{denoised_results_synthetic_setting} (a), we have very limited amount of projections (i.~e., only 2) and a low dose which leads to a low \gls*{snr}. In settings like this where data are scarce and full of noise, a natural choice is to average the two available projections. Although the denoised result is not close to the noise-free reference, our method outperforms the simple average of the two available projections in terms of \gls*{ssim}~\cite{ssim}, \gls*{mse}, \gls*{psnr} and \gls*{ffl}. In the very high dose cases (i.~e., {700, 1000, 1500, 2000}), we show consistently better results than simple averaging in terms of \gls*{ffl}. For example, the quantitative comparison in \autoref{denoised_results_synthetic_setting} (b) shows that the \gls*{ffl} of our method is roughly 4 times better than simple averaging. The quantitative comparisons of these two figures show that our physically-based Noise2Noise-like model is effective.

To compare with the original Noise2Noise~\cite{noise2noise} model and Topaz-Denoise~\cite{topaz_denoiser}, we present our results from the setting (100, 50) denoted as Ours-50 in \autoref{synthetic_denoising_comparison}. The figure also shows the results of a trained Noise2Noise model, \texttt{Topaz P} and \texttt{Topaz T}. The Noise2Noise model does not learn to denoise even if we have a large amount of available data (i.~e., 5000). Therefore, we further tested the performance of the Topaz model with our synthetic dataset. The results show that the pretrained model (i.~e., \texttt{Topaz P}) generalizes well by filtering out high-frequency noise and preserves the particles, but it is worse than our method because the boundaries of particles are blurred. Our trained Topaz model (i.~e., \texttt{Topaz T}) filters out most of the noise and recognizes the positions of the particles. However, it does not preserve the shape of the particles, which makes it unusable. On the contrary, our method can work in the setting with less data (i.~e., (100, 50)), which shows that our method is advantageous over training a Noise2Noise model~\cite{noise2noise} and a Topaz-Denoise~model~\cite{topaz_denoiser} from scratch.

It is noteworthy that the Noise2Noise model and \texttt{Topaz T} are trained with data settings (100, 5000) and (100, 4500) respectively. It is not realistic nor viable to obtain such large datasets in real-world scenarios, which means these two models are not usable. However, \texttt{Topaz P} generalizes well in our setting, so we use it for data augmentation in the real-world scenario to generate more data for our denoising method.

For more results of our method under synthetic settings, we refer the reader to the supplementary material.

\begin{figure}[t]
\begin{center}
\includegraphics[width=\linewidth]{figures/denoising_imgs/synthetic/Synthetic-comparison-3.png}
\caption{Comparison of denoising the synthetic data: images in the top row show original noisy projection (NP) and clean ground truth projection (GT) and our denoising result, in bottom row are denoising results using the trained Noise2Noise result (Noise2Noise), pretrained Topaz-Denoise~(\texttt{Topaz P}), and Topaz-Denoise~trained on our noisy projections (\texttt{Topaz T}).}
\vspace{-2.5em}
\label{synthetic_denoising_comparison}
\end{center}
\end{figure}

\begin{table*}[t]
\centering\caption{SSIM, PSNR and FFL comparison} \label{SSIM_PSNR_FFL_comparison}
\resizebox{\textwidth}{!}{%
\begin{tabular}{cccccc|ccccc|ccccc}
\multicolumn{6}{c}{\textbf{SSIM}} & \multicolumn{5}{c}{\textbf{PSNR}} & \multicolumn{5}{c}{\textbf{FFL}}\\
\toprule
\multicolumn{1}{c}{}      & \texttt{Noisy}    & \texttt{Topaz}    & \texttt{Ours(T)}        & \texttt{Ours(N)}         & \texttt{Ours(T+N)}      & \texttt{Noisy}    & \texttt{Topaz}    & \texttt{Ours(T)}        & \texttt{Ours(N)}        & \texttt{Ours(T+N)}      & \texttt{Noisy}    & \texttt{Topaz}    & \texttt{Ours(T)}        & \texttt{Ours(N)}   & \texttt{Ours(T+N)}      \\
\midrule
\multicolumn{1}{r|}{\texttt{Noisy}} & 1.0000   & 0.6017   & 0.4139   & 0.7395   & \textbf{0.6547}   & 361.200  & 65.8250  & 65.8270  & 70.1500  & \textbf{66.1000}  & 0.0000   & 0.0096   & 0.0082   & 0.0032   & \textbf{0.0095}   \\
\multicolumn{1}{r|}{\texttt{Topaz}} & 0.6017   & 1.0000   & 0.6894   & 0.7184   & \textbf{0.8651}   & 65.8250  & 361.200  & 72.2000  & 68.9700  & \textbf{75.0400}  & 0.0096   & 0.0000   & 0.0012   & 0.0045   & \textbf{0.0008}   \\
\multicolumn{1}{r|}{\texttt{Ours(T)}}     & 0.4139   & 0.6894   & 1.0000   & 0.7797   & 0.8001   & 65.8270  & 72.2000  & 361.200  & 71.6200  & 76.2800  & 0.0082   & 0.0012   & 0.0000   & 0.0020   & 0.0001   \\
\multicolumn{1}{r|}{\texttt{Ours(N)}}     & 0.7395   & 0.7184   & 0.7797   & 1.0000   & 0.9098   & 70.1500  & 68.9700  & 71.6200  & 361.200  & 71.6700  & 0.0032   & 0.0045   & 0.0020   & 0.0000   & 0.0029   \\
\multicolumn{1}{r|}{\texttt{Ours(T+N)}}   & 0.6547   & 0.8651   & 0.8001   & 0.9098   & 1.0000   & 66.1000  & 75.0400  & 76.2800  & 71.6700  & 361.200  & 0.0095   & 0.0008   & 0.0001   & 0.0029   & 0.0000   \\
\bottomrule
\end{tabular}
}
\end{table*}

\begin{figure*}[t]
\begin{center}
\includegraphics[width=\textwidth]{figures/reconstruction-and-dvr-comparisson-small.png}
\caption{Tomographic reconstructions from the tilt series (left) and DVR renderings (right).}
\vspace{-1em}
\label{fig:reconstruction-and-dvr-comparison}
\end{center}
\end{figure*}

In the experiments using real micrographs, we use the tilt series from Yao \textit{et~al.}~\cite{Sai2020:CryoET-Pipeline}, which show SARS-CoV-2 virions. The dose in the dataset is 320 electrons per nm\textsuperscript{2} per projection. Therefore, without data augmentations, the data setting is (320, 1). As discussed in Section~\ref{denoising_section}, our method cannot work in a scenario with only single projection, so in our experiments, we apply a number of data augmentations. In the setting \texttt{Ours(T)}, the data augmentation with \texttt{Topaz P} is used, yielding one more projection, which means the data setting is (320, 2). In the setting denoted as \texttt{Ours(N)}, neighboring projections of a projection are aligned and used, so the data setting is (320, 3). In the setting \texttt{Ours(T+N)}, \texttt{Topaz P} and aligned neighboring projections are used, so the data setting is (320, 6). The case of \texttt{Ours(T+N)} matches exactly our \autoref{denoising_optimization} while in the case of \texttt{Ours(T)}, $P_n=\{P\}$, and in the case of \texttt{Ours(N)}, preprocessing with Topaz-Denoise~is not used (i.~e., $P_{\train}=\normalize(P_n)$). In the real setting, we do not have noise-free ground truth projections to compare against, so we first present results of the qualitative comparison.

On the top of \autoref{denoising_with_real_data}, we can see \texttt{Topaz} has block-like artifacts and \texttt{Noisy} contains more noise while others have minor visual difference. However, if we zoom-in the images, as shown in the bottom of ~\autoref{denoising_with_real_data}, we can see (1) \texttt{Topaz} has more blurring artifacts, (2) \texttt{Noisy} has the most noise, and (3) \texttt{Ours(T+N)} is better than \texttt{Ours(N)} and \texttt{Ours(T)} in that it has less high-frequency noise while preserving clear boundaries. This makes our method \texttt{Ours(T+N)} the best among these five because of its good balance between denoising and preserving boundaries from the view of quality.

\begin{figure}[t]
\begin{center}
\includegraphics[width=\linewidth]{figures/RealDataDenoising-composite.png}
\caption{\texttt{Topaz}: the real projection denoised by \texttt{Topaz P}; \texttt{Noisy}: the real noisy projection; Top: Full views; Middle: Zoom-in of the first row. In the bottom row is another example of our denoising method, \texttt{Ours(T+N)} contains less noise while preserving boundaries.}
\vspace{-2em}
\label{denoising_with_real_data}
\end{center}
\end{figure}

We also present quantitative results in Table~\ref{SSIM_PSNR_FFL_comparison}. Based on the bottom images from \autoref{denoising_with_real_data}, we can conclude that \texttt{Topaz} contains the least noise while \texttt{Noisy} contains the most noise. With this conclusion, the values in bold in the table show that (1) \texttt{Ours(T+N)} is the closest one to \texttt{Topaz} and (2) \texttt{Ours(T+N)} has a moderate amount of noise compared to \texttt{Ours(T)} and \texttt{Ours(N)}. This means \texttt{Ours(T+N)} is a good compromise between removing noise and preserving details, echoing with our qualitative comparison.

In addition to the example on the top two rows of \autoref{denoising_with_real_data}, we give another example of denoising real data. From the last row in \autoref{denoising_with_real_data}, we can draw similar conclusions. \texttt{Topaz} has the least noise, but compromised boundaries while \texttt{Ours(T+N)} has a moderate amount of noise but also clearer boundaries and more visibility of spikes.

The application of such denoising are also useful for the tomographic reconstruction. In left part of  \autoref{fig:reconstruction-and-dvr-comparison} we show the mid-slice from the volumes reconstructed using SIRT reconstruction technique from IMOD~\cite{kremer1996computer} from original noisy projections (left), from denoised projections using \texttt{Topaz} (middle), and from denoised projections using our approach \texttt{Ours(T+N)} (right). One can clearly see that \texttt{Topaz} and our approach are better that reconstruction from the original noisy projections. Moreover, our approach retains more details which are smudged out by \texttt{Topaz}.

Additionally, in the right part of \autoref{fig:reconstruction-and-dvr-comparison}, we show \gls*{dvr} rendering of the same volumes using the same rendering parameters apart from the threshold, which was manually adapted for each case. One can see that \texttt{Topaz} cleans the volume extensively, but also removes parts of the virions and removes the high frequency details from the structures. This is not the case for results using our approach. While there is more noise present in the volume, the high-frequency details are retained. Both approaches are clearly better choice than rendering the original noisy volume.

\subsection*{More Comparisons of Denoising Synthetic Projections}
Here we present more comparisons of results of different synthetic settings to show the effectiveness and usefulness of our denoising method.

When the number of samples is small, for example, 2 or 5. Our results are close to the averages of the datasets visually. Yet, our results are better in terms of quantitative metrics as shown in \autoref{fig:num_2_dose_200} and \autoref{fig:num_5_dose_200} in which the results of data settings of (200, 2) and (200, 5) are presented.

When the dose is high (i.~e., greater than 700 electrons per nm\textsuperscript{2} per projection), our method consistently outperforms the averages of the datasets by a large margin in terms of \gls*{ffl}, which indicates our denoising method is much more correct in frequency domain. Here we present the results of the settings (700, 10) and (1500, 10). In \autoref{fig:num_10_dose_700} we can see the \gls*{ffl} is almost 50\% better than the average of the dataset and in \autoref{fig:num_10_dose_1500} the margin is more than 400\%.

Along with the results in \autoref{denoised_results_synthetic_setting}, we can conclude that despite the margins are different with different data settings, our method is advantageous over simple averaging consistently. This is especially beneficial when the data are too scarce and too noisy to train neural networks.

\begin{figure}[H]
\begin{center}
\includegraphics[trim={3.5cm 2cm 3.5cm 2cm},clip=true, width=\linewidth]{figures/denoising_imgs/synthetic/comparison_ref_sample_num=2_dose=200_normalized_metrics.png}
\caption{Denoising visual and quantitative comparison of synthetic data with data setting (200, 2)}
\label{fig:num_2_dose_200}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[trim={3.5cm 2cm 3.5cm 2cm},clip=true,width=\linewidth]{figures/denoising_imgs/synthetic/comparison_ref_sample_num=5_dose=200_normalized_metrics.png}
\caption{Denoising visual and quantitative comparison of synthetic data with data setting (200, 5)}
\label{fig:num_5_dose_200}
\end{center}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[trim={3.5cm 2cm 3.5cm 2cm},clip=true,width=\linewidth]{figures/denoising_imgs/synthetic/comparison_ref_sample_num=10_dose=700_normalized_metrics_unavg.png}
\caption{Denoising visual and quantitative comparison of synthetic data with data setting (700, 10)}
\label{fig:num_10_dose_700}
\end{center}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[trim={3.5cm 2cm 3.5cm 2cm},clip=true,width=\linewidth]{figures/denoising_imgs/synthetic/comparison_ref_sample_num=10_dose=1500_normalized_metrics_unavg.png}
\caption{Denoising visual and quantitative comparison of synthetic data with data setting (1500, 10)}
\label{fig:num_10_dose_1500}
\end{center}
\end{figure}

\subsection{Denoising}
\label{denoising_related_work}
When denoising \gls*{cryo-em} micrographs, clean references are not available, so many approaches have been proposed to specifically denoise \gls*{cryo-em} micrographs, addressing this problem.

Researchers have proposed various analytical filters for denoising, such as simple averaging, Gaussian filters, and many filters in the frequency space that enhance or damp signal responses. Additionally, physics-based techniques were proposed, such as \gls*{nad}~\cite{NAD}, which takes advantage of a diffusion model to denoise 3D volumes. For a comprehensive review on filters used in \gls*{cryo-em}, we refer the reader to the review~\cite{cryo_em_filtering}.

In the era of \gls*{dl}, many neural networks are proposed to denoise \gls*{cryo-em} micrographs. Noise2Noise~\cite{noise2noise} is a general denoising approach that does not require clean data and has been adapted in CryoCARE~\cite{cryo-care} and Topaz-Denoise~\cite{topaz_denoiser}. In CryoCARE~\cite{cryo-care}, a neural network is trained as a denoiser with multiple micrographs that come from multiple consecutive exposures or multiple aligned micrographs in a tilt series or movie frames captured during a long exposure. In Topaz-Denoise~\cite{topaz_denoiser}, the researchers trained the model with a larger dataset containing thousands of micrographs acquired in different conditions, which enables better generalizability of their neural networks. Besides Noise2Noise-based approaches, \gls*{nt2c}~\cite{nt2c} trains two neural denoisers and a \gls*{gan}-based noise synthesizer to decouple noise from micrographs that are generated using a \gls*{cryo-em} simulator based on real and synthetic data.

The implicit reconstruction approach by Kniesel~\textit{et~al.}~\cite{kniesel2022stem} also incorporates a normalizing flows neural network to model the noise distribution in the \gls*{stem} data in a signal-dependent manner, allowing them to separate the noise from the clean reconstruction.

Our approach is similar to Noise2Noise~\cite{noise2noise} in the way that our model is required to output similar noisy images (i.e., noisy \gls*{em} projections). However, in our simulator, we do not involve neural networks but a physical model.

\section{TEM Implementation}
To accelerate the \gls*{tem} simulator, we developed the \gls*{gputem} simulator that exploits parallelism of the GPU hardware using the CUDA API~\cite{cuda}. A large portion of the simulator pipeline is implemented in CUDA, while there are some  stages where we resort to the CPU implementation. We make a parallel version of the phantom generator and a per-slice parallelized version the electron-particle interaction algorithm.

For the phantom generator, the most time-consuming task is computing the scattering potential of biological structures (particles). The scattering potential of one particle is calculated by the sum of the scattering potential of all atoms that form the particle. In the original \gls*{tem}~\cite{rullgaard2011simulation} simulator, this task is performed sequentially.
To overcome this bottleneck in our \gls*{gputem}, we make a CUDA kernel for computing scattering potential for one atom, so the scattering potential of all atoms in the particle is computed in parallel at the atom level.

To calculate the electron-specimen interaction of each slice, we first determine which structures of the specimen are within the slice. Then we calculate the scattered electron wave after it interacts with these particles. This task is also performed sequentially in the original \gls*{tem}~\cite{rullgaard2011simulation} simulator. We make another CUDA kernel in our \gls*{gputem} simulator to compute the scattered electron wave after it interacts with one structure, so the electron-specimen interaction of each slice is computed in parallel at the particle level.

To implement our \gls*{difftem} simulator, we use Pytorch~\cite{paszke2019pytorch} and Taichi~\cite{hu2019taichi, hu2019difftaichi}. Pytorch is a well-known machine learning framework for parallel computing with GPU acceleration. With Pytorch, all gradients of each component in the simulator are calculated automatically and can be used for optimization with gradient descent. However, the \gls*{tem} simulator uses the multi-slice method for computing the electron-specimen interaction, which is a sequential algorithm and hard to be expressed in Pytorch operators. Therefore, to speed up and differentiate this computation, we use Taichi~\cite{hu2019taichi, hu2019difftaichi} and implement the calculation using Taichi kernels.

Our \gls*{difftem} simulator is differentiable, so all continuous parameters of different parts of the simulator can be optimized or learned from the real data via backpropagation and gradient descent. We demonstrate this capacity through two examples: (1) \gls*{mtf} parameters estimation and (2) denoising. From the pipeline, the final projection is obtained by adding Poisson noise to the noise-free projection and applying \gls*{mtf} afterward. Therefore, to estimate \gls*{mtf}, we perform backpropagation from real projections and use gradient descent for optimization. To denoise the real projections with gradient-based optimizations, we estimate the gradients of Poisson noise by using the method proposed by Schulman \textit{et~al.}~\cite{Schulman2015:GE-SCG} and apply gradient descent.

\subsection{Performance}
We demonstrated the scalability of our GPU-TEM in two aspects: (1) number of atoms and (2) number of virion instances. The performance of each simulator for a single instance each virus type is shown in~\autoref{table:performance_atoms}.

\begin{table}[H]
\small
\caption{Performance on w.r.t. number of atoms}
\label{table:performance_atoms}
\centering
\begin{tabular}{lccc}
\toprule
\textsc{Model} & \textsc{GPU-\gls*{tem}} & \textsc{\gls*{tem}~\cite{rullgaard2011simulation}}\\
\midrule
{\gls*{tmv}} & {4s} & {45s}\\
{ZIKA} & {1m46s} & {1m51s} \\
{SARS-CoV-2} & {} & {} \\
\bottomrule
\end{tabular}
\end{table}

From the~\autoref{table:performance_atoms}, the performance on single instance of ZIKA virus of our GPU-TEM is similar to the baseline~\cite{rullgaard2011simulation}, so we measured the performance on many instances of ZIKA, the results are shown in~\autoref{table:performance_instances}.
As in~\autoref{table:performance_instances}, our GPU-TEM performs better for scenes containing high number of instances.

\begin{table}[H]
\small
\caption{Performance on many instances of ZIKA}
\label{table:performance_instances}
\centering
\begin{tabular}{lccc}
\toprule
\textsc{No. of instances} & \textsc{GPU-\gls*{tem}} & \textsc{\gls*{tem}~\cite{rullgaard2011simulation}}\\
\midrule
{1} & {\textbf{1m 46s}} & {1m 51s} \\
{5} & {3m 1s} & {2m 9s} \\
{10} & {3m 10s} & {2m 16s} \\
{50} & {\textbf{3m 30s}} & {3m 34s} \\
{100} & {3m 48s} & {2m 16s} \\
{500} & {\textbf{5m 42s}} & {21m 17s} \\
{1000} & {\textbf{6m 22s}} & {35m 27s} \\
{5000} & {\textbf{18m 9s}} & {4h11m } \\
\bottomrule
\end{tabular}
\end{table}

\subsection{MTF parameters estimation}
\begin{figure*}[t]
\begin{center}
\includegraphics[width=\textwidth]{figures/Graphs.pdf}
\caption{\gls*{mtf} parameters estimation experiment results}
\label{fig:mtf_parameters_exp}
\end{center}
\end{figure*}

The effect of \gls*{mtf} is hard to see in the real \gls*{cryo-em} images because of low \gls*{snr} so we applied both ground truth \gls*{mtf} and predicted \gls*{mtf} on the noise-free image of synthetic scene and performed evaluation. The details of the synthetic data, the number of background images, and ground truth values of \gls*{mtf} parameters $(a, b, c, \alpha, \beta)$ are shown in~\autoref{table:desc_mtf_exp}.

\begin{table}[ht]
\small
\caption{Details of experiments for estimating \gls*{mtf} parameters}
\label{table:desc_mtf_exp}
\centering
\begin{tabular}{lccc}
\toprule
\textsc{Subject} & \textsc{\gls*{mtf} paramters} & \textsc{No. of images}\\
\midrule
{3 \gls*{tmv}} & {0.7, 0.2, 0.1, 10, 40} & {5, 10, 20, 50}\\
{$\approx 10$ SARS-CoV-2} & {0.6, 1.24, -0.84, 0.97, 0.97} & {48} \\
{Synthetic \gls*{mtf}} & {-2.86, 2.113, -0.1318, 57.288, 25.429} & {50} \\
\bottomrule
\end{tabular}
\end{table}

To evaluate the effect of the \gls*{mtf} on the noise-free image, we calculated the contribution of each frequency for the noise-free image and determined which frequency band yields the largest contribution by summing the magnitude for each band. The frequency ($\frac{1}{pixel}$) bands we used are $\left\{[0.0, 0.02), [0.02, 0.04), \ldots, [0.48, 0.5)\right\}$. If the predicted \gls*{mtf} aligns well with the most contributing frequency band of the ground truth \gls*{mtf}, their effects on the noise-free images are similar. We additionally calculated the mean and maximum difference of values for both images. If two \gls*{mtf} align with the ground truth \gls*{mtf} in the same degree, the better one is the \gls*{mtf} has the mean  and maximum differences values are smaller.

The correlation between the number of background images and the performance of each loss function on synthetic data is shown in \autoref{fig:mtf_parameters_exp}(a, b). \autoref{fig:mtf_parameters_exp}(a, b) shows that the \gls*{mse} performance decreases when the number of background images is greater than $10$, the gap between the predicted \gls*{mtf} and the ground truth \gls*{mtf} increases, the mean and maximum differences also increase. With $50$ background images, the images applied predicted \gls*{mtf} is blurring, and hard to see the details of each \gls*{tmv} virions. On the contrary, \gls*{ffl} performance improves when increasing the number of background images. The predicted \gls*{mtf} with $5$ background images is more aligned with the ground truth \gls*{mtf} than $10, 20$ background images. However, the gap between it the ground truth \gls*{mtf} is larger than the gap between the \gls*{mtf} estimated from $50$ background images. The image applied the predicted \gls*{mtf} with $50$ has the smallest mean and maximum differences values. The comparison of the best \gls*{mse} and \gls*{ffl} results is shown in~\autoref{fig:mtf_parameters_exp}(a, b). It is clear that \gls*{ffl} is better choice than \gls*{mse}.

We also investigated how the change of \gls*{ffl} $\alpha$ parameter affects the \gls*{mtf} estimation for different setting. The results are shown in  \autoref{fig:mtf_parameters_exp}(c, d). The best predicted \gls*{mtf} for \gls*{tmv}, SARS-CoV-2, and synthetic setting are obtained by using $\alpha = 2.0$, $\alpha = 0.5$, $\alpha = 2.0$, for \gls*{ffl} respectively. The best $\alpha$ value for each experiment is different because of the settings (electron dose, sample preparation, etc.) of different experiments are not the same, the distribution of frequencies in background images is changed. The results show that our system works effectively for different types of \gls*{mtf} curve. The first type is the real \gls*{mtf} (\gls*{tmv} and SARS-CoV-2) gives high modulation in the most contributing frequency band, so the details of virions are not blurry too much. The second type is the synthetic \gls*{mtf}, the modulation decreases in the most contributing frequency band, so the image is applied the \gls*{mtf} is more blurry.
Our predicted \gls*{mtf} for each case is aligned well with the ground truth in the most contributing frequency band, and give the similar effect on the noise free image.

\subsection{Denoising}
We compare the denoising capabilities of our approach in synthetic and real settings qualitatively and quantitatively.

\label{sec:denoising_synth_discussion}

\begin{figure}[H]
\begin{center}
\includegraphics[trim={3.5cm 0 3.5cm 0},clip=true,width=\linewidth]{figures/denoising_imgs/synthetic/comparison_ref_sample_num=2_dose=100_normalized_metrics.png}
\caption{Denoised result comparison of synthetic data setting (100, 2)}
\label{num_2_dose_100}
\end{center}
\end{figure}

Using a synthetic data with ground truth in the setting (100, 2)~, we have very limited amount of projections (i.~e., only 2) and a low dose which leads to a low \gls*{snr}. In settings like this where data are scarce and full of noise, a natural choice is to average the two available projections. Although the denoised result is not close to the noise-free reference, our method outperforms the simple average of the two available projections in terms of \gls*{ssim}~\cite{ssim}, \gls*{mse}, \gls*{psnr} and \gls*{ffl}. In the case of very high doses (i.~e., {700, 1000, 1500, 2000}), we show consistently better results than simple averaging in terms of \gls*{ffl}. For example, Fig.~\ref{num_50_dose_2000} shows that the \gls*{ffl} of our method is roughly 4 times better than simple averaging. The quantitative comparisons of these two figures show that our physics-based Noise2Noise-like model is effective.

\begin{figure}[H]
\begin{center}
\includegraphics[trim={3.5cm 0 3.5cm 0},clip=true,width=\linewidth]{figures/denoising_imgs/synthetic/comparison_ref_sample_num=50_dose=2000_normalized_metrics_unavg.png}
\caption{Denoised result comparison of synthetic data setting (2000, 50)}
\label{num_50_dose_2000}
\end{center}
\end{figure}

To compare with the original Noise2Noise~\cite{noise2noise} model and Topaz-Denoise~\cite{topaz_denoiser}, we present the our result of the setting (100, 50) denoted as Ours-50 in \autoref{synthetic_denoising_comparison}, in which the results of a trained Noise2Noise model, Topaz P and Topaz T are also shown. The Noise2Noise model cannot work even if we have a large amount of data available (i.~e., 5000). Therefore, we further tested the performance of Topaz models with our synthetic dataset. As we can see, the pretrained model (i.~e., \verb|Topaz P|) generalizes well in that it filters out high-frequency noise and preserves the particles, but it is worse than our method because the boundaries of particles are blurred. Our trained Topaz model (i.~e., \verb|Topaz T|) is better in that it filters out most of noise and recognizes the positions of the particles. However, it does not preserve the shape of the particles, which makes it unusable. On the contrary, our method can work in the setting with less data (i.~e., (100, 50)), which shows that our method is advantageous over training a Noise2Noise~\cite{noise2noise} model and a Topaz-Denoise~\cite{topaz_denoiser} model from scratch.

It is noteworthy that the Noise2Noise model and \verb|Topaz T| are trained with data settings (100, 5000) and (100, 4500) respectively. It is not realistic nor viable to obtain such large datasets in real-world scenarios, which means these two models are not usable. However, \verb|Topaz P| generalizes well in our setting, so we will use it to do data augmentation in the real-world scenario to generate more data for our denoising method.

For more results of our method under synthetic settings, we refer the reader to the \fl{our supplementary materials}.

\begin{figure}[H]
\begin{center}
\includegraphics[width=\linewidth]{figures/denoising_imgs/synthetic/Synthetic-comparison-3.png}
\caption{Comparison of denoising the synthetic data: images in the top row show original noisy projection (NP) and clean projection (GT), in bottom row are denosing results using Trained Noise2Noise result (Noise2Noise), pretrained Topaz Denoise (Topaz P), and Topaz Denoise trained on our images (Topaz T).}
\label{synthetic_denoising_comparison}
\end{center}
\end{figure}

\label{denoising_real_results}

In the setting of using the real data, we do not have noise-free ground truth projections to compare against to, so we first present results for qualitative comparisons.

\begin{figure}[H]
\begin{center}
\includegraphics[width=\linewidth]{figures/denoising_imgs/real_data_compare/RealDataDenoising.png}
\caption{Topaz: the real projection preprocessed by Topaz P; Noisy: the real noisy projection; Top: full views; Bottom: Zoom-in}
\label{denoising_with_real_data}
\end{center}
\end{figure}

\begin{table*}[t]
\centering\caption{SSIM, PSNR and FFL comparison} \label{SSIM_PSNR_FFL_comparison}
\resizebox{\textwidth}{!}{%
\begin{tabular}{cccccc|ccccc|ccccc}
\multicolumn{6}{c}{\textbf{SSIM}} & \multicolumn{5}{c}{\textbf{PSNR}} & \multicolumn{5}{c}{\textbf{FFL}}\\
\toprule
\multicolumn{1}{c}{}      & Noisy    & Topaz    & Ours (T)        & Ours (N)         & Ours (T+N)      & Noisy    & Topaz    & Ours (T)        & Ours (N)        & Ours (T+N)      & Noisy    & Topaz    & Ours (T)        & Ours (N)   & Ours (T+N)      \\
\midrule
\multicolumn{1}{r|}{Noisy} & 1.0000   & 0.6017   & 0.4139   & 0.7395   & \textbf{0.6547}   & 361.200  & 65.8250  & 65.8270  & 70.1500  & \textbf{66.1000}  & 0.0000   & 0.0096   & 0.0082   & 0.0032   & \textbf{0.0095}   \\
\multicolumn{1}{r|}{Topaz} & 0.6017   & 1.0000   & 0.6894   & 0.7184   & \textbf{0.8651}   & 65.8250  & 361.200  & 72.2000  & 68.9700  & \textbf{75.0400}  & 0.0096   & 0.0000   & 0.0012   & 0.0045   & \textbf{0.0008}   \\
\multicolumn{1}{r|}{Ours (T)}     & 0.4139   & 0.6894   & 1.0000   & 0.7797   & 0.8001   & 65.8270  & 72.2000  & 361.200  & 71.6200  & 76.2800  & 0.0082   & 0.0012   & 0.0000   & 0.0020   & 0.0001   \\
\multicolumn{1}{r|}{Ours (N)}     & 0.7395   & 0.7184   & 0.7797   & 1.0000   & 0.9098   & 70.1500  & 68.9700  & 71.6200  & 361.200  & 71.6700  & 0.0032   & 0.0045   & 0.0020   & 0.0000   & 0.0029   \\
\multicolumn{1}{r|}{Ours (T+N)}   & 0.6547   & 0.8651   & 0.8001   & 0.9098   & 1.0000   & 66.1000  & 75.0400  & 76.2800  & 71.6700  & 361.200  & 0.0095   & 0.0008   & 0.0001   & 0.0029   & 0.0000   \\
\bottomrule
\end{tabular}
}
\end{table*}
In \autoref{denoising_with_real_data}, we can see \verb|Topaz| has block-like artifacts and \verb|Noisy| contains more noise while the other techniques have few differences. However, if we zoom-in the images, as shown in the bottom of ~\autoref{denoising_with_real_data}, we can see (1) \verb|Topaz| has more blurring artifacts, (2) \verb|Noisy| has the most noise, and (3) \verb|T+N| is better than \verb|N| and \verb|T| in that it has less high-frequency noise while preserving clear boundaries. This makes our method \verb|T+N| the best among these five because of its good balance between denoising and preserving boundaries from the view of quality.

We also present quantitative results in Table~\ref{SSIM_PSNR_FFL_comparison}. Based on the bottom images from \autoref{denoising_with_real_data}, we can conclude that \verb|Topaz| contains the least noise while \verb|Noisy| contains the most noise. With this conclusion, the values in bold in the three tables show that (1) \verb|T+N| is the closest one to \verb|Topaz| besides \verb|Topaz| itself and (2) \verb|T+N| has a moderate amount of noise compared to \verb|T| and \verb|N|. This means our method is a good compromise between removing noise and preserving details, echoing with our qualitative comparison.

Besides the example in \autoref{denoising_with_real_data}, we give another example of denoising real data. From \autoref{denoising_with_real_data_zoomin_another}, we can draw the similar conclusions. \verb|Topaz| has the least noise, but compromised boundaries while \verb|T+N| has moderate amount of noise but also clearer boundaries and more visibility of spikes.

\begin{figure}[H]
\begin{center}
\includegraphics[width=\linewidth]{figures/denoising_imgs/real_data_compare/RealDataDenoising-2.png}
\caption{Denoising with real data zoom-in another example}
\label{denoising_with_real_data_zoomin_another}
\end{center}
\end{figure}

\fl{TODO: add results of reconstruction}

\subsection{MTF Parameters Estimation}

Because \gls*{mtf} parameters are independent of the specimen, we can optimize these parameters by using a few crops of background from the real projections. These crops are $M \times N$ pixels in size. They are considered as the labelled projections $I_{l_1}, I_{l_2}, \ldots, I_{l_n}$. The predicted projections $I_{p_1}, I_{p_2}, \ldots, I_{p_n}$ are created using our simulator with the same configuration for the sample, electron beam and optical system. The specimen of the predicted projection contains only the background. The predicted projections have the same size as the labelled projections. Both labelled projections and predicted projections are normalized to the range $[0, 1]$. Applying \gls*{mtf} is performed in the frequency domain which is motivation for us to choose a loss function that can capture the differences of the predicted and the labelled projection in both spatial and frequency domain. In this work, we choose \gls*{ffl}~\cite{jiang2021focal} combined with \gls*{mse} by summing them into the loss function.

The computational graph for estimating the \gls*{mtf} parameters is shown in \autoref{fig:compute_graph_of_detector_parameter_estimation_noisy}.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=\linewidth]{figures/detector_parameter_estimation_compute_graph_noisy.pdf}
\caption{Compute Graph of Detector Parameter Estimation: The \gls*{difftem} first receives noisy training data, computes the loss with simulated noisy data, applies the predicted \gls*{mtf}, and performs gradient descent for optimizing \gls*{mtf} parameters.}
\vspace{-1.5em}
\label{fig:compute_graph_of_detector_parameter_estimation_noisy}
\end{center}
\end{figure}

The optimization can be formally described as:
\begin{equation*}
\underset{\Theta_D}{\argmin}\ \mathbb{E}(L(\mtf_{\Theta_D}(\poisson(\tilde{I})),\  \nd)),
\end{equation*}
where $\Theta_D$ is the set of \gls*{mtf} parameters, $\tilde{I}$ a noise-free projection without \gls*{mtf}, $\nd$ is the noisy dataset, and $L$ is the loss function. $\nd$ can be the synthetic or real-world data. We conducted experiments for both.

The \gls*{ffl} loss is defined through a weight matrix and a frequency distance matrix.
$F_{l_1}, F_{l_2}, \ldots, F_{l_n}$ are 2D discrete Fourier transforms of labelled projections $I_{l_1}, I_{l_2}, \ldots, I_{l_n}$. $F_{p_1}, F_{p_2}, \ldots, F_{p_n}$ are 2D discrete Fourier transforms of predicted projections $I_{p_1}, I_{p_2}, \ldots, I_{p_n}$.

The \gls*{ffl} of labelled projection $I_{l_i}$ and predicted projection $I_{p_i}$ is:
\begin{equation}\label{ffl}
    \ffl(F_{l_i}, F_{p_i}) = \mean(w(F_{l_i}, F_{p_i}) * \left|F_{l_i} - F_{p_i}\right|^{2})
\end{equation}
where $\left|F_{l_i} - F_{p_i}\right|$ is frequency distance between the labelled projection $I_{l_i}$ and the predicted projection $I_{p_i}$, $w(F_{l_i}, F_{p_i}) = Normalize(\left|F_{l_i} - F_{p_i}\right|^{\alpha})$ is weighted matrix, $\alpha$ is the scaling factor for flexibility. The bigger $\alpha$ is, the more focus on optimization for hard frequencies (frequency components that are hard to synthesize). We conducted an experiment with different values for $\alpha$ and chose the best one.

The \gls*{mtf} parameters are optimized using the Adam optimizer~\cite{Kingma2014} $(\beta_{1} = 0.8, \beta_{2} = 0.992, lr = 0.1)$. We also experimented with optimization using only \gls*{mse} loss. The evaluation shows that combining \gls*{ffl} and \gls*{mse} performs best for the task.

\begin{algorithm}[tb]
\small
\begin{algorithmic}
\caption{\gls*{mtf} Parameters Optimization}
\label{mtf_optimization}
\State $I_{l_n} \gets Noisy\ Projections$
\State $I_{l_n} \gets \normalize(I_{l_n})$
\State $\Theta_D \gets random\ values$
\State $\tilde{I} \gets noise\verb|-|free\ projection\ without\ \mtf$
\State $\alpha \gets Learning\ Rate$
\State $IT \gets\ iterations$
\State $i \gets 1$
\While{$i\leq IT$}
    \State $I_{p_n} \gets \mtf(\poisson(\tilde{I}))$
    \State $I_{p_n} \gets \gets \normalize(I_{p_n})$
    \State $loss \gets \mse(I_{p_n}, I_{l_n}) + \ffl(I_{p_n}, I_{l_n})$
    \State $\Theta_D\gets \Theta_D - \alpha \frac{\partial loss}{\partial \Theta_D}$
    \State $i \gets i+1$
\EndWhile
\State \Return $\Theta_D$
\end{algorithmic}
\end{algorithm}

\subsection{Denoising}
\label{denoising_section}

As shown in \autoref{fig:electron_microscope_schematic_tem}, in the forward computation of our simulator, the input of quantization is a noise-free projection without \gls*{mtf} and the output of our simulator is a noisy projection.

Since our simulator is differentiable, we can do inverse computations via gradient descent, which means we can re-formulate the problem of denoising into an inverse problem. That is, given noisy data and a differentiable compute pipeline, to recover the noise-free projections without \gls*{mtf}, which are in the upstream of our forward computation. The optimization can be formulated as:
\begin{equation*}
\underset{\tilde{I}}{\argmin}\ \mathbb{E}(\slf(\mtf(Q(\tilde{I})),\  \rd)),
\end{equation*}
where $\tilde{I}$ is the recovered denoised projection without \gls*{mtf}, $\slf$ the surrogate loss function devised from the work of stochastic compute graph~\cite{Schulman2015:GE-SCG}, $\rd$ the real data (e.~g., noisy SARS-CoV-2 projections), and $L$ a loss function. To get a denoised projection $I_d$ from $\tilde{I}$, we bypass the quantization and apply \gls*{mtf} directly, which means $I_d=\mtf(\tilde{I})$.

The compute flow of the optimization of denoising is shown in \autoref{compute_graph_of_denoising}. The $\tilde{I}$ (i.~e., \verb|Projection without MTF|) is initialized with random values. The quantization models the shot noise of a detector. Given a Poisson distribution and its parameter $\lambda$, we can calculate the log likelihood of a sampled value $x\sim \poisson(\lambda)$, which is denoted as $\LL(x,\lambda)=ln(P(x|\poisson, \lambda))$ where $ln$ is the natural logarithm. After sampling with quantization and calculating log likelihoods, the results with shot noise $I_q$ (i.~e. \texttt{Noisy Projection without MTF} in \autoref{compute_graph_of_denoising}) are applied with \gls*{mtf}, which yields the predicted noisy projections $I_n$.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=\linewidth]{figures/denoising_compute_graph.pdf}
\caption{Compute graph of denoising: The projection without \gls*{mtf} is the optimization target. Log likelihood values are used in our surrogate loss function along with the comparison loss that computes the difference between noisy projection predictions and noisy projections in a dataset.}
\vspace{-1.5em}
\label{compute_graph_of_denoising}
\end{center}
\end{figure}

Inspired by the surrogate loss functions by Schulman \textit{et~al.}~\cite{Schulman2015:GE-SCG}, we take the surrogate loss function as Eq.~\ref{surrogate_loss_denoising} where $\logl$ equals $\LL(I_q, \tilde{I})$ and $\mean(\cdot)$ is the operator taking the mean value of all values in a tensor.

\begin{equation}\label{surrogate_loss_denoising}
    \slf(I_n, \rd, \logl) = \mean((\rd-I_n)^2*\logl).
\end{equation}

The whole optimization algorithm is shown in \autoref{denoising_optimization}. The normalization trick in Line 5 ensures consistent and standardized contrast and thus stabilizes our optimization. Note that due to the scarcity of real noisy projections, we perform data augmentations to generate meaningful training data. For one micrograph $P$ in a tilt series, we align its neighboring micrographs, which are denoted as $P_{-}$ and $P_{+}$, using perspective-based image alignment in OpenCV~\cite{opencv_library} because we assume that neighboring micrographs have minor differences compared to $P$, which is a reasonable assumption when the tilt angle increment is small, for example, 2-3 degrees. Therefore, in \autoref{denoising_optimization}, $P_n = \{P_{-}, P, P_{+}\}$. Furthermore, we preprocess $P_n$ with Topaz-Denoise~\cite{topaz_denoiser} to generate training data that are less noisy. The data augmentation is optional in the synthetic setting because we have sufficient amount of noisy training data, but in the real setting, it is necessary because (1) the available raw data is scarce (i.~e., only one), (2) the amount of information from only one projection with a low \gls*{snr} is limited, and (3) as discussed in Section~\ref{denoising_related_work}, our method is essentially a physical-model-based Noise2Noise~\cite{noise2noise}, which means it cannot be trained with only one projection. In the denoising optimization, we use Adam~\cite{Kingma2014} as the optimizer with $(\beta_1=0.8,\beta_2=0.992,lr=1.0)$.

\begin{algorithm}[tb]
\small
\begin{algorithmic}[1]
\caption{Denoising Optimization}
\label{denoising_optimization}
\State $P_n \gets Noisy\ Aligned\ Projections$
\State $\alpha \gets Learning\ Rate$
\State $IT \gets\ iterations$
\State $P_t \gets \topazOp(P_n)$
\State $P \gets \normalize(P_n\cup P_t)$
\State $\tilde{I} \gets random\ values$
\State $i \gets 1$
\While{$i\leq IT$}
    \State $I_q \gets Q(\tilde{I})$
    \State $\logl\gets \LL(I_q, \tilde{I})$
    \State $I_n \gets \mtf(I_q)$
    \State $\loss \gets \slf(I_n, P, ll)$
    \State $\tilde{I}\gets \tilde{I} - \alpha \frac{\partial \loss}{\partial \tilde{I}}$
    \State $i \gets i+1$
\EndWhile
\State $I_d=\mtf(\tilde{I})$
\end{algorithmic}
\end{algorithm}
\subsection{Performance evaluation}
We compared the performance of the original TEM simulator and our GPU-TEM. In the phantom generator we used the model of SARS-CoV-2 virus modelled by Nguyen~\textit{et~al.}~\cite{nguyen2021modeling}, ZIKA virus \gls*{pdb}id 5IRE and \gls*{tmv} PDBid 2OM3. The details of each model are shown in~\autoref{table:desc_virus_model}. All performance experiments were done on a Windows system with one processor Intel(R) Xeon(R) CPU E5-2687W v4 3GHz and one GPU NVIDIA GeForece RTX 3090.

\begin{table}[H]
\small
\caption{Details of the models used in the performance evaluation.}
\label{table:desc_virus_model}
\centering
\begin{tabular}{lccc}
\toprule
\textsc{Virions}& \textsc{No. of pdbs} & \textsc{No. of atoms}  \\
\midrule
{\gls*{tmv}} & {1} & {12505} \\
{ZIKA} & {1} & {21149} \\
{SARS-CoV-2} & {17} & {$\approx 13$ millions} \\
\bottomrule
\end{tabular}
\end{table}

\subsection{MTF parameters estimation}
For estimating \gls*{mtf} parameters, we performed experiments on both synthetic and real data. The synthetic data is created from the envelope of SARS-CoV-2 virion and \gls*{tmv} virion. The real dataset is provided by Yao~\textit{et~al.}~\cite{Sai2020:CryoET-Pipeline} consists of 20 tilt series containing several SARS-CoV-2 virions each, but we only use cropped sub-image containing background for optimization. The real data was acquired using a Titan Krios microscope with K3 detector from Gatan, operated at a voltage of $300$ kV.
The ground truth of \gls*{mtf} parameters for real data was obtained by performing curve-fitting from \gls*{mtf} curve of Gatan K3 detector\footnote{https://www.gatan.com/K3\#resources}. The ground truth values of \gls*{mtf} parameters for synthetic data \gls*{tmv} virion is provided by Rullg{\aa}rd~\textit{et~al.}~\cite{rullgaard2011simulation} in their simulation and calibrated with their real data which is the micrograph of \gls*{tmv} virions measured by a Philips CM 200 FEG \gls*{tem} microscope. We conducted experiments for each loss function (\gls*{mse}, \gls*{mse} + \gls*{ffl}) with different number of background images of synthetic data to see the correlation between performance of each loss and number of images. Then we conducted the experiments on the real data with the number of images similar to the number that gives the best performance on synthetic data. We also conducted experiments with different \gls*{ffl} loss $\alpha$ values:  $(0.0, 0.5, 1.0, 2.0)$. Beside, to demonstrate the generalization of the estimation method, we also experimented with synthetic \gls*{mtf} that has the shape is different with real \gls*{mtf}.
All estimating \gls*{mtf} parameters experiments were conducted on Linux system with four GPUs NVIDIA RTX A6000.

\subsection{Denoising}
First, we conducted experiments using a synthetic dataset generated with our simulator given an atomistic model (i.~e., \gls*{tmv} model from the research by Sachse \textit{et~al.}~\cite{SACHSE2007812}). In this setting, we obtain a clean ground truth data that are synthesized by bypassing the quantization.

In simulation and in real electron microscopy process, the dose, measured in electrons per nm\textsuperscript{2} per projection, is a major factor that determines \gls*{snr}. We tested our denoising method with synthetic datasets with applying varying doses. The higher the dose is, the higher \gls*{snr} is. For generating a tilt series of 60 projections, 100 electrons/nm\textsuperscript{2} per projection is a typical low dose while 200 is a typical high dose. We tested datasets of doses \{100, 200, 300, 500, 700, 1000, 1500, 2000\} because it is possible to use a very high dose (e.~g., 1000) if the number of total projections is low and thus the total dose stays the same. Moreover, the number of available projections of the same scene (i.~e., same tilt angle and specimen) also affects denoising quality. We tested our method on the numbers of available projections of \{2, 5, 10, 25, 50\}. We denote varying data settings with pairs, i.e., dose, the number of available projections.

To compare with the original Noise2Noise~\cite{noise2noise} model and Topaz-Denoise~\cite{topaz_denoiser}, we trained a Noise2Noise model and a Topaz-Denoise~model with our synthetic data. The data setting for training Noise2Noise model is (100, 5000) while it is (100, 4500) for training Topaz-Denoise~model. We also use pretrained Topaz-Denoise~model to contrast the results.

In the experiments using real micrographs, we use the tilt series from Yao \textit{et~al.}~\cite{Sai2020:CryoET-Pipeline}, which is SARS-CoV-2. The dose in their data collection is 320 electrons/nm\textsuperscript{2} per projection. Therefore, without data augmentations, the setting is (320, 1). As discussed in Section~\ref{denoising_section}, our method cannot work in a scenario with only single projection, so in our experiments, we apply a number of data augmentations. In the setting \texttt{Ours(T)}, the data augmentation with Topaz P is used, yielding one more projection, which means the setting is (320, 2). In the setting denoted as \texttt{Ours(N)}, neighboring projections of a projection are aligned and used, so the data setting is (320, 3). In the setting \texttt{Ours(T+N)}, Topaz P and aligned neighboring projections are used, so the data setting is (320, 6). The case of \texttt{Ours(T+N)} matches exactly our \autoref{denoising_optimization} while in the case of \texttt{Ours(T)}, $P_n=\{P\}$, and in the case of \texttt{Ours(N)}, preprocessing with Topaz-Denoise~is not used (i.~e., $P_{train}=Normalize(P_n)$). The results of the three cases (i.~e., \texttt{Ours(T)}, \texttt{Ours(N)} and \texttt{Ours(T+N)}) will be presented along with the projection preprocessed by \texttt{Topaz P} and the noisy projection in Section~\ref{denoising_real_results}.

\subsection{Contribution}
This paper presents several technical contributions that and visualization applications:
\begin{itemize}[noitemsep,topsep=0pt]
    \item The overall architecture of the first differentiable microscopy simulator that is scalable to render realistic microscopy simulations within seconds or minutes.
    \item Learning detector parameters from real projections, to generate synthetic micrographs with the same detector characteristics that lead to similar visual appearance of the synthetic and the real projection.
    \item Utilization of deterministic and stochastic differentiation for various stages of the microscopy pipeline.
    \item Generation of large micrographs \emph{on par} with physical dimensions and complexity of real-world microscopy.
    \item A new denoising method which learns from noisy projections to create a noise-free version of real-world projections. Subsequent tomographic reconstruction then results in a noise-free volume that can be used for 3D visualization.
\end{itemize}
\section{Introduction}
\input{includes/introduction}

\section{Related work}
\input{includes/related-work}

\section{Architecture overview}
\input{includes/archiecture-overview}

\section{Learning components}
\input{includes/learning-components}

\section{Reverse-Engineering MTF}
\input{includes/MTF}

\section{Scalable TEM Simulation}
\input{includes/GPUTEM}

\section{Denoising Micrographs}
\input{includes/denoising}

\section{Conclusion}
\input{includes/conclusion}

\section*{Acknowledgment}

The authors would like to thank Sai Li and his team at School of Life Sciences, Tsinghua University, China for sharing the SARS-CoV-2 \gls*{cryo-em} data for this work.

\clearpage
\bibliographystyle{unsrt}


\section{Introduction}
Atomizing sprays have a wide range of applications in industrial and environmental fields (e.g. fuel injection~\cite{Zhao1999}, agricultural sprays~\cite{Law2001}, pharmaceutical sprays~\cite{Maa1999}, and air-sea interaction~\cite{Fuentes2010}). Consequently, the body of research on this topic is vast and spans multiple centuries~\cite{Savart1833,Rayleigh1878}. Researchers have made significant progress in understanding many aspects of the atomization process, however, the mechanisms that drive instabilities in coherent liquid structures and ultimately cause breakup are still widely unresolved. This is, in part, due to limitations in experimental data collection from atomizing flows. Visual analysis of atomizing jets is difficult because of the scale and speed at which these systems develop. Furthermore, the development of droplets creates an opaque cloud, which blocks view of the liquid core and severely limits the ability to study primary instability and breakup mechanics. Recent advances in ballistic~\cite{Linne2005,Sedarsky2010,Slangen2016} and x-ray~\cite{Wang2006,Heindel2018,Li2021} imaging of atomizing flows have attempted to remedy visualization issues. These methods provide snapshots of the atomization process and lack important temporal information. Because of challenges with current experimental methodology, numerical simulations, and particularly high-fidelity simulations that resolve the relevant time and length scales, have been developed and provide an alternative method to study the physics of atomization.

Advances in numerical methods and computational efficiency within the past decade have greatly expanded the capabilities of numerical simulations. Researchers are now able to simulate multiphase flows with exceptional accuracy. Numerically simulating multiple phases is not a trivial process, however. These systems involve a wide range of topological scales in fluid and turbulent structures in addition to jumps in physical qualities at the phase interfaces. Significant effort has been and is still invested in improving numerical simulation methodology to more accurately and efficiently simulate atomizing systems. Accordingly, most high-fidelity simulation studies on atomization to date are focused on the development of numerical methods and testing against experimental results rather than using the simulations to produce novel atomization data. But, with the pending maturity of the field, several research groups have been able to gather some insightful statistics from high-fidelity atomization simulations.

Following Gorokhovski and Herrmann's review of atomization simulations~\cite{Herrmann2008}, the field developed rapidly. Among the first high-fidelity simulations of atomizing liquid were Shinjo and Umemura~\cite{Shinjo2010} and Desjardins and Pitsch~\cite{Desjardins2010} who simulated round and planar liquid jets, respectively. These works provided a foundation for future simulation studies by identifying numerous breakup mechanisms and providing qualitative descriptions of previously unclassified processes. Following these pioneering studies, various other works were published which utilized high-fidelity simulations to expand upon existing descriptions and identify novel breakup mechanisms. Several of these groups incorporated vorticity measurements into their post-processing to strengthen descriptions of breakup by integrating this additional dimension to their analysis \cite{Shinjo2011b,Jarrahbashi2014,Watanabe2014,Jarrahbashi2016,Li2016,Zandian2017,Zandian2018}. Additionally, many researchers focus their efforts on comparing linear-stability theory~\cite{Reitz1987,Beale1999} and high-fidelity simulations~\cite{Shinjo2011b,Fuster2013,Deshpande2015,Trujillo2018,Krolick2018}. These studies process visualizations of the atomizing jet to analyze and track wavelengths as they become more unstable, ultimately leading to breakup. Several recent studies have obtained turbulence statistics from throughout the simulation using constrained spatial and/or temporal domains to limit data-processing requirements~\cite{Hasslberger2019,torregrosa2020,Salvador2018}.

The numerical studies mentioned above have greatly improved understanding of the atomization process. They developed robust qualitative descriptions of breakup events and utilized bulk quantities to elucidate some of the underlying physics of liquid breakup. However, despite the recent advances in atomization simulations, there remain significant gaps in our understanding of the process. This is primarily due to a lack of available methodology for obtaining many relevant statistics from high-fidelity simulations. In particular, current research methods have not been able to extract local, temporally continuous, quantitative data on liquid breakup. The data include information on droplet shape and size characteristics, local flow field data in both the liquid and gas phases, and information on how these values evolve temporally and spatially throughout a spray system. An appreciable challenge obtaining this data is derived from the massive size of resultant data sets; these can be hundreds of terabytes or even larger. Parsing such data sets for relevant information is not practical and is often impossible.

While numerical simulations are becoming more efficient, it remains computationally expensive and requires the use of high-performance computing. Reduced-order atomization models aim to provide a viable alternative to high-fidelity simulations at a fraction of the computational cost. Some prominent examples of atomization models include the Taylor analogy breakup (TAB) model \cite{TAB1987}, the Pilch-Erdman (PE) model \cite{Pilch1987}, the Kelvin-Helmholtz-Rayleigh-Taylor (KH-RT) model \cite{Reitz1987,Teclllwlvgy1987,Patterson1998}, the Eulerian-Lagrangian spray and atomization (ELSA) model \cite{Wang2011}, the bag type breakup (BTB) model \cite{Wang2014}, the multi-mode breakup (MMB) model \cite{Wang2015}, and the modified Taylor analogy breakup (MTAB) model \cite{Sula2020}. These atomization models try to predict how and when a liquid structure breaks apart without fully resolving the physics. As mentioned above, however, there is little information on the mechanisms that control breakup events and even less information on the process of atomization. Model developers would greatly benefit from local data sampled from breakup events throughout the simulation.

The present work aims to address the above obstacles through the improvement of a framework, first introduced in Rubel and Owkes~\cite{RubelOwkes2019}, which extracts relevant data from breakup and coalescence events in atomization simulations. Through the use of this tool, many researchers would gain access to extensive and easily accessible data sets from atomization simulations. Using this data, they could build on previous work, introduce high-quality statistical analyses, and significantly advance the field. Additionally, the structure geometry and flow field can be assessed to quantify the details of the breakup conditions. This new information will advance our understanding of the atomization process and provide statistics with which to inform and improve reduced-order atomization models.

\section {Methods}
The droplet physics extraction tool, originally proposed in Rubel and Owkes \cite{RubelOwkes2019}, is improved and applied in this work. The tool is implemented into high-fidelity atomization simulations to gather data from discrete events within these simulations coinciding to liquid breakup and coalescence. This section will outline the methodology for this process in detail.

\subsection{Computational Platform}
The proposed extraction tool can be applied to any high-fidelity Navier-Stokes solver, given the solver incorporates two identification numbers used to determine breakup and coalescence events, which will be defined in the following sections. This work employs the NGA computational platform \cite{Desjardins2008a,Desjardins2008b,Owkes2014,Owkes2017}. NGA solves the two-phase formulations of the Navier-Stokes mass and momentum conservation equations, defined as
\begin{equation}
    \frac{\partial \rho_{\phi}}{\partial t}+\nabla\cdot(\rho_{\phi}\bm{u}_{\phi})=0
\end{equation}
and
\begin{equation}
    \frac{\partial \rho_{\phi} \bm{u}_{\phi}}{\partial t}+\nabla\cdot(\rho_{\phi}\bm{u}_{\phi}\otimes \bm{u}_{\phi})= -\nabla p_{\phi} +\nabla\cdot(\mu_{\phi}[\nabla \bm{u}_{\phi} +\nabla \bm{u}^{T}_{\phi}]) +\rho_{\phi}\bm{g}
\end{equation}
respectively. In these equations, $\rho_{\phi}$ is density, $\bm{u}_{\phi}$ is the velocity vector, $p_{\phi}$ is the pressure, $\mu_{\phi}$ is the dynamic viscosity, $t$ is time, $\bm{g}$ is the gravitational acceleration vector and $\phi$ is the phase indicator with either $\phi = g$ (gas) or $\phi = l$ (liquid). These equations are solved on a staggered Cartesian mesh, in which scalar values, such as pressure, are stored at cell centers and velocities at the cell faces. Time is discretized with an iterative Crank-Nicolson formulation and a semi-implicit correction is applied on each sub-iteration \cite{Choi1994}. Away from the phase interface, mass, momentum, and any other scalars are transported with conservative, high-order finite difference operators~\cite{Desjardins2008a}. Near the interface, an un-split semi-Lagrangian geometric volume-of-fluid (VOF) method is utilized to ensure conservative transport of mass and momentum \cite{Owkes2014,Owkes2017}. The interface is reconstructed using a piece-wise linear interface reconstruction (PLIC) \cite{Youngs1982} with interface normal vectors computed  using the efficient least-squares VOF interface reconstruction algorithm (ELVIRA) \cite{Pilliod2004}. The pressure Poisson equation is solved utilizing the ghost fluid method \cite{Fedkiw1999} and a black box solver \cite{Dendy1982}. Interface curvature is computed with the adjustable curvature evaluation scale (ACES) method \cite{Owkes2018}. NGA is fully parallelized with message passing interface (MPI) and scales well to tens-of-thousands of cores~\cite{Desjardins2010}.

\subsection{Breakup and Coalescence Event Identification}
The coalescence and split identification processes operate through the implementation of two identification numbers, which are integers unique to every independent liquid structure within a simulation. These values are referred to as the structure identification number $\mathcal{S}$ and the liquid identification number $\mathcal{L}$. Every liquid structure in the simulation domain is tagged with each of these values. The process in which $\mathcal{S}$ and $\mathcal{L}$ are assigned, transported, and reassigned is the basis for identifying liquid splitting and merging events within a simulation.

\subsubsection{Identification Numbers} \label{id_num}
Liquid identification numbers ($\mathcal{L}$) are transported with the liquid and provide a time history of where liquid moves.  Therefore, the $\mathcal{L}$ values are not reassigned at every timestep. They are persistent through time and are only reassigned after a split event creates a new liquid structure or a coalescence event destroys a liquid structure. $\mathcal{L}$ values are transported with the liquid they occupy, so $\mathcal{L}$ transport must be tied into the simulation's liquid transport scheme. In our implementation, $\mathcal{L}$ fluxes are calculated in the same manner as the semi-Lagrangian liquid fluxes in NGA's geometric VOF scheme~\cite{Owkes2014,Owkes2017}. A list ($\mathcal{L}_\mathrm{Flux}$) of all the $\mathcal{L}$ that flux into a cell during a timestep is constructed. This is done by first constructing lists at each cell face with signed $\mathcal{L}$ values indicating the direction. After lists associated with each cell face are computed, they are compiled, to create a list of all $\mathcal{L}$'s within a cell. These lists are looped through in a later section of the code to identify coalescence. Fig.~\ref{fig:LID_flux} displays a visual representation of this process. See \ref{appendixa} for the $\mathcal{L}$ transport and coalescence identification algorithm.

Structure identification numbers ($\mathcal{S}$) are assigned at every timestep following liquid transport. The values of $\mathcal{S}$ are not consequential, provided that every structure is assigned a unique value. Independent structures are identified and tagged using a band-growth algorithm first described by Herrmann~\cite{Herrmann2010}. In the present application, the algorithm operates as follows. To begin, $\mathcal{S}$ is assigned on the entire domain to be zero. Then, the domain is looped through to find the first untagged liquid node, meaning the liquid volume fraction is greater than zero and $\mathcal{S}$ = 0. This point is then tagged with an $\mathcal{S}$ and the cells around that node are looped through until the entire continuous structure is identified and assigned the $\mathcal{S}$ value. A logical statement determines whether an adjacent cell will be assigned the same $\mathcal{S}$ value as the present cell. The statement requires that: 1) there must be liquid in that cell, 2) it must be untagged ($\mathcal{S} = 0$), and  3) the $\mathcal{L}$ of the adjacent cell must equal that of the current cell. The third statement ensures that liquid structures do not re-merge immediately following breakup as structures advance and naturally inhabit adjacent cells. See Section~\ref{fictitious_e} for a more detailed description of these criteria. Herrmann's work also described the parallelization of this process, which is utilized in the extraction tool. This requires special consideration be given to communicating $\mathcal{S}$ between processors. Many other methods exist to assign the structure identification number. This type of tagging problem is known as Connected Component Labeling, e.g.~\cite{Hendrickson2020,Hoshen1976}.

\begin{figure}[htbp]
 \centering
 \begin{tikzpicture}
    \draw[step=4cm,thin] (-1.9,-1.9) grid (5.9,5.9);
    \draw (0,0) -- (0,4);
    \fill [fill={myblue}] (-1.2,2) circle (1.2);
    \fill [fill={myellow}] (1.4,-1.4) circle (1.2);
    \node [text width=1.5cm,align=center] at (-1.2,2) {\large$\mathcal{L}=1$};
    \node [text width=1.5cm,align=center] at (1.4,-1.4) {\large$\mathcal{L}=2$};
    \filldraw[gray, opacity = 0.2] (0,4) -- (-1,3) -- (-1,-1) -- (0,0);
    \draw [help lines] (0,4) -- (-1,3) -- (-1,-1) -- (0,0);
    \draw [help lines] (-1,-1) -- (3,-1) -- (4,0);
    \filldraw[gray, opacity=0.2] (-1,-1) -- (3,-1) -- (4,0) -- (0,0);
    \node [text width=5cm,align=center] at (2,2) {\Large $\displaystyle \mathcal{L}_\mathrm{Flux}\!\!\!=\!\!\![1,2]$};
    \draw[->] (4,0) -- (5,1);
    \draw[->] (0,0) -- (1,1);
    \draw[->] (0,4) -- (1,5);
  \end{tikzpicture}
 \caption{An example of two droplets fluxing into the same cell. Gray regions are semi-lagrangian flux regions~\cite{Owkes2014}, which carry $\mathcal{L}$ into the center cell due to the velocity shown with vectors. The $\mathcal{L}$'s of the two droplets are stored in the list $\mathcal{L}_\mathrm{Flux}$.}
 \label{fig:LID_flux}
\end{figure}

\subsubsection{Breakup Identification}
Liquid structure splitting or breakup is identified when two different structures (two unique $\mathcal{S}$'s) have the same $\mathcal{L}$. This $\mathcal{L}$ is persistent through time, indicating that these two droplets were part of the same structure at the previous timestep. The logical statement, described in Section~\ref{id_num}, loops over adjacent cells to connect liquid structures. Two structures become independent of each other when they are not in adjacent cells i.e., there is a full cell separating them. This is when the $\mathcal{S}$ assignment portion of the code assigns two structures with unique $\mathcal{S}$'s. The two structures have the same $\mathcal{L}$, indicating that a breakup occurred. When a breakup is identified, data are extracted from the event, and a new $\mathcal{L}$ is assigned to the smaller structure(s) created by the breakup event. See $t=0$ and $t=1$ in Fig~\ref{fig:splitmerge} for a visual representation of this process.

\subsubsection{Coalescence Identification}\label{coalescence}
Coalescence of liquid structures is identified using the $\mathcal{L}_{flux}$ lists mentioned in Section~\ref{id_num}. More than one unique $\mathcal{L}$ in a cell indicates that two structures are merging together. When coalescence is identified, data are extracted from the event, and the new merged structure takes the smaller $\mathcal{L}$ of the two parent structures.

\begin{figure}[htbp]
 \centering
 \includegraphics[width=1\textwidth]{merge_split_process.PNG}
 \caption{An example of two droplets undergoing breakup and coalescence. Two droplets start on the left,  the top droplet breaks into two droplets.  The two droplets are identified and the split event is processed.  The bottom two droplets coalesce, the merge is identified and the merge event is processed. Note that the $\mathcal{S}$ can vary between timesteps, but the $\mathcal{L}$ is persistent until a split or merge causes the value to change.}
 \label{fig:splitmerge}
\end{figure}
\subsubsection{Addressing Fictitious Events} \label{fictitious_e}
An issue from the previous work, described in Rubel and Owkes~\cite{RubelOwkes2019} was the occurrence of numerous fictitious merge and split events. These occurred when the tool identified a structure that broke into two structures and then the two structures coalesced. The algorithm identified this process and repeated many times, when in fact the liquid structures never re-merged following breakup.
The error stemmed from the original $\mathcal{S}$ assignment algorithm. Droplets in adjacent cells were automatically treated as the same structure, which resulted in recently split liquid structures being incorporated back into their parent droplet when they were in neighboring cells. To avoid an accumulation of non-physical results, the authors opted to run the split identification portion of the tool only every 10 timesteps. The present work changed the parameters that define a coalescence event to prevent droplets from merging together immediately after breakup. This allows the tool to run every timestep, bolstering it's accuracy and usefulness. See Fig.~\ref{fig:updates} for a visualization of this issue and the result following the updates in this work.
\begin{figure}
  \centering
    \begin{tikzpicture}

      \node [text width=1.5cm,align=center] at (3.8,0) {$t = 3$};
      \node [text width=1.5cm,align=center] at (3.8,2.4) {$t = 2$};
      \node [text width=1.5cm,align=center] at (3.8,4.8) {$t = 1$};
      \node [text width=1.5cm,align=center] at (3.8,7.2) {$t = 0$};
      \node [text width=4cm,align=center] at (0.2,10.4) {Original Algorithm~\cite{RubelOwkes2019}};
      \node [text width=4cm,align=center] at (7.4,10.4) {Updated Algorithm};

      \draw (-2.6,8.4) -- (3,8.4) -- (3,10.2) -- (-2.6,10.2) -- cycle;
      \node at (0.2,9.6) {VOF $> 0$};
      \node at (0.2,9) {$\mathcal{S} = 0$};
      \draw (4.6,8.4) -- (10.2,8.4) -- (10.2,10.2) -- (4.6,10.2) -- cycle;
      \node at (7.4,9.8) {VOF $> 0$};
      \node at (7.4,9.3) {$\mathcal{S} = 0$};
      \node at (7.4,8.8) {$\mathcal{L}_i$ = $\mathcal{L}_{i+1}$};

      \draw [step=0.8cm,gray,thin] (-2.6,-1.1) grid (3,8.3);
      \fill [{myblue}] (-2.6,-0.6) rectangle (0,0.6);
      \filldraw [{myblue}] (0,-0.6) -- (0.6,0) -- (0,0.6) -- cycle;
      \node [text width=1.5cm,align=center] at (-1.8,0.3) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (-1.8,-0.3) {$\mathcal{L}=1$};
      \filldraw [{myellow}] (1.8,0) -- (2.4,-0.6) -- (3,0) -- (2.4,0.6) -- cycle;
      \node [text width=1.5cm,align=center] at (2.4,0.2) {$\mathcal{S}=2$};
      \node [text width=1.5cm,align=center] at (2.4,-0.2) {$\mathcal{L}=2$};

      Row
      \fill [{myblue}] (-2.6,1.8) rectangle (0,3);
      \filldraw [{myblue}] (0,1.8) -- (0.3,2.4) -- (0,3) -- cycle;
      \node [text width=1.5cm,align=center] at (-1.8,2.7) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (-1.8,2.1) {$\mathcal{L}=1$};
      \filldraw [{myblue}] (1.4,2.4) -- (1.6,2) -- (2.4,2) -- (2.6,2.4) -- (2.4,2.8) -- (1.6,2.8) -- cycle;
      \node [text width=1.5cm,align=center] at (2,2.6) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (2,2.2) {$\mathcal{L}=1$};

      \fill [{myblue}] (-2.6,4.2) rectangle (-0.8,5.4);
      \filldraw [{myblue}] (-0.8,4.2) -- (-0.2,4.8) -- (-0.8,5.4) -- cycle;
      \node [text width=1.5cm,align=center] at (-1.8,5.1) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (-1.8,4.5) {$\mathcal{L}=1$};
      \filldraw [{myellow}] (1,4.8) -- (1.6,4.2) -- (2.2,4.8) -- (1.6,5.4) -- cycle;
      \node [text width=1.5cm,align=center] at (1.6,5.0) {$\mathcal{S}=2$};
      \node [text width=1.5cm,align=center] at (1.6,4.6) {$\mathcal{L}=2$};

      \fill [{myblue}] (-2.6,6.6) rectangle (-0.8,7.8);
      \filldraw [{myblue}] (-0.8,6.6) -- (-0.5,7.2) -- (-0.8,7.8) -- cycle;
      \node [text width=1.5cm,align=center] at (-1.8,7.5) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (-1.8,6.9) {$\mathcal{L}=1$};
      \filldraw [{myblue}] (0.6,7.2) -- (0.8,6.8) -- (1.6,6.8) -- (1.8,7.2) -- (1.6,7.6) -- (0.8,7.6) -- cycle;
      \node [text width=1.5cm,align=center] at (1.2,7.4) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (1.2,7.0) {$\mathcal{L}=1$};

      \draw [step=0.8cm,gray,thin] (4.6,-1.1) grid (10.2,8.3);
      \fill [{myblue}] (4.6,-0.6) rectangle (7.2,0.6);
      \filldraw [{myblue}] (7.2,-0.6) -- (7.8,0) -- (7.2,0.6) -- cycle;
      \node [text width=1.5cm,align=center] at (5.4,0.2) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (5.4,-0.2) {$\mathcal{L}=1$};
      \filldraw [{myellow}] (9,0) -- (9.6,-0.6) -- (10.2,0) -- (9.6,0.6) -- cycle;
      \node [text width=1.5cm,align=center] at (9.6,0.2) {$\mathcal{S}=2$};
      \node [text width=1.5cm,align=center] at (9.6,-0.2) {$\mathcal{L}=2$};

      \fill [{myblue}] (4.6,1.8) rectangle (7.2,3);
      \filldraw [{myblue}] (7.2,1.8) -- (7.5,2.4) -- (7.2,3) -- cycle;
      \node [text width=1.5cm,align=center] at (5.4,2.6) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (5.4,2.2) {$\mathcal{L}=1$};
      \filldraw [{myellow}] (8.6,2.4) -- (8.8,2) -- (9.6,2) -- (9.8,2.4) -- (9.6,2.8) -- (8.8,2.8) -- cycle;
      \node [text width=1.5cm,align=center] at (9.2,2.6) {$\mathcal{S}=2$};
      \node [text width=1.5cm,align=center] at (9.2,2.2) {$\mathcal{L}=2$};

      \fill [{myblue}] (4.6,4.2) rectangle (6.4,5.4);
      \filldraw [{myblue}] (6.4,4.2) -- (7,4.8) -- (6.4,5.4) -- cycle;
      \node [text width=1.5cm,align=center] at (5.4,5.0) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (5.4,4.6) {$\mathcal{L}=1$};
      \filldraw [{myellow}] (8.2,4.8) -- (8.8,4.2) -- (9.4,4.8) -- (8.8,5.4) -- cycle;
      \node [text width=1.5cm,align=center] at (8.8,5.0) {$\mathcal{S}=2$};
      \node [text width=1.5cm,align=center] at (8.8,4.6) {$\mathcal{L}=2$};

      \fill [{myblue}] (4.6,6.6) rectangle (6.4,7.8);
      \filldraw [{myblue}] (6.4,6.6) -- (6.7,7.2) -- (6.4,7.8) -- cycle;
      \node [text width=1.5cm,align=center] at (5.4,7.4) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (5.4,7.0) {$\mathcal{L}=1$};
      \filldraw [{myblue}] (7.8,7.2) -- (8,6.8) -- (8.8, 6.8) -- (9,7.2) -- (8.8,7.6) -- (8,7.6) -- cycle;
      \node [text width=1.5cm,align=center] at (8.4,7.4) {$\mathcal{S}=1$};
      \node [text width=1.5cm,align=center] at (8.4,7.0) {$\mathcal{L}=1$};
    \end{tikzpicture}
    \caption{The diagram on the left shows the issue with fictitious merging and splitting events in the original algorithm~\cite{RubelOwkes2019}. The diagram to the right shows how the split and merge identifier works in the present work. Notice that at $t = 2$ in the updated code the structures do not merge back together. This is because the criterion for coalescence was changed to require the structures to exist within the same cell, not adjacent cells.}
    \label{fig:updates}
\end{figure}

As stated above, in Section~\ref{coalescence}, coalescence now occurs when two structures flux into the same cell. This methodology requires that special care be given to the $\mathcal{S}$ assignment algorithm to ensure that separate structures in adjacent cells are not combined. To address this, the $\mathcal{L}$ of droplets was incorporated into the band-growth algorithm. Fig.~\ref{fig:updates} illustrates this updated $\mathcal{S}$ assignment process. At $t = 0$ a droplet is breaking up so that there is liquid in adjacent cells. Since a split has not been identified yet, the $\mathcal{L}$ of both cells is the same, resulting in a continuous structure. Then, at $t = 1$ the structure on the right moves, leaving a full cell between it and the other structure. The cells adjacent to the blue structure are void of liquid, so the $\mathcal{S}$ value is only assigned to that structure. Since the droplet on the right becomes independent, it is assigned a new $\mathcal{S}$, a split is processed and it is also assigned a new $\mathcal{L}$. Finally, at $t = 2$, previously, the two droplets would have merged back together as they occupy neighboring cells.  But the structures remain separate with the new algorithm, because the $\mathcal{L}$ of the adjacent droplet does not equal the $\mathcal{L}$ of the blue droplet, they remain separate structures and will not coalesce unless they enter the same cell.
\label{sec:fict_events}

\subsection{Data Extraction}
Following identification of breakup or coalescence in the simulation, data are extracted from that location and time. This is intended to provide access to the local conditions of these events. The original work extracted data from the timestep immediately following breakup. We expanded the temporal sampling range, allowing the tool to extract data from the timestep immediately preceding breakup in addition to the timestep following breakup. This improvement is intended to provide insight into the conditions which lead to breakup and to the types of structures which the conditions produce. To highlight some potential uses of the tool, the present work extracted the locations of the events, the gas and liquid velocities, droplet volume, and droplet shape information.

For the purposes of this study, liquid and gas velocities are calculated as volume averaged velocities,
\begin{align}
    \bm{U}_\mathrm{liquid} &= \frac
    {\sum_{i=1}^{N_{s}}{\bm{u}_{i} \mathcal{V}_{\mathrm{cell},i} ~\alpha_i}}
    {\sum_{i=1}^{N_{s}}{\mathcal{V}_{\mathrm{cell},i} ~\alpha_i}}\\
    \bm{U}_\mathrm{gas} &= \frac
    {\sum_{i=1}^{N_{s}}{\bm{u}_{i} \mathcal{V}_{\mathrm{cell},i} ~(1-\alpha_i)}}
    {\sum_{i=1}^{N_{s}}{\mathcal{V}_{\mathrm{cell},i} ~(1-\alpha_i)}}
\end{align}
respectively. These values are calculated by considering all cells that contain a liquid structure and all cells adjacent to those ($N_s$ total cells). In the expressions, $\bm{u}_{i}$ is the velocity vector, $\mathcal{V}_{\mathrm{cell},i}$ is the cell volume, and $\alpha_i$ is the liquid volume fraction in the ${i}^\mathrm{th}$ cell. This method for calculating the velocities is preliminary and does not fully capture the local flow field dynamics. Future work focused on extracting topological data from the flow field will work to better quantify these values.

\subsection{Graphical Database}
The primary goal of the tool is to make data from atomization simulations more accessible. Special consideration was given to ensuring that extracted data was stored in an efficient and easily queriable format. We opted to store the atomization data in a Neo4j graphical database. Graph databases are commonly used in the corporate sector for companies to create connections between users and products through paths, which can reveal patterns in group dynamics and buying trends. The same principles can be applied to atomization simulations through the construction of paths which connect droplets through breakup and coalescence events. Storing atomization data in a graph database format allows researchers to analyze the evolution of liquid as it moves from the liquid core to small droplets in an atomizing flow. In this work, the Neo4j database is used to store the atomization data. Data from the simulation in this work is initially written to a CSV file. Then, rather than simply storing the data as a disconnected list of events, the CSV file is uploaded to Neo4j. The graph format presents an easy way to connect related events and develop a droplet ancestry. Additionally, storing data in graph databases provide a unique way to visualize atomization data and a novel method for studying the atomization process.
\begin{figure}[htbp]
 \centering
 \includegraphics[width=1\textwidth]{neo4j_full2.png}
 \caption{Example of data points in Neo4j graphical database. Red nodes represent some droplets which broke up 4 times throughout the simulation. Other colors represent intermediate breakups (blue = primary, yellow = secondary, green = tertiary). Lines connect related droplets, i.e., droplets which split from each other. The liquid core is present at the center of the image. Each node contains relevant statistics from the breakup event.}
 \label{fig:Neo4j}
\end{figure}

Fig.~\ref{fig:Neo4j} displays how data are stored in a Neo4j database. This image represents some of the data from the simulation described in Section 3. The nodes are droplets produced from breakup events. The breakup events are represented by the lines connecting the nodes. The node colors represent the number of times they have broken up. The liquid core is in the center of the image and the red droplets represent the sixth breakup events. This image is made up of all the droplets which broke up 6 times during the first 2.3 $\mu$s of the simulation, and all the intermediate droplets between the coherent liquid core and these sixth breakup droplets. Within each node is stored the extracted data from the point and time in the simulation in which that droplet broke up and became independent. So, this database is saturated with data, but organized in a logical and accessible format.

The graph database can be efficiently edited, reorganized, and parsed with the Cypher Query Language. It is syntactically similar to SQL, but designed to specifically query nodes and relationships and the paths formed with these components. As mentioned above, data is imported to Neo4j via a CSV file, organized in such a way that each row is an independent liquid breakup or coalescence event. These rows are imported into Neo4j as droplet nodes and breakup and coalescence relationships are created. See \ref{appendixb} for a detailed description of the process to import data into Neo4j. Following data import, Cypher can be used to further organize the data and/or query the data to analyze the atomization process. Below are some examples of the capabilities of Cypher in the present application.

\begin{lstlisting}
// Rename the liquid core
MATCH(n:droplet)     **// Find nodes with the "droplet" label**
WHERE n.Event='None' **// Node with Event = none is the liquid core**
CALL apoc.refactor.rename.label("droplet","core",n)  **// Rename node**
YIELD committedOperations
RETURN committedOperations

// Identify and rename primary droplets
MATCH (n:droplet)
WHERE n.OldLID = 1   **// Find droplets which broke off of the liquid core**
WITH collect(n) as p **// Compile a list of nodes matching criteria**
CALL apoc.refactor.rename.label("droplet","primary",p)
YIELD committedOperations
RETURN committedOperations

// Rename secondary,tertiary, etc. (repeat until no "droplet" nodes remain)
MATCH (n:droplet),(d:primary)
WHERE n.OldLID = d.NewLID **// Find "droplet" nodes that split from "primary"**
WITH collect(n) as p
CALL apoc.refactor.rename.label("droplet","secondary",p)
YIELD committedOperations
RETURN committedOperations
\end{lstlisting}

This portion of code first identifies the core by finding the only node which has an event not equal to either breakup or coalescence. Then, it renames this droplet "core". Note: the "Awesome Procedures on Cypher" (APOC) library must be enabled to access the renaming features. Next, primary droplets are identified by finding all the droplet nodes which previously had a $\mathcal{L}$ (LID) equal to one (the liquid core has $\mathcal{L} = 1$). These droplets are renamed "primary". After primary droplets are identified, secondary droplets can be identified using a similar process. The algorithm loops through all droplets and looks for nodes which previously had $\mathcal{L}$ equal to the $\mathcal{L}$ of a primary droplet. This process can then be repeated for each subsequent breakup event until no more "droplet" nodes remain. This is very useful in order to analyze the evolution of droplets as they break up further. This will be analyzed in more detail in the Results section.

\begin{lstlisting}
// Collect all breakup paths of droplets which break up six times
MATCH (n:sixth), (c:core), p = shortestPath((c)-[:Split*]->(n))
RETURN [d in nodes(p)| d.Volume]
\end{lstlisting}

The "shortestPath" function is used to extract the paths between specific nodes. In the above example, the shortest path from the liquid core to droplets that broke up six times is queried. Then, the volume of each droplet in that path is returned. This feature provides a simple method to extract useful relational data from the breakup process.

In addition to the Neo4j graphical database system and the Cypher Query Language, we utilized Python for analysis of results. Both programs have benefits and drawbacks, which led to the utilization of both. Neo4j allowed us to build paths very easily between droplets which break up from each other or coalesce together. From this information, we can easily output a CSV which displays how statistics evolve with breakup events or coalescence events. Following the output of these CSV files, it is simple to analyze the data using Python's data science libraries. Additionally, multiple options exist to query or create a graph databases through Python. These include the official Neo4j driver for Python and the py2neo Python library. These prove to be very powerful tools, because they allow users to combine python loops, if statements, and data science libraries with the unique organizational system of Neo4j. The py2neo library was used in the present work to build the time-series plots in section~\ref{sec:time}. See appendix~\ref{appendixc} for an example python script querying a Neo4j database.

\section{Application of Extraction Tool on Diesel Jet}

\subsection{Simulation Setup}
To test the utility of the updated tool we ran a simulation inspired by a diesel injector.  The simulation is the same as that described in Rubel and Owkes~\cite{RubelOwkes2019} and consists of a turbulent liquid injection into a quiescent air.  The turbulence was computed from a preliminary turbulent pipe-flow simulation and the velocity field was stored and used as the inlet boundary condition for the liquid jet.  Table~\ref{Table:3Djetnum} provides the parameters of the simulated jet. The simulation was run on 160 processors on the Hyalite High-Performance Computing Cluster at Montana State University. The dimensions of the computational mesh are $\mathcal{N}_x\times\mathcal{N}_y\times\mathcal{N}_z=1024\times128\times128$. Note that the resolution of this simulation is not fine enough to accurately capture small scale interface features or the smallest scales of turbulence. However, this simulation is sufficient to demonstrate the efficacy of the proposed tool. Future work will focus on applying the tool to a higher resolution simulation to more accurately extract information on physical phenomena.

\begin{table}[htbp]
\begin{center}
\begin{tabular}{ |c|c|c| }
 \hline
 Number & Definition & Value \\
 \hline
 Bulk Reynolds number & $\rho_{l}U_{\mathrm{jet}}D_{\mathrm{jet}}/\mu_{l}$ & 25,000 \\
 Bulk Weber number & $\rho_{l}U_{\mathrm{jet}}^2D_{\mathrm{jet}}/\sigma_{l}$ & 10,000 \\
 Density ratio & $\rho_{l}/\rho_{g}$ & 40 \\
 Viscosity ratio & $\mu_{l}/\mu_{g}$ & 1.67\\
 Domain length & $L_{x}/D_{\mathrm{jet}}$ & 60 \\
 Domain widths & $L_{y,z}/D_{\mathrm{jet}}$ & 7.5 \\
 Cells across diameter & $D_{\mathrm{jet}}/\Delta x$ & 17.06 \\
 CFL number & $|u|_{\mathrm{max}}\Delta t/\Delta x$ & 0.4 \\
 \hline
\end{tabular}
\end{center}
\caption{Non-dimensional parameters used in the diesel jet simulation.}
\label{Table:3Djetnum}
\end{table}

\begin{figure}[htbp]
 \centering
 \includegraphics[width=1\textwidth]{Jet_white_2.png}
 \caption{Rendering of the liquid jet run in simulation. $\mathcal{L}$ values at final timestep represented by colors. Rendered using VisIt~\cite{HPV:VisIt}.}
 \label{fig:jet}
\end{figure}

\subsection{Addressing Fictitious Events Confirmation}
Section~\ref{sec:fict_events} discusses an update in the coalescence and breakup identification algorithm, which is intended to prevent the occurrence of fictitious events. In the original work~\cite{RubelOwkes2019}, the fictitious events coincided with every breakup within the simulation, which forced the authors to institute a delay in the breakup identification portion of the code. This meant they only identified breakup events every 10 timesteps. In the present work, breakup identification is performed on every timestep. It is reasonable in complex flow systems that droplets may breakup and then re-merge with their parent droplet, so the algorithm is not intended to prevent this phenomena from occurring altogether.

To confirm that the updated method is producing reasonable results, we utilized Neo4j to identify all the droplets which split and re-merged with their parent droplet. Then, the difference in time between the coalescence event and the initial breakup event was calculated.
17.93\% of all droplets re-merged with their parent droplet, with the majority of those being primary droplets coalescing with the liquid core. This is a well documented mechanism~\cite{Shinjo2010} and expected in a round jet injected into quiescent gas. Analyzing only secondary droplets, we found that 9.31\% of secondary droplets re-merged with their parent droplets. The fictitious events were documented to occur within one or two timesteps of breakup, but in the present simulation, the average time between coalescence and breakup is $61.7$ $\mu$s which translates to about 300 timesteps. Given these statistics, we conclude that the updates made to the tool in the present work remedied the fictitious breakup and coalescence issue addressed in the original work and the identified events in the present work are consistent with physical events in the spray.

\subsection{Secondary Atomization Analysis}
Common descriptions of breakup regimes within an atomizing system involve discussion of primary atomization, i.e. droplets splitting from the liquid core, and secondary atomization, i.e., all further breakup. As seen in Fig.~\ref{fig:sec_hist}, many droplets are created after multiple breakup events.  Further analysis shows that 53.6\% of the total number of droplets formed in the test diesel jet broke up three or more times, with 17.0\% percent breaking up at least five times. At later times in the simulation, this percentage becomes more pronounced. The yellow curve with triangle markers in Fig.~\ref{fig:sec_hist} is the distribution of droplets at the final timestep. 67.7\% of drops present at this timestep broke up 3 or more times. Thus, a majority of the breakup in the simulation occurred after secondary breakup. This indicates that the final spray field is heavily dependent on mechanisms within these later atomization regimes, which are all generally grouped in secondary atomization.

\begin{figure}[htbp]
 \centering
 \includegraphics[width=1\textwidth]{secondary_atomization_analysis.png}
 \caption{The total number of droplets associated with each breakup stage from all breakup events within the diesel jet simulation.}
 \label{fig:sec_hist}
\end{figure}

\subsection{Coalescence Analysis}
A trend extracted from the test case simulation is the prevalence of coalescence events. This matches findings in Prakash et al.~\cite{Prakash2019}, who identified a significant number of coalescence events within their simulation. The present simulation produced 76,608 breakup events and 57,908 coalescence events. Table~\ref{Table:merge} displays the average number of coalescence events for a given number of breakup events. These data were calculated using a path finding algorithm in Neo4j. The algorithm finds the shortest path of breakup events from droplets at the given breakup stage back to the liquid core and counts the number of coalescence events along that path. The average number of coalescence events shown in the table is the average amongst all droplets that have broken up the associated number of times.  The data in Table~\ref{Table:merge} indicates that the proportion of coalescence events increases as the jet evolves and more breakup occurs. This makes intuitive sense, because as the jet evolves, a high droplet density cloud develops, which leads to more droplet collisions.

\begin{table}[htbp]
\begin{center}
\begin{tabular}{ |c|c|c| }
 \hline
 \thead{Breakup Stage} & \thead{Average Number of\\Coalescence Events} & \thead{Ratio of Coalescence\\to Breakup}\\
 \hline
 Secondary & 1.077 & 0.5385\\
 \hline
 Third & 1.899 & 0.6330\\
 \hline
 Fourth & 2.898 & 0.7245\\
 \hline
 Fifth & 3.974 & 0.7948\\
 \hline
 Sixth & 5.241 & 0.8735\\
 \hline
 Seventh & 6.823 & 0.9747\\
 \hline
 Eighth & 8.276 & 1.035\\
 \hline
\end{tabular}
\end{center}
\caption{The average number of coalescence events between primary breakup and the breakup stage is listed in the second column. The third column displays the ratio of coalescence events to breakup events along breakup paths.}
\label{Table:merge}
\end{table}

\subsection{Local Flow Field Statistics}
Another use for this tool is to better understand how the local flow field is affecting the atomization process. This is information that is exceedingly difficult, if not impossible, to obtain through experimental methods. Moreover, gathering a sufficient sample of statistics on these local flow mechanisms has remained nonviable for those conducting numerical simulations. This work extracted preliminary data on the local flow field surrounding breakup events and calculated the resultant local Weber number, defined as $\mathrm{We}_\mathrm{local}={\rho_g U_s^2 L}/{\sigma}$, where $\rho_g$ is the gas velocity, $U_s$ is the slip velocity (defined as $U_s = |U_\mathrm{liquid}-U_\mathrm{gas}|$), $L$ is the characteristic length (in this case, it is the equivalent spherical diameter of the droplet prior to break up), and $\sigma$ is the surface tension coefficient.
A probability density function (PDF) of the logarithm of $\mathrm{We}_\mathrm{local}$ number is developed from every breakup within the simulation (Fig.~\ref{fig:weber_pdf}). The values are centered roughly about $\mathrm{We}_\mathrm{local}=0.37$. These are very small values of the Weber number, which likely indicates that aerodynamic forces are not major factors in breakup in this system. This is logical, because the simulation is liquid injected into quiescent air.

\begin{figure}[!htbp]
 \centering
 \includegraphics[width=1\textwidth]{Weber_pdf_all.png}
 \caption{A probability density function of the local Weber Number associated with every breakup in the simulation.}
 \label{fig:weber_pdf}
\end{figure}

\subsection{Atomization Evolution}
The evolution of local droplet characteristics and global jet development throughout the atomization process are important attributes, which can help elucidate the underlying physics of atomization. Understanding how droplets change from first breakup off of the liquid core to a final droplet in dilute spray could provide useful information to atomization model developers seeking to not only describe the final spray formation, but also the intermediate spray development. The tool provided in this work enables researchers to extract droplet statistics from throughout a simulation and analyze the evolution of these droplets and the system as a whole.

\subsubsection{Droplet Shape and Size Evolution}
Fig.~\ref{fig:diam_pdf} displays a probability density function of the equivalent spherical diameter of droplets as a function of the number of breakup events. A reasonable trend toward smaller and more uniformly sized droplets is seen, which provides confidence in the values extracted by the tool.
Fig.~\ref{fig:diam_change} displays the change in diameter of droplets between breakup stages. Notice that not all values are negative. This further confirms the prevalence of coalescence events within these systems. Positive values in the figure indicate that a significant portion of droplets undergo coalescence and increase in size between breakup events. Additionally, notice the tendency of droplets to change size less as more breakup events occur, with the PDFs becoming more narrow and centered around 0 $\mu$m.

\begin{figure}[!htbp]
 \centering
 \includegraphics[width=1\textwidth]{diameter_pdf.png}
 \caption{A probability density function of droplet diameter as a function of number of breakups}
 \label{fig:diam_pdf}
\end{figure}

\begin{figure}[!htbp]
 \centering
 \includegraphics[width=1\textwidth]{change_in_diameter.png}
 \caption{A probability density function displaying the change in diameter between breakup events.}
 \label{fig:diam_change}
\end{figure}

\subsubsection{Time Evolution} \label{sec:time}
Because the extraction tool samples from every timestep in which breakup occurs. Thus, once the jet develops, breakup occurs on every timestep, providing high resolution time-series statistics describing evolution of the jet. Fig.~\ref{fig:drop_vol_evo} shows an example of some of the time dependent statistics that were extracted in the present work. These plots describe the rate of droplet development in the simulation. It is clear from this plot that secondary atomization quickly dominates droplet production in this system. Future work will focus on improving the extracted quantities and further explore how atomizing systems develop through time.

\begin{figure}[!htbp]
 \centering
 \includegraphics[width=1\textwidth]{Drops_Created1.png}
 \caption{The cumulative number of primary, secondary, and total droplets. Note: secondary droplets in this image represent the classical definition of secondary droplets (i.e. all breakup following primary).}
 \label{fig:drop_vol_evo}
\end{figure}

\section{Conclusions}
An extraction tool for numerical simulations of atomization was developed and tested. Improvements upon the methodology originally proposed by Rubel and Owkes~\cite{RubelOwkes2019} were made. Improvements include 1) an updated algorithm to stop fictitious identification of merge and split events, 2) enabling a wider temporal sampling range, and 3) the extraction of more atomization statistics.

A diesel-type jet was simulated and droplet statistics were extracted. The tool displayed utility for extracting data relevant to atomization model developers as well as providing previously inaccessible information on the underlying mechanisms of atomization. The data extracted are preliminary, but show promise for the tool's utility in future high-resolution studies.

Future work will focus on the implementation of the tool into high-resolution atomization simulations to extract useful data to help better quantify atomization mechanisms on global and local scales. Special effort will be given to aiding reduced-order model developers to develop relevant statistics for creating more accurate models. Additionally, we will focus on the improvement of sampling methods, particularly the local flow field to attempt to extract not only velocity magnitudes, but also topological flow data to better understand the small-scale turbulence that effects liquid breakup.

\section{Acknowledgements}
This material is based upon work supported by the National Science Foundation under Grant No.~1749779.
\vskip 4mm
Computational efforts were performed on the Hyalite High Performance Computing System, operated and supported by  University Information Technology Research Cyberinfrastructure at Montana State University.

\bibliographystyle{unsrt}

\section{}
\label{}

\bibliographystyle{elsarticle-num}