A Learning Approach to Optical Tomography

We describe a method for imaging 3D objects in a tomographic configuration implemented by training an artificial neural network to reproduce the complex amplitude of the experimentally measured scattered light. The network is designed such that the voxel values of the refractive index of the 3D object are the variables that are adapted during the training process. We demonstrate the method experimentally by forming images of the 3D refractive index distribution of cells.

We describe a method for imaging 3D objects in a tomographic configuration implemented by training an artificial neural network to reproduce the complex amplitude of the experimentally measured scattered light. The network is designed such that the voxel values of the refractive index of the 3D object are the variables that are adapted during the training process. We demonstrate the method experimentally by forming images of the 3D refractive index distribution of cells.
The learning approach to imaging we describe in this paper is related to adaptive techniques in phased antenna arrays 1 and inverse scattering 2, 3 . In the optical domain an iterative approach was demonstrated by the Sentenac group 4, 5 who used the coupled dipole approximation 6 for modelling light propagation in inhomogeneous media (a very accurate method but computationally intensive) to simulate light scattering from small objects (1µm × 0.5µm) in a point scanning microscope configuration. Very recently an iterative optimization method was demonstrated 7 for imaging 3D objects using incoherent illumination. Our method relies on digital holography 8,9 to record the complex amplitude of the field. We use the Beam Propagation Method (BPM) 10,11 to model the 1 arXiv:1502.01914v1 [physics.optics] 5 Feb 2015 scattering process and the error back propagation method 12 to train the system. At the end of the training process the network discovers a 3D index distribution that is consistent with the experimental observations. We experimentally demonstrate the technique by imaging polystyrene beads and HeLa and hTERT-RPE1 cells.

Experimental Setup
A schematic diagram of the experimental setup is shown in Figure 1. It is a holographic tomography system 13 , in which the sample is illuminated with multiple angles and the scattered light is holographically recorded. Several variation of the holographic tomography system have been demonstrated before [14][15][16][17] . The optical arrangement we used is most similar to the one described by Choi et al. 14 . The samples to be measured were prepared by placing polystyrene beads and cells between two glass cover slides. The samples were illuminated with a continuous collimated wave at 561nm at 80 different angles. The amplitude and phase of the light transmitted through the sample was imaged onto a 2D detector where it was holographically recorded by introducing a reference beam. These recordings constitute the training set with which we train the computational model that simulates the experimental setup. We construct the network using the BPM. The inhomogeneous medium (beads or cells) is divided into thin slices along the propagation direction (z). The propagation through each slice is calculated as a phase modulation due to the local transverse index variation followed by propagation in a thin slice of a homogenous medium having the average value of the index of refraction of the sample.

Methodology
A schematic description of the BPM simulation is shown in Figure 2. The straight lines connecting any two circles represent multiplication of the output of the unit located in the l-th layer of the network at x = n 1 δ, y = m 1 δ by the discretized Fresnel diffraction kernel e jπ[(n 2 where n l and m l are integers and λ is the wavelength of light. δ is the sampling interval in the transverse coordinates (x, y) whereas δ z is the sampling interval along the propagation direction z. The circles in the diagram of Figure 2 perform a summation of the complex amplitude of the signals converging to each circle and also multiplication of this sum by e j(2π∆nδzz)/λ . ∆n(x, y, z) is the unknown 3D index perturbation of the object.
In the experiments the network has 420 layers with ∆n(x, y, z) being the adaptable variable.
In contrast with a conventional neural network, the output of the layered structure in Figure 2 is a linear function of the input complex field amplitude. However, the dependence of the output is nonlinearly related to ∆n(x, y, z). The BPM can be trained using steepest descent exactly as the back propagation algorithm in neural networks [18][19][20] . Specifically, the learning algorithm carries out the following minimization: In the above expression E k (∆n) is the current prediction of the BPM network for the output when the system is illuminated with the k-th beam and M k (∆n) is the actual measurement obtained by the optical system. ∆n indicates the estimate for the index perturbation due to the object. The term S(∆n) is a sparsity constraint 21 to enhance the contrast while τ is a parameter that can be tuned to maximize image quality by systematic search. The positivity constraint takes advantage of the assumption that the index perturbation is real and positive. The optimization is carried out iteratively by taking the derivative of the error with respect to each of the adaptable parameters following steepest descent is the error, α is a constant and the change in ∆n is proportional to the error and its derivative. This is achieved efficiently via a recursive computation of the gradient, which is the back propagation part of our learning algorithm.

Results
We first tested the system with polystyrene beads encapsulated between two glass slides in immersion oil. The sample was inserted in the optical system of Figure 1 and 80 holograms were recorded by illuminating the sample at 80 distinct angles uniformly distributed in the range -45 degrees to +45 degrees. The collected data is the training set for the 420-layer BPM network which simulates a physical propagation distance of 30µm and transverse window 37µm × 37µm (δ x = δ y = 72nm). The network was initialized with the standard filtered back projection reconstruction algorithm (Radon transform) 22 and the resulting 3D images before and after 100 iterations are shown in Figure 3. The final image produced by the learning algorithm is an accurate repro- duction of the bead shape.
A sample of a HeLa cell was also prepared and the same procedure was followed to obtain a 3D image. The results are shown in Figure 4 where the error function is plotted as a function of iteration number. In this instance, the system was initialized with a constant but nonzero value (∆n = 0.007). Also shown in Figure 4 are the results obtained when the system was initialized with the Radon reconstruction from the same data. After 100 iterations both runs yield essentially identical results. Notice that the error in the final image (after 100 iterations) is significantly lower than to the error of the Radon reconstruction. This is also evident by visual inspection of the 6 images in Figure 4 where the artifacts due to the missing cone 23 and diffraction 14 are removed by the learning process.
In general, optical 3D imaging techniques rely on the assumption that the object being imaged does not significantly distort the illuminating beam. This is assumed for example in Radon or diffraction/holographic tomography. In other words, these 3D reconstruction methods rely on the assumption that the measured scattered light consists of photons that have only been scattered once before they reach the detector. The BPM, on the other hand, allows for multiple forward scattering events. The only simplification is that reflections are not taken into account; these could eventually be incorporated in the network equation without fundamentally altering the approach described in this paper. Since biological tissue is generally forward scattering, BPM can be a good candidate to model propagation of thick biological samples and this may be the most significant advantage of the learning approach. To demonstrate this point, we prepared two glass slides with a random distribution of hTERT-RPE1 cells (immortalized epithelial cells from retina) on each slide. When we attach the two slides together, we can find locations where two cells are aligned in z, one on top of the other. Figure 5 (a)-(e) shows the image of such a stack of two cells produced with a direct inversion using the Radon transform. Figure 5 (f)-(j) shows the same object imaged with the proposed learning algorithm. The learning method was able to distinguish the two cells where the Radon reconstruction merged the two into a single pattern due to the blurring in z which is a consequence of the missing cone.
In conclusion, we have demonstrated a neural-network based algorithm to solve the optical