Year: 2021

Abstract:
3D snapshot microscopy enables fast volumetric imaging by capturing a 3D volume in a single 2D camera image and performing computational reconstruction. Fast volumetric imaging has a variety of biological applications such as whole brain imaging of rapid neural activity in larval zebrafish. The optimal microscope design for this optical 3D-to-2D encoding is both sampleand task-dependent, with no general solution known. Deep learning based decoders can be combined with a differentiable simulation of an optical encoder for end-to-end optimization of both the deep learning decoder and optical encoder. This technique has been used to engineer local optical encoders for other problems such as depth estimation, 3D particle localization, and lensless photography. However, 3D snapshot microscopy is known to require a highly non-local optical encoder which existing UNet-based decoders are not able to engineer. We show that a neural network architecture based on global kernel Fourier convolutional neural networks can efficiently decode information from multiple depths in a volume, globally encoded across a 3D snapshot image. We show in simulation that our proposed networks succeed in engineering and reconstructing optical encoders for 3D snapshot microscopy where the existing state-of-the-art UNet architecture fails. We also show that our networks outperform the state-of-the-art learned reconstruction algorithms for a computational photography dataset collected on a prototype lensless camera which also uses a highly non-local optical encoding.

heading: 1 INTRODUCTION
text: Volumetric imaging has been a valuable technique for measuring neural activity at single neuron resolution across the brain [2, 22, 36]. We are particularly motivated by the problem of imaging neuronal activity across the whole brain of transparent model organisms such as larval zebrafish (Danio rerio) at higher volume rates, to enable richer studies of neural dynamics and signal propagation. Existing approaches to whole-brain imaging involve sequential acquisition of stacks of 2D images [2], typically resulting in 0.5Hz - 2Hz volume capture rates, limited by the number of 2D depth plane acquisitions. However, since neural activity occurs at orders of magnitude faster timescales, faster volume imaging rates can provide understanding of neural computations with significantly more detail.
3D snapshot microscopy techniques can significantly speed up the imaging of a 3D sample volume using non-local optical encodings to combine information corresponding to different depth planes in the volume into a single 2D image, allowing for fast volume acquisition at the frame rate of the camera sensor followed by computational decoding to reconstruct the 3D volume [5, 3, 4, 7, 13, 16, 18, 17, 20, 27, 26, 30, 38, 37, 33]. The optimal optical design for a 3D snapshot microscope is unknown in general [5], and likely depends on the specific sample structure. Newly developed programmable microscopes with up to 106 free parameters, e.g. pixels on a spatial light modulator (SLM), enable the implementation of a rich space of microscope encodings, and present the possibility of direct optimization of snapshot microscope parameters specifically for particular classes of samples and imaging tasks.
Deep learning-based decoders combined with differentiable simulations of optical encoders have enabled end-to-end optimization of both the deep learning decoder and the optical encoder. This endto-end optical engineering approach has recently been applied to depth from defocus [6, 10, 14, 35], lensless photography [13], or particle localization microscopy [23]. However, all such applications
either use far fewer optical parameters or operate on relatively compact or local encodings due to small filters in their deep networks [6, 23, 35, 10, 14, 34, 13]. There are two barriers to optimizing more parameters and producing global versus local encodings: (1) the typical convolutional networks with small filters used by previous methods for computational reconstruction have a prior to decode locally, and therefore limit optical designs to local encodings, unlike the global PSFs used in snapshot microscopy and (2) the optimization of optical parameters involves backpropagation through both the deep learning reconstruction and a differentiable optical simulation, which is more expensive computationally for non-local optical encodings. We focus on two computational optics problems that demonstrate a superior deep network architecture that can successfully reconstruct from images encoded by non-local optical encoders and with orders of magnitude more optical parameters than previously attempted. To address the computational burden of simulating and optimizing non-local encoders, we also developed a parallel, multi-GPU optical simulation framework.
Problem statement We define an optical encoder as Mφ parameterized by φ and a computational decoder as Rθ parameterized by θ. We wish to develop a reconstruction network architecture for decoding non-local encodings in 2D images c produced by Mφ and also to enable the end-to-end optimization of φ to produce non-local encodings. We’ll consider two concrete applications: (1) 3D snapshot microscopy using an SLM-based programmable microscope (Figure 1A), where we wish to optimize both the reconstruction network Rθ and the optical encoding Mφ based on the SLM parameters in order to image 3D samples v and reconstruct 3D volumes v̂, and (2) lensless photography (Figure 4), where we wish to optimize only the reconstruction network Rθ to produce images of natural scenes imaged by a lensless camera M implementing a fixed optical encoding.

heading: 1.1 OUR CONTRIBUTIONS
text: 1. We developed a large-scale, parallel, multi-GPU differentiable wave optics simulation of a programmable microscope, based on a 4f optical model with a phase mask (φ) implemented using a spatial light modulator (SLM), described further in Appendix A.1. SLMs φ can have over 106 optimizable parameters.
2. We collected a large dataset of high resolution 3D confocal volumes of zebrafish larvae for the purpose of end-to-end sample-specific optimization of the parameters of the programmable 3D snapshot microscope.
3. We introduce an efficient new reconstruction network architecture for decoding from nonlocal optical encoders using very large global convolutions implemented via Fourier convolutions.
4. We show that our networks outperform the state-of-the-art deep networks for both volume reconstruction and microscope parameter optimization, and for image reconstruction from lensless photography [20].
5. Our method enables, for the first time, direct end-to-end optimization of highly non-local optical encoders in the space of SLM pixels with over 106 parameters. In simulation, we demonstrate the potential for significant improvements in imaging resulting from samplespecific engineering of optical encoders.

heading: 1.2 PRIOR WORK
text: Neural network architectures for computational imaging have all used convolution layers with small filters. End-to-end optimization of optical encoders have largely been performed with UNet-based architectures [6, 10, 14, 35, 13], with one method using ResNet-based architectures [23]. Such optimization has always led to local optical encodings. End-to-end optimization has never been attempted for large-field of view 3D snapshot microscopy due to the difficulty of simulating and reconstructing from non-local encoders. However, small filter convolutional deep networks have been used in a non-end-to-end manner to reconstruct volumes from 3D snapshot microscopes designed using microlens arrays [33, 34]. For photography and MRI, using deep learning combined with unrolling iterations of traditional deconvolution algorithms provides the benefits of fast amortized inference and higher quality reconstructions due to learning of structural priors [8, 9, 20, 34].
In our work, we demonstrate that convolution layers with large filters implemented efficiently in the Fourier domain enable the end-to-end learning of highly non-local optical encoders for 3D snapshot microscopy. Large convolution filters have been shown to be helpful for other computer vision
applications such as semantic segmentation and salient object detection [25, 21], and the Fourier domain parameterization of small filters has been described previously [28].
A pioneering strategy in 3D snapshot microscopy has been light field microscopy [16], which employs a microlens array at the microscope’s image plane to create subimages encoding both amplitude and phase of light [16, 1]. A variety of microlens-array-based light field microscopes have been used to perform whole-brain imaging [16, 27, 37, 36, 30, 26, 7, 12, 38]. [38] optimizes the placement of microlenses, but not in an end-to-end manner. Despite variation in design, microlens-based microscopes have, to various degrees, three main limitations that can be improved: 1) blocking or scattering of light between microlenses, causing light inefficiency, 2) not making use of all pixels on the camera to encode a 3D sample, leading to inefficient compression and suboptimal reconstructions v̂, and 3) a fixed optical encoding scheme.
An alternative to using microlenses is to implement a coded detection strategy using a phase mask or diffuser to spread light broadly across the camera sensor [3, 4, 5, 13, 17, 18, 38]. The designed phase masks can be implemented either by manufacturing a custom optical element or using a programmable SLM [10, 14, 35, 13, 38, 23]. Using a programmable element allows different microscope parameters to be used for different samples.

heading: 2 METHODS
text: We show our network architecture and an overview of autoencoder training both the microscope parameters φ and reconstruction network parameters θ in Figure 1A. The programmable microscope is simulated by a differentiable implementation of a wave-optics model of light propagation. We have selected a programmable microscope design based on pupil-plane phase modulation with a programmable spatial light modulator, for which imaging is well-approximated by a computationallyefficient convolution [11]. A detailed description of our simulation is provided in Appendix A.1. For lensless photography, there is no optical simulation because the images have been collected on a real camera.

heading: 2.1 FOURIERNET AND FOURIERUNET ARCHITECTURES FOR RECONSTRUCTION FROM OPTICAL ENCODERS
text: Because images created by optical encoders can potentially encode signals from the incoming light field to any location in the camera image in a non-local manner, it is essential that a reconstruction network have global context. Existing multi-scale architectures such as the UNet [29] can achieve global context, but at the expense of many computation layers with small convolution kernels which we believe have a local information bias which is inappropriate for computational optics. In this paper, we introduce relatively shallow architectures for computational imaging, which rely on convolution layers with very large filters, implemented efficiently in the Fourier domain [28].
FourierNet We propose a simple three layer convolutional network architecture with very large global convolutions at the very first layer, followed by two standard local convolutional layers (Figure 1B). We define a global convolution as a convolution with kernel size equal to the input image. Such a convolution achieves global context in a single step but is computationally expensive. Global convolutions are implemented more efficiently in the Fourier domain, yielding a speed up of two orders of magnitude. Due to the use of Fourier convolutions to enable global context, we call our architecture FourierNet. In contrast to a typical UNet which can contain many tens of convolution layers, the FourierNet is only three layers deep, which requires backpropagation through fewer layers compared to a typical UNet with the same receptive field.
FourierUNet We also propose a multi-scale variant of the FourierNet by bringing together elements of the multi-scale UNet and the single-scale FourierNet. Here, we take advantage of the fact that down-sampling in image space corresponds to a simple cropping operation in the Fourier domain, resulting in a band-limited computation of a feature map. We efficiently implement multi-scale global Fourier convolutions (Figure 1C) to replace the encoding/“analysis” pathway of a UNet. We then use the standard decoding/“synthesis” pathway of the UNet to combine the multi-scale features into a single 3D volume reconstruction (Appendix A.3, A.4). Thus we can study whether multi-scale features or global context is more important for decoding non-local optical encoders.
Fourier domain convolutions It is well-known that large kernel convolutions can be implemented more efficiently in the Fourier domain [24, 19, 32]. A naive implementation of global convolution requires O(N2) operations, where N is number of pixels in both the image and the kernel. An alternative global convolution implementation is to Fourier transform the input x and convolution kernel w, perform element-wise multiplication in Fourier space, and finally inverse Fourier transform, requiring onlyO(N logN) operations [19, 32]. Following [28], we store and optimize the weights in Fourier space W . This over-parameterization costs 8× the memory of an equivalent real valued large filter but saves the computational cost of Fourier transforming the real-valued weights (Appendix A.3). Thus a Fourier convolution is defined:
Re{F−1 {W F {x}}} (1) For image and kernel sizes of 256× 256 pixels, our implementation leads to nearly 500× speedup: standard PyTorch convolution takes 2860ms, while Fourier convolution takes 5.92ms on a TITAN X.
Multi-scale Fourier domain convolutions It is well-known [24] that downsampling corresponds to cropping in the Fourier domain. Thus the Fourier convolution can be extended to efficiently produce multi-scale feature representations in one step (Figure 1C). We define our multi-scale Fourier convolution as {
Re{F−1 {W1 crop1 [c]}}, . . . ,Re{F−1 {Wn cropn [c]}} }
(2) where subscript denotes scale level (higher subscript indicates lower spatial scale/more cropping in Fourier space) and we precompute c := F {c} once.

heading: 2.2 PHYSICS-BASED AUTOENCODER FOR SIMULTANEOUS ENGINEERING OF MICROSCOPE ENCODER AND RECONSTRUCTION NETWORK DECODER
text: We describe the imaging process as the following transformation from the 3D light intensity volume of the sample v to the 2D image formed on the camera c:
µc = Mφ(v) (3) c = max ([µc + √ µc ] , 0), ∼ N (0, 1) (4)
where Mφ denotes the microscope parameterized by a 2D phase mask, φ. This phase mask φ describes the 3D-to-2D encoding of this microscope model completely. A Poisson distribution
with mean rate µc describes the physics of photon detection at the camera, but sampling from this distribution is not differentiable. We approximate the noise distribution with a rectified Gaussian. We include details on Mφ in Appendix A.1 [11]. Jointly training reconstruction networks and microscope parameters involves image simulation, reconstruction, then gradient backpropagation to update the reconstruction network and microscope parameters. Details on parallelization, planewise reconstruction networks, and planewise sparse gradients are provided in Appendix A.2.
Our loss function computes the normalized mean squared error (MSE) LHNMSE between the high pass filtered sample and reconstruction weighted by the normalized MSE LNMSE between the original sample and reconstruction (to include low frequency content in the objective). Formally, our loss function L(v, v̂) is defined:
L(v, v̂) = LHNMSE(v, v̂) + βLNMSE(v, v̂) (5) LHNMSE(v, v̂) = E [ (H(v)−H(v̂))2 ] E(H(v)2) , LNMSE(v, v̂) = E [ (v − v̂)2 ] E(v2) (6)
where H(·) denotes high pass filtering and E(·) denotes computing the mean. Both loss terms are normalized as shown to reduce variance in L, which can otherwise take on large magnitude values and cause training instability. For our experiments, we set the weight β for the LNMSE term to 0.1.

heading: 3 RESULTS
text: Larval Zebrafish Dataset High resolution volumes of transgenic larval zebrafish whole brains expressing nuclear-restricted GCaMP6 calcium indicator in all neurons were imaged using a confocal microscope. These images are representative of brain-wide activity imaging. We train on 58 different zebrafish volumes (which we augment heavily) and test on 10 held-out volumes. For all experiments, we downsample the high resolution confocal data to (1.0 µm z, 1.625 µm y, 1.625 µm x). We created 4 datasets from these data called Type A, B, C, and D corresponding to imaging different fields of view. Full specifications for these datasets are in Table 4 (Appendix A.2). For Figure 2, we restrict the field of view to (200 µm z, 416 µm y, 416 µm x) with a tall cylinder cutout of diameter 193 µm and height 200 µm and image with 256× 256 pixels on the simulated camera sensor. We call this setting Type D. Figure 3 and Table 2 show our larger experiments with 512× 512 simulated camera pixels, with a field of view of (250 µm z, 832 µm y, 832 µm x) and sample types Type A, Type B, and Type C.
DiffuserCam Lensless Mirflickr Dataset We also test reconstruction performance on experimental computational photography data1 from [20] (Figure 4). This is a dataset constructed by displaying RGB color natural images from the MIRFlickr dataset on a monitor and then capturing diffused images by the DiffuserCam lensless camera. The dataset contains 24,000 pairs of DiffuserCam and ground truth images. The goal of the dataset is to learn to reconstruct the ground truth images from the diffused images. As in [20], we train on 23,000 paired diffused and ground truth images, and test on 999 held-out pairs of images.

heading: 3.1 FOURIERNETS OUTPERFORM UNETS FOR ENGINEERING NON-LOCAL OPTICAL ENCODERS
text: We compare optimizing microscope parameters φ with two neural networks: 1) using our FourierNet with 2D convolutions (FourierNet2D) and 2) using a vanilla UNet with 2D convolutions (UNet2D).
1Publicly available: https://waller-lab.github.io/LenslessLearning/dataset.html
Training is a two step process in which a planewise reconstruction network with fewer parameters and 2D convolutions is used during optimization of φ, then once φ is fixed a more powerful reconstruction network with 3D convolutions is used. We find this scheme achieves better reconstruction quality because the reconstruction network does not need to constantly adapt to a changing optical encoding (Appendix A.3). We train these microscopes and reconstruction networks on Type D, where samples are tall cylindrical cutouts of zebrafish with diameter 193 µm and height 200 µm. Samples are imaged on a camera with 256× 256 pixels (Figure 2). FourierNet2D has 2 convolution layers (with 99.8% of its kernel parameters in the initial Fourier convolution layer), while UNet2D has 32 convolution layers (kernel parameters approximately uniformly distributed per layer). FourierNet2D and UNet2D both have ∼ 4× 107 parameters; FourierNet3D has ∼ 6× 107 parameters vs. ∼ 108 for UNet3D. Architecture details are in Appendix A.3, A.4.
Simulated camera images (Figure 2) show that the UNet microscope does not make sufficient use of camera pixels, producing only a single view of the sample. We believe this is due to a local information prior in the small kernels of UNets. The FourierNet microscope uses more camera pixels and performs better than the UNet microscope for reconstruction (Figure 2, quantified in Table 1).
Times in Table 1 are for only the forward pass on a single GPU; one training iteration on 8 GPUs takes ∼0.4 seconds for FourierNet3D and ∼0.8 seconds for UNet3D (Appendix A.3). Both reconstruction networks must reconstruct from images that have a compressed encoding of 3D information, but the FourierNet2D is clearly more effective than the UNet2D at optimizing this encoding.

heading: 3.2 FOURIERNETS OUTPERFORM UNETS FOR 3D SNAPSHOT MICROSCOPY VOLUME RECONSTRUCTION
text: We can determine which architecture is better for volume reconstruction by choosing fixed microscope parameters and varying the architecture. In Table 1, we compare results using a FourierNet with 3D convolutions (FourierNet3D), a FourierUNet with 3D convolutions (FourierUNet3D), and a vanilla UNet with 3D convolutions (UNet3D). Architecture details are in Appendix A.3, A.4.
Reconstruction results in Table 1 compare normalized MSE LHNMSE between the high pass filtered sample and high pass filtered reconstruction, the multiscale structural similarity MS− SSIM between the sample and reconstruction, and finally the peak signal-to-noise ratio PSNR. We also visualize reconstruction results for a volume in the head of a zebrafish in Figure 2. The best reconstruction networks are (equally) FourierNet3D and FourierUNet3D. The UNet3D reconstruction networks (using either microscope) fall significantly behind the FourierNet3D/FourierUNet3D reconstruction networks in all metrics.

heading: 3.3 ENGINEERED OPTICAL ENCODING DEPENDS ON SAMPLE STRUCTURE
text: To explore the effect of sample structure and size on optimized φ and the resulting reconstruction performance, we optimized φs for three different sample types: 1) Type A, samples with a short cylinder cutout of 386 µm diameter and 25 µm height, 2) Type B, samples with a tall cylinder cutout of 386 µm diameter and 250 µm height, and 3) Type C, samples without any cutout of dimension (250 µm z × 832 µm y × 832 µm x) (Table 2). All sample types were imaged with 512 × 512 pixels on the simulated camera. We then tested reconstruction performance on all combinations of optimized microscopes and samples, as shown in Table 2 and visualized in Figure 3 for Type B. We include architecture details in Appendix A.3, A.5.
We see in Table 2 that for all sample types, highest performance is achieved using the microscope optimized for that particular sample. These results show that our optimization process produces optical encoders that are more optimal for a specific type of sample over others.

heading: 3.4 FOURIERNETS OUTPERFORM STATE-OF-THE-ART FOR RECONSTRUCTING NATURAL IMAGES CAPTURED BY DIFFUSERCAM LENSLESS CAMERA
text: We compare our FourierNet architecture to the best learned method from [20] using an unrolled ADMM and a denoising UNet, as well as to a vanilla UNet from [20]. Architecture details are in Appendix A.3, A.6. We can see that both of our methods visually outperform the methods from [20] in Figure 4. Table 3 shows that our methods outperform the others on all metrics. In [20], a combined loss using MSE and learned perceptual loss (LPIPS) is suggested to improve visual reconstruction quality. Table 3 shows that training our FourierNet on both MSE and LPIPS results in the lowest
LPIPS loss across all methods, but does not offer better visual performance than our MSE-only model as seen in Figure 4. FourierUNet does not show improved results over FourierNet.

heading: 3.5 GLOBAL RECEPTIVE FIELD IS MORE IMPORTANT THAN MULTISCALE FEATURES
text: UNets are effective because: (1) features are computed at multiple scales and (2) large receptive fields are achieved in few layers. FourierUNets allowed us to decouple these two explanations because the receptive field is global in a single layer. We see on both our microscopy dataset (which does not have multiscale structure) and the lensless photography dataset (which does have multiscale structure) that the FourierUNet does not improve upon the FourierNet. Thus we see that it is more important for decoding from non-local optical encoders to have a global receptive field than multi-scale features.

heading: 4 DISCUSSION
text: Summary We have presented FourierNets for decoding from highly non-local optical encodings. We have demonstrated that FourierNets enable end-to-end optimization of non-local optical encodings
where UNets fail for simulated 3D snapshot microscopy, and also that FourierNets beat state-of-the-art decoders for lensless photography. FourierNets provide fast amortized inference and reconstruction times many orders of magnitude faster than traditional iterative reconstruction algorithms [20] and higher quality reconstructions due to effective learning of sample structural priors. Generally, our global kernel architecture using Fourier convolutions could be applied to other problems where global integration of features is necessary (though we have focused on computational microscopy and photography where such global mixing of information dominates).
Limitations Our engineered optical encoders have not yet been experimentally tested on a programmable microscope. This will require measuring and accounting for system aberrations and calibration to implement the engineered optical encodings. However, the parameters of even an imperfectly implemented optical encoder can be perfectly measured and used to train reconstruction networks to near optimal performance. Our claim that sample type-specific optical encoders perform better on their respective samples has only been evaluated in simulation. Further, imaging a different sample class (i.e. with different spatiotemporal statistics) will require retraining at least the reconstruction network, and ideally also the optical encoder. Our autoencoder training style requires simulation of imaging with gradients, which can be memory expensive and limits the total field of view. Our training approach requires multi-GPU parallelization (described in Appendix A.2), and in practice we optimized microscopes using 8 Quadro RTX 8000 GPUs or 4 RTX 2080 Ti GPUs for the larger and smaller experiments, respectively. Training a microscope and decoder requires at least two weeks of total training. Our Fourier convolution layers require a fixed input size to truly be global, though this could be addressed by resampling.
Reproducibility We train on 8 Quadro RTX 8000 GPUs for the largest experiments, and have described our pre-processing, training, and testing procedures in Appendix A.2, A.3, A.4. We will also release our simulation software, training scripts, and our datasets upon publication.

heading: A APPENDIX
text: A.1 FORWARD SIMULATION OF PROGRAMMABLE 3D SNAPSHOT MICROSCOPE
Here we describe our wave optics simulation of the microscope Mφ, which we model as a 4f system [11]. The 4f optical system consists of two lenses, the first spaced one focal length from the object plane and the second spaced one focal length away from one focal length beyond the first lens (Figure 5). In between these two lenses, we can place a phase mask to manipulate the light field before passing through the second lens and forming an image on the camera sensor.
We are concerned here with fluorescence microscopy, meaning that the sources of light that we image are individual fluorescent molecules, which we can model as point emitters. Because these molecules emit incoherent light, the camera sensor in effect sums the contributions of each point emitter. In order to model such an imaging system, we first need to address modeling a single point emitter’s image on the camera.
We can analytically calculate the complex-valued light field one focal length after the first lens (which we call the pupil plane) due to a point source centered at some plane z (where z is a distance from the object plane z = 0). If the point source were centered (x = 0, y = 0) in the object focal plane z = 0, we would have a plane wave at the pupil plane, but for the more general case of a point source at an arbitrary plane z relative to the object plane z = 0, we can analytically calculate the complex-valued light field entering the pupil plane:
upoint(k; z) = exp
[ i2πz √(n λ )2 − ||k||22 ] (7)
where upoint is the incoming light field entering the pupil due to a point source centered in the plane at z, k ∈ R2 denotes frequency space coordinates of the light field in the pupil plane, n is the refractive index, and λ is the wavelength of light [11].
In this pupil plane, we can then apply a phase mask φ to the light field, which is modeled as a multiplication of upoint(k; z) and eiφ(k), the complex phase of the pupil function. The light field exiting the pupil is therefore described by
upupil(k; z) = upoint(k; z)p(k) (8)
where p(k) is the pupil function, composed of an amplitude a(k) and phase φ(k):
p(k) = a(k)eiφ(k) (9)
a(k) =
{ 1 ||k||2 ≤ NA λ
0 ||k||2 > NA λ
(10)
where NA is the numerical aperture of the lens [11].
The light field at the camera plane can then be described by a Fourier transform [11]:
ucamera = F {upupil} (11)
The camera measures the intensity of this complex field:
s(x; z) = |ucamera(x; z)|2 (12) where x ∈ R2 denotes spatial coordinates in the camera plane [11]. We can call this intensity s the point response function (PRF). If the shape of the PRF is translationally equivariant in x, meaning that moving a point source in-plane creates the same field at the camera, just shifted by the corresponding amount, then we call this PRF a point spread function (PSF). Note that moving the point source in z will not give the same shape, which allows our system to encode depth information through the PSF [5].
In order to avoid edge effects during imaging, we simulate the PSF at a larger field of view, then crop and taper the edges of the PSF: staper = crop[s] t (13) where t is a taper function created by taking the sigmoid of a distance transform divided by a width factor controlling how quickly the taper goes to 0 at the edges and denotes elementwise multiplication. We intentionally simulate a larger field of view than the sample in order to avoid edge artifacts. The purpose of the crop[·] is to cut the PSF to the correct field of view. The purpose of the tapering is to remove artifacts at the edges of the cropped PSF. After we compute this cropped and tapered PSF, we also downsample staper to the size of the data v in order to save memory.
Imaging is equivalent to the convolution of the incoming light field volume intensity v and the cropped and tapered PSF staper for a given plane. At the camera plane, the light field intensity is measured by the camera sensor. Therefore, we can describe the forward model as the following convolution and integral over planes:
µc(x) = ∫∫ v(τx; z)staper(x− τx; z) dτx dz (14)
We then model shot noise of the camera sensor to produce the final image c, for which the appropriate model is sampling from a Poisson distribution with a mean of µc [11]:
c ∼ Poisson (µc) (15) However, because we cannot use the reparameterization trick to take pathwise derivatives through the discrete Poisson distribution, we instead approximate the noise model with a rectified Gaussian distribution:
∼ N (0, 1) (16) c ≈ max ([µc + √ µc ] , 0) (17)
We now turn our attention to selecting the number of pixels used in the phase mask, i.e. the number of parameters for Mφ. We first need to determine the pixel size for Nyquist sampling the image plane with an objective of a given NA (numerical aperture). For a given pixel size ∆x, we know that in frequency space coordinates we will have a bandwidth of 1∆x , spanning − 1 2∆x to 1 2∆x . Because we must have
||k||2 ≤ NA
λ (18)
we know that the Nyquist sampling pixel size is given by
∆x∗ = λ
2NA . (19)
Therefore, in the image plane, for a desired field of view L we must have at least
N∗ = L
∆x∗ (20)
pixels. The discretization in the pupil plane will be the same, which means we will need to have at least N∗ pixels in the pupil plane to achieve the appropriate light field in the image plane. For our settings of NA = 0.8, λ = 0.532µm, and L = 823µm, we have ∆x∗ = 0.3325µm and N∗ = 2476 pixels. Thus, a reasonable choice is ∆x = 0.325µm andN = 2560 pixels. Note that these simulation parameters are independent of the camera pixels; we have determined only how many pixels must be used in the phase mask in order to ensure our PSF can occupy the full field of view. The camera sensor can sample the field at the image plane at an independent pixel size.
A.2 TRAINING PSFS AND VOLUME RECONSTRUCTION NETWORKS
Given a simulation of imaging, we can define two modes of autoencoder training: (1) jointly training the phase mask parameters φ and weak reconstruction networks in order to learn a good PSF for a particular class of samples (i.e. samples with the same spatiotemporal statistics), and (2) training a stronger reconstruction network only with a fixed, pre-trained φ.
Definition of terms For both cases of training, the general framework is to simulate imaging using confocal volumes of pan-neuronal labeled larval zebrafish, reconstruct from the simulated image, then update the reconstruction network and, if desired, the microscope parameters. We will define the microscope parameters as φ and the reconstruction network parameters as θ for any reconstruction network Rθ(c) where Rθ maps 2D images to 3D volume reconstructions. For our training algorithms listed below, we also define: D our dataset, v a ground truth volume, v̂ a reconstructed volume, L a computed loss, zs a list of z plane indices that will be imaged/reconstructed, αφ the learning rate for the microscope parameters, αθ the learning rate for the reconstruction network parameters, and β the weight of the non-high pass filtered component of the loss. When selecting a random ground truth volume, we also perform random shift, rotation, flip, and brightness augmentations.
Microscope simulation parameters When simulating the zebrafish imaging, we use a wavelength of 0.532 µm for all simulations. The NA of our microscope is 0.8. The refractive index n is 1.33. We downsample all volumes to (1.0 µm z, 1.625 µm y, 1.625 µm x). We use a taper width of 5 for all
simulations, and simulate the PSF at 50% larger dimensions in x and y. The resolution of the camera (for all zebrafish datasets) is also (1.625 µm y, 1.625 µm x).
Initialization of φ For Type A, B, and D, we initialize φ to produce a PSF consisting of 6 pencil beams at different locations throughout the depth of the volume, with the centers of these beams arranged in a hexagonal pattern in x and y. Because our optimizations generally find PSFs with many pencils, we find that initializing with such a pattern helps to converge to a more optimal PSF (data not shown).
For Type C, we instead initialize with a single helix spanning the depth of the volume (the “Potato Chip” from [5]), which seems to find a local minimum for φ that produces a PSF with more pencils (and therefore views in the camera image).
Data settings and augmentation for Type A, B, C, D Using our total 58 training zebrafish volumes and 10 testing zebrafish volumes (imaged through confocal microscopy), we crop in four different ways to create four different datasets. For training volumes, we crop from random locations from each volume as a form of augmentation. For testing, we crop from the same location. Physically, these crops correspond to either placing a circular aperture before light hits the 4f system or changing the illumination thickness in z, because samples would be illuminated from the side in a real implementation of this microscope. We model these by cropping cylinders (or cubes if there is no aperture) of different diameters and heights. We show details for all types Type A, B, C, D in Table 4,
where the diameter of the cylinder is labeled “Aperture Diameter” and the illumination thickness is labeled “Height”.
We augment our volumes during training by taking random locations from these volumes, randomly flipping the volumes in both z and y, and also randomly rotating in pitch, yaw, and roll. Most importantly, we also randomly scale the brightness of our samples and add random background levels which serve to adjust the signal-to-noise ratio (SNR) of the resulting simulated images. The only exception to these augmentations is Type C, where we set all the volumes to the same in-plane vertical orientation (while still applying rotation augmentations in pitch and roll).
Parallelizing imaging and reconstruction Furthermore, because this simulation can become too expensive in memory to fit on a single device, we generally perform the simulation, reconstruction, and loss calculation in parallel for both training modes. Therefore, any variable that has a s subscript refers to a list of chunks of that variable that will be run on each device. A j superscript indicates a particular chunk for GPU j. For example zs is a list of plane indices to be imaged/reconstructed, and zjs is the j
th chunk of plane indices that will be imaged/reconstructed on GPU j. We denote parallel for any operations that are performed in parallel and scatter for splitting data into chunks and spreading across multiple GPUs. Imaging can be cleanly parallelized: chunks of a PSF and sample can be partially imaged on multiple GPUs independently because the convolution occurs per plane, then finally all partial images can be summed together onto a single GPU. The reconstructions can similarly take the final image and reconstruct partial chunks (as well as calculate losses on partial chunks) of the volume independently per device. We implicitly gather data to the same GPU when computing sums ( ∑ ) or means (E). The functions parallel image and compute PSF follow the definitions above in equations 14 and 12. In the algorithms shown here, parallel image applies the same convolution described above in equation 14.
Sparse gradients and downsampling We additionally support training and reconstructing only some of the planes for imaging and potentially a different partial group of planes during reconstruction, as a way to sparsely compute gradients for optimization of θm and save memory. The planes not imaged with gradients can still contribute to the image (without their gradients being tracked) in order to make the problem more difficult for the reconstruction network. Over multiple iterations, this can become equivalent to the more expensive densely computed gradient method, essentially trading training time for memory. An additional memory saving measure not written in the algorithms is to compute the PSF at a high resolution, then downsample the PSF using a 2D sum pool to preserve total energy in order to reduce memory usage when performing the imaging and reconstruction. We denote with no gradient tracking to show an operation without gradients.
A.3 IMPLEMENTATION DETAILS
Fourier convolution details Our Fourier convolution uses complex number weights, implemented as two channels of real numbers. Furthermore, in order to prevent the convolution from wrapping around the edges, we have to pad the input to double the size. The size of the weight must match the size of this padded input. This means that the number of parameters for our Fourier convolution implementation is 8× the number of parameters required for a global kernel in a spatial convolution (though the Fourier convolution is significantly faster). We do this to save an extra padding and Fourier operation, trading memory for speed. Because the simulation of imaging requires more memory than the reconstruction network, we found this to be an acceptable tradeoff.
Common network details All convolutions (including Fourier convolutions) use “same” padding. For FourierUNets and vanilla UNets, downsampling and upsampling is performed only in the x and y
Algorithm 1: Parallel PSF engineering by joint training of reconstruction network and phase mask. Microscope Mφ parameters are φ, reconstruction network Rθ parameters are θ, dataset is D, learning rates for φ and θ are αφ and αθ respectively, plane indices to image and reconstruct from zs, and weight for LNMSE is β. Input :Mφ,φ, αφ,Rθ,θ, αθ,D, zs, β
1 for v ∈ D do // select plane indices to be imaged with and without
gradients 2 zs,gradient, zs,no gradient ← select planes(zs)
// move sample planes to be imaged with gradients to multiple GPUs
3 vs,gradient ← scatter(v, zs,gradient) // move sample planes to be imaged without gradients to
multiple GPUs 4 vs,no gradient ← scatter(v, zs,no gradient) // compute PSF with gradients on multiple GPUs 5 ss,gradient ← parallel(compute PSF(Mφ, zjs) for zjs in zs,gradient) // compute partial image with gradients on multiple GPUs 6 cgradient ← parallel image(ss,gradient,vs,gradient) // compute PSF without gradients on multiple GPUs 7 with no gradient tracking 8 ss,no gradient ← parallel(compute PSF(Mφ, zjs) for zjs in zs,no gradient) 9 end
// compute partial image without gradients on multiple GPUs 10 with no gradient tracking 11 cno gradient ← parallel image(ss,no gradient,vs,no gradient) 12 end
// compute full image by summing partial images onto one GPU 13 c← ∑ [cgradient, cno gradient] // select plane indices to be reconstructed 14 zs,reconstruct ← select planes(zs)
// move sample planes that will be reconstructed to multiple GPUs
15 vs,reconstruct ← scatter(v, zs,reconstruct) // compute mean of high passed sample for loss normalization 16 µH(v) ← E(H(vs,reconstruct)2) // compute mean of sample for loss normalization 17 µv ← E(v2s,reconstruct) // move reconstruction networks to multiple GPUs 18 Rθ,s ← scatter(Rθ) // compute reconstruction and loss on multiple GPUs 19 L← parallel reconstruct/loss(c,vs,reconstruct,Rθ,s, µH(v), µv, β) // compute gradients for all parameters 20 gθ ← ∇θL 21 gφ ← ∇φL // update all parameters 22 θ←ADAM(αθ,θ, gθ) 23 φ←ADAM(αφ,φ, gφ) 24 end
dimensions (we do not downsample or upsample in z because there could potentially not be enough planes to do so). We train all networks using the ADAM optimizer with all default PyTorch parameters except the learning rate, which we always set to 10−4 for the reconstruction network parameters θ and 10−2 for the phase mask parameters φ.
Algorithm 2: Parallel training a reconstruction network given a pre-trained phase mask. Microscope Mφ parameters are φ (phase mask), reconstruction network Rθ parameters are θ, dataset is D, learning rates for φ and θ are αφ and αθ respectively, plane indices to image and reconstruct from are zs, and weight for LNMSE is β. Input :Mφ,φ, αφ,Rθ,θ, αθ,D, zs, β // compute PSF without gradients on multiple GPUs
1 with no gradient tracking 2 sno gradient ←parallel(compute PSF(Mφ, zjs) for zjs in zs) 3 end 4 for v ∈ D do
// select plane indices to be imaged without gradients 5 zs,no gradient ← select planes(zs)
// move sample planes to be imaged without gradients to multiple GPUs
6 vs,no gradient ← scatter(v, zs,no gradient) // move necessary PSF planes to multiple GPUs 7 ss,no gradient ← scatter(sno gradient, zs,no gradient) // compute image without gradients on multiple GPUs
8 with no gradient tracking 9 c← parallel image(ss,no gradient,vs,no gradient)
10 end // select plane indices to be reconstructed 11 zs,reconstruct ← select planes(zs) // move sample planes that will be reconstructed to multiple
GPUs 12 vs,reconstruct ← scatter(v, zs,reconstruct) // compute mean of high passed sample for loss normalization 13 µH(v) ← E[H(vs,reconstruct)2] // compute mean of sample for loss normalization 14 µv ← E[v2s,reconstruct] // move reconstruction networks to multiple GPUs 15 Rθ,s ← scatter(Rθ) // compute reconstruction and loss on multiple GPUs 16 L← parallel reconstruct/loss(c,vs,reconstruct,Rθ,s, µH(v), µv, β) // compute gradients for reconstruction networks only 17 gθ ← ∇θL // update reconstruction network parameters only 18 θ←ADAM(αθ,θ, gθ) 19 end
Algorithm 3: Parallel imaging. PSF planes on multiple GPUs are ss, sample planes on multiple GPUs to be imaged are vs. Input :ss,vs Output :c // compute images in parallel on multiple GPUs, then sum to
single GPU 1 c← ∑ [parallel(convolve(sjs,v j s) for (s j s,v j s) in (ss,vs))] 2 return c
Normalization We use input scaling during both training and inference in order to normalize out differences in the brightness of the image and prevent instabilities in our gradients. This means we divide out the median value of the input (scaled by some factor in order to bring the loss to a reasonable range) and then undo this scaling after the output of the network. This effectively linearizes our reconstruction networks, meaning a scaling of the image sent to the network will exactly scale the output by that value. We also find this is a more effective and simpler alternative to
Algorithm 4: Parallel reconstruction/loss calculation. Camera image is c, sample planes on multiple GPUs are vs, reconstruction networks on multiple GPUs are Rθ,s, mean for LHNMSE normalization is µH(v), mean for LNMSE normalization is µv, and weight for LNMSE is β.
Input :c,vs,Rθ,s, µH(v), µv, β Output :L // compute reconstruction and loss in parallel on multiple GPUs
1 v̂s ← concatenate(parallel(Rjs(c) forRjs inRθ,s)) 2 Ls ← parallel(E[(H(v j s)−H(v̂ j s))
2] µH(v) + β E[(vjs−v̂ j s) 2] µv for (v̂js,v j s) in (v̂s,vs))
// compute mean of scattered losses on single GPU 3 L← E[Ls] 4 return L
using a BatchNorm on our inputs. We continue to use BatchNorm between our convolution layers within the reconstruction network [15], which is effectively InstanceNorm in our case where batch size is 1 [31].
Planewise network training logic When we train PSFs by optimizing φ, we train separate reconstruction networks per plane. This allows us to flexibly compute sparse gradients across different planes from iteration to iteration, as described in Appendix A.2. In order to do this, we create placeholder networks on any number of GPUs, then copy the parameters stored on CPU for each plane’s reconstruction network to a network on the GPU as needed during a forward pass. After calculating an update with the optimizer, we copy the parameter values back to the corresponding parameter on CPU.
Training times We optimize our Type D microscopy experiments on 8 RTX 2080 Ti GPUs. For these, we can compare training times for the different network architectures. One training iteration (including microscope simulation, reconstruction, backpropagation, and parameter update) takes ∼0.6 seconds for FourierNet2D and ∼1.3 seconds for UNet2D when optimizing both φ and θ. One training iteration takes ∼0.4 seconds for FourierNet3D, ∼0.7 seconds for FourierUNet3D, and ∼0.8 seconds for UNet3D when only optimizing θ. More details are found in Tables 5 and 6.
Table 5: Type D experiment training times
Network Optimizing # parameters # train steps Train step time (s) Total time (h)
FourierNet2D θ,φ ∼ 4.2× 107 106 ∼ 0.8 ∼ 222 FourierNet3D θ ∼ 6.3× 107 106 ∼ 0.4 ∼ 111 FourierUNet3D θ ∼ 8.4× 107 106 ∼ 0.7 ∼ 194 UNet2D θ,φ ∼ 4.0× 107 106 ∼ 1.3 ∼ 361 UNet3D θ ∼ 1.0× 108 106 ∼ 0.8 ∼ 222
A.4 DETAILS FOR FOURIERNETS OUTPERFORM UNETS FOR ENGINEERING NON-LOCAL OPTICAL ENCODERS AND 3D SNAPSHOT MICROSCOPY VOLUME RECONSTRUCTION
For our experiments in Sections 3.1 and 3.2, we use 40 planes at 5µm resolution in z and therefore 40 reconstruction networks to train PSFs. When training reconstruction networks only to produce the higher quality reconstructions, we use 96 planes at 1µm resolution in z (chosen so that the planes actually span 200 µm in z). We train in both settings without any sparse planewise gradients, meaning we image and reconstruct all 40 or all 96 planes, respectively. We show details of all datasets used for training reconstructions in Table 4.
We show the details of our FourierNet2D architecture for training PSFs in Table 7 and our FourierNet3D architecture for training reconstruction networks in Table 8. We also show details for training times for both training PSFs and for training more powerful reconstruction networks in Table 5. We
Table 6: Type A, B, C experiment training times
Network Optimizing # parameters Type # train steps Train step time (s) Total time (h)
FourierNet2D θ,φ ∼ 1.7× 108 A 5.8× 105 ∼ 1.1 ∼ 177 FourierNet3D θ (fixed φ for A) ∼ 3.4× 108 A ∼ 2.6× 105 ∼ 1.6 ∼ 116 FourierNet3D θ (fixed φ for A) ∼ 3.4× 108 B ∼ 1.3× 105 ∼ 1.6 ∼ 58 FourierNet3D θ (fixed φ for A) ∼ 3.4× 108 C ∼ 1.3× 105 ∼ 1.6 ∼ 58 FourierNet2D θ,φ ∼ 1.7× 108 B 5.8× 105 ∼ 1.1 ∼ 177 FourierNet3D θ (fixed φ for B) ∼ 3.4× 108 A ∼ 1.2× 105 ∼ 1.6 ∼ 53 FourierNet3D θ (fixed φ for B) ∼ 3.4× 108 B 106 ∼ 1.6 ∼ 444 FourierNet3D θ (fixed φ for B) ∼ 3.4× 108 C ∼ 5.0× 105 ∼ 1.6 ∼ 222 FourierNet2D θ,φ ∼ 1.7× 108 C 5.8× 105 ∼ 1.1 ∼ 177 FourierNet3D θ (fixed φ for C) ∼ 3.4× 108 A ∼ 3.4× 105 ∼ 1.6 ∼ 151 FourierNet3D θ (fixed φ for C) ∼ 3.4× 108 B ∼ 3.4× 105 ∼ 1.6 ∼ 151 FourierNet3D θ (fixed φ for C) ∼ 3.4× 108 C ∼ 3.7× 105 ∼ 1.6 ∼ 164
trained all networks for Type D for the same number of iterations (more than necessary for PSFs to meaningfully converge)2.
The architecture of FourierUNet3D is 4 scales, with a cropping factor of 2 per scale in the encoding path and an upsampling factor of 2 in the decoding path. For each scale, we perform a Fourier convolution in the encoding path producing 480 feature maps, which are concatenated with the incoming feature maps of the decoding convolutions at the corresponding scale (just as in a normal UNet). In the decoding path, we use 3D convolutions with kernel size (3, 5, 5), producing 12 3D feature maps each. There are two such convolutions per scale. Note that this requires we reshape the 2D feature maps from the Fourier convolutions to 3D. This is followed by a 1x1 convolution
2Training times are approximate, and actual total time was longer due to checkpointing/snapshotting/validation of data and/or differences in load on the clusters being used.

heading: 2 (128, 128) (2, 2) (60, 1, 128, 128)

heading: 3 (64, 64) (2, 2) (60, 1, 64, 64)

heading: 4 (32, 32) (2, 2) (60, 1, 32, 32)
text: producing the 3D reconstruction output. We show a diagram of this architecture in Figure 1C, and details of this architecture in Table 9.
For our UNet2D, each encoding convolution produced 24 feature maps (except the first scale, for which the first convolution produced 12 feature maps and the second convolution produced 24 feature maps). Each decoding convolution produced 24 feature maps, but took an input of 48 feature maps where 24 feature maps were concatenated from the corresponding encoding convolution at that scale. At the end of the UNet2D, a (1, 1) convolution reduced the 24 final feature maps to 1 feature map. This single feature map is interpreted as the final output of the network, i.e. the reconstructed plane. UNet2D requires many more feature maps per plane and more layers than FourierNet, because these are necessary in order for the network to be able to integrate information from a larger field of view. We show the details of our UNet2D architecture in Table 10.
The architecture of the vanilla UNet3D is also 4 scales, with a max pooling factor of 2 per scale in the encoding path and an upsampling factor of 2 in the decoding path. Each scale of the encoding path produces 480 2D feature maps. These are concatenated to the incoming feature maps of the decoding convolutions at the corresponding scale, again with a reshape from 2D to 3D. Each scale of the decoding path produces 48 3D feature maps. Again, this is followed by a 1x1 convolution producing the 3D reconstruction output. All convolutions are in 3D with a kernel size of (5, 7, 7), with the z dimension being ignored for the encoding path because the input is 2D. We show the details of our UNet3D architecture in Table 11.
A.5 DETAILS FOR ENGINEERED OPTICAL ENCODING DEPENDS ON SAMPLE STRUCTURE
For our experiments in Section 3.3, we use 64 planes at 1µm resolution in z and therefore 64 reconstruction networks to train PSFs. When training reconstruction networks only to produce the higher quality reconstructions, we use 128 planes at 1µm resolution in z (chosen so that the planes actually span 250 µm in z). We train in the reconstruction only setting without any sparse planewise gradients, meaning we image and reconstruct all 128 planes. However, when training a PSF we image and reconstruct 40 planes at a time with gradient per iteration (spread across 8 GPUs). These 40
planes are chosen randomly at every iteration from the 64 total possible planes, making potentially separate draws of planes for imaging and reconstruction. We show details of all datasets used for training reconstructions in Table 4.
We show the details of our FourierNet2D architecture for training PSFs at the larger field of view in Type A, B, C in Table 12 and our FourierNet3D architecture for training reconstruction networks at the larger field of view in Type A, B, C in Table 13. There are no other networks used for these larger field of view experiments. We also show details for training times for both training PSFs and for training more powerful reconstruction networks in Table 6. All PSFs in these networks were trained for the same number of iterations. However, reconstruction networks for some of these experiments were only trained for as long as necessary to converge (with some exceptions where we attempted longer training to check for performance gains with long training periods). Generally, we observed that performance for such reconstruction networks does not meaningfully change with many more iterations of training3.
3Training times are approximate, and actual total time was longer due to checkpointing/snapshotting/validation of data and/or differences in load on the clusters being used.
A.6 DETAILS FOR FOURIERNETS OUTPERFORM STATE-OF-THE-ART FOR RECONSTRUCTING NATURAL IMAGES CAPTURED BY DIFFUSERCAM LENSLESS CAMERA
We performed no augmentations for this set of trainings reconstructing RGB color images of natural scenes from RGB diffused images taken through a DiffuserCam [20]. We modified our FourierNet2D architecture to create the FourierNetRGB architecture and our FourierUNet2D architecture to create
Table 13: FourierNet3D detailed architecture (8 GPUs)
Layer type Kernel size Stride Notes Shape (C, D, H, W) InputScaling - - scale: 0.01 (1, 1, 512, 512) FourierConv2D (512, 512) (2, 2) - (80, 1, 512, 512) LeakyReLU - - slope: -0.01 (80, 1, 512, 512) BatchNorm2D - - - (80, 1, 512, 512) Reshape2D3D - - - (5, 16, 512, 512) Conv3D (11, 7, 7) (1, 1, 1) - (5, 16, 512, 512) LeakyReLU - - slope: -0.01 (5, 16, 512, 512) BatchNorm3D - - - (5, 16, 512, 512) Conv3D (11, 7, 7) (1, 1, 1) - (1, 16, 512, 512) ReLU - - - (1, 16, 512, 512) InputRescaling - - scale: 0.01 (1, 16, 512, 512)
Table 14: DLMD experiment training times. Superscripts denote loss function: 1 MSE, 2 MSE+LPIPS.
Network Optimizing # parameters # train steps Train step time (s) Total time (h)
FourierNetRGB1 θ ∼ 1.6× 107 2.2× 105 ∼ 0.43 ∼ 26 FourierNetRGB2 θ ∼ 1.6× 107 1.1× 105 ∼ 0.47 ∼ 14 FourierUNetRGB1 θ ∼ 7.1× 107 2.5× 105 ∼ 3.3 ∼ 229 Le-ADMM-U2 [20] θ ∼ 4.0× 107 - - - UNet2 [20] θ ∼ 1.0× 108 - - -
the FourierUNetRGB architecture, outlined in Table 15 and Table 16 respectively. Training details are shown in Table 14. Because these reconstructions are of 2D images only and required no microscope simulation, we were able to use a batch size of 4 images per iteration.

Review summaries:
This paper introduces an end-to-end optical encoder and deep learning decoder optimization for 3D snapshot microscopy. The main challenge is that the decoder needs to be able to handle the global PSFs used in the optical encoder in 3D snapshot microscope. This is difficult to be achieved in conventional UNet architectures as the global context requires many layers that is computationally expensive. The authors use a Fourier convolutional network that achieves the global context in a single layer with a much lower computational cost. The method is demonstrated on simulated data of one volume with four different fields of view. The authors also applied the technique to lensless imaging. 

This paper introduces the use of Fourier convolutional neural networks to model both image encoding and reconstruction processes. These neural networks learn multiplicative weights in the Fourier domain, enabling efficient global convolutions (i.e., convolutions with kernels the size of the entire image), which are particularly important for imaging applications with a wide PSF. These networks compare favorably to UNets for the task of joint 3D snapshot microscope parameter optimization and image reconstruction. These networks also outperform the state-of-the-art architectures for a lensless imaging reconstruction task.

This paper introduces a convolutional neural network architecture for designing phase patterns to display on a spatial light modulator in the imaging path of a microscope, as well as decoding the resulting image into a 3D volume. The authors propose to use Fourier features in the decoding stage to enable features that are more strongly nonlocal than in regular (primal-domain) convolution than what they claim would be possible to achieve using "vanilla" UNet/ResNet-like architectures. They demonstrate their approach for snapshot 3D microscopy as well as the reconstruction of light fields from input captured using lensless system (with randomized optical elements). 

This paper proposes a neural network architecture that applies a series of convolution kernels on the Fourier domain to focus on global information of the captured images. The authors make use of the Fourier-plane phase modulation and apply this idea to single-shot 3D microscopic imaging using the Fourier-plane phase modulation for 3D microscope and lensless imaging. Also, this is further extended to the end-to-end optimization of the phase modulation and the reconstruction network. An extensive synthetic evaluation is conducted to validate the effectiveness of the proposed method over existing non-Fourier neural networks in the context of computational imaging.