# High-resolution limited-angle phase tomography of dense layered objects using deep neural networks

^{a}3D Optics Laboratory, Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139;^{b}Microsystems Technology Laboratories, Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139;^{c}BioSystems and bioMechanics (BioSyM) Interdisciplinary Research Group, Singapore-MIT Alliance for Research and Technology (SMART) Centre, Singapore 117543, Singapore

See allHide authors and affiliations

Edited by Terrence J. Sejnowski, Salk Institute for Biological Studies, La Jolla, CA, and approved August 19, 2019 (received for review December 14, 2018)

## Significance

We demonstrate that it is possible to use deep neural networks to produce tomographic reconstructions of dense layered objects with small illumination angle as low as 10 °. It is also shown that a DNN trained on synthetic data can generalize well to and produce reconstructions from experimental measurements. This work has application in the field of X-ray tomography for the inspection of integrated circuits and other materials studies.

## Abstract

We present a machine learning-based method for tomographic reconstruction of dense layered objects, with range of projection angles limited to

Tomography is the quintessential inverse problem. Since the interior of a 3-dimensional (3D) object is not accessible noninvasively, the original insight of tomographic approaches was to illuminate through from multiple angles of incidence and then process the resulting projections to reconstruct the interior slice by slice (1⇓–3). In the simplest case, when diffraction is negligible and the illumination is collimated, as is generally permissible to assume for X-ray attenuation (4⇓⇓–7) and electron scattering (8⇓–10) in the far field and for features of size ∼1 μm and above, the object’s interior is represented by its Radon transform (11) of line integrals along straight parallel paths. The interior of the volume is then reconstructed by use of the Fourier-slice theorem for the Radon projections. On the other hand, if the X-ray beam is not collimated but spherical, then the slice-by-slice approach is no longer applicable and full volumetric reconstruction is required (12, 13). Even when the object is available for observation from the full

Additional challenges occur when the inverse problem of interest is to reconstruct in 3D the index of refraction, rather than the attenuation. If the object features are large enough compared to the wavelength, such that diffraction may still be neglected, and the index variations through the object volume are relatively small, then each projection may be modeled as a set of Fermat integrals of phase delay along approximately straight lines. The phase integrals may be obtained, for example, using holographic interferometry (15, 16) or transport of intensity (17). For smaller-sized features and still assuming weak scattering (first-order Born approximation), the projection integrals are instead obtained along curved paths on the surface of the Ewald sphere, a method referred to as diffraction tomography (18, 19). By decoupling the problem into 2 parts, first phase projection retrieval, followed by tomography, these approaches enjoy the benefit of using the advanced algorithms in the 2 respective research fields. However, there is also the danger that errors generated independently during each step may amplify each other. Finally, when strong scattering may no longer be neglected, all 2-step approaches become questionable because the interpretation of the first step as line integrals is no longer valid.

Generally, ill-posed inverse problems are solved by regularized optimization. If f is the object and g the measurement, then the object estimate **1** if a set of basis functions where the object class is sparse is a priori known. Alternatively, if a database of representative objects is available, then these examples may be used to learn the optimal set of basis functions as a dictionary (28, 29).

Rapid recent developments in the field of machine learning, and deep neural networks (30) (DNNs) in particular, have provided an additional set of tools and insights for inverse problems. It may be shown (31, 32) that recurrent or unfolded multistage DNN architectures are formally equivalent to the iterative solution to the inverse problem in Eq. **1** where the prior Φ need no longer be known or depend on sparsity; instead, examples guide the discovery of the prior through the DNN training process. Simpler learning architectures, where g is fed to the DNN directly or after first passing through a preprocessor, have been used for retrieval of phase from intensity (33⇓⇓–36); 3D holographic reconstruction (37⇓–39); superresolution photography (40⇓–42) and microscopy (43); imaging through scatter (44⇓⇓–47); and imaging under extremely low light conditions in the 3 contexts computational ghost imaging (48), consumer-camera photography (49), and phase retrieval (50).

Multistage DNN architectures have been shown to yield high-quality reconstructions in numerous Radon tomography configurations (32, 51⇓⇓⇓⇓–56). Recently, Nguyen et al. (57) used the inverse Radon transform for optical tomography with a single-stage DNN intended to partially correct for the assumption of line integrals breaking down.

In this paper, we apply a Fourier-based beam propagation method (BPM) (58) as a preprocessing step immune from any Radon assumptions. The strongly scattering object is illuminated by a parallel beam under a limited angular range of

Large datasets, typically consisting of more than 5,000 examples, are generally required for DNN training. That is feasible in many cases through spatial light modulators (33, 57). However, in many cases of interest spatial light modulators have insufficient space-bandwidth product or are unavailable (e.g., in X-rays); and alternatives to generate physical specimens are expensive or restricted due to proprietary processes. Instead, our approach is to train the DNN on purely synthetic data with the rigorous BPM forward model and then use a physical test specimen (phantom) to test the reconstruction quality with well-calibrated ground truth in experiments.

We chose to design our phantom as emulating the 3D geometry of integrated circuits (ICs). These would normally be inspected with X-rays, so we scaled up the feature dimensions in the phantom for visible wavelengths. The advantage of this choice is that IC layouts provide strong geometrical priors, e.g., Manhattan geometries, and our phantom also exhibited large spatial gradients and refractive index contrast to strengthen scattering. Thus, our methodology is directly applicable to all cases of tomography at optical wavelengths, e.g., 3D-printed specimen characterization and identification and biological studies in cells and tissue with moderate scattering properties. In each case, testing the ground truth would require the fabrication or the accurate simulation of different phantoms meeting the corresponding priors.

There is also value in the study of emulating X-ray inspection of ICs at visible wavelengths, as extensive outsourcing by the IC industry has created a growing concern that the ICs delivered to the customer may not match the expected design and that malicious features may have been added (59). However, in our emulation the phase contrast of the features against the background and the Fresnel number are both higher than typical corresponding IC configurations even at soft X-ray wavelengths.

One advantage of deep learning for inverse problems is speed. Solving **[1]** separately for each instance of g is computationally intensive, and training a DNN is even more so. This is because both operations are iterative, and the latter is run on large datasets. On the other hand, once the DNN has learned the inverse map from the preprocessed version of g to **1** with the forward operator H itself consisting of an expensive computational procedure in each iteration. In our approach, the preprocessing performed prior to the DNN is the most time-consuming operation, and therefore we aim at simplifying the preprocessing step as much as we can, i.e., tolerate a crude approximation, and leave it to the DNN to correct it. In our case, the execution time is 51 s (out of which only 300 ms are taken by the DNN, with the rest being the preprocessor), while learning tomography (60, 61), which is based on a similar gradient descent algorithm, takes 212 s and yields inferior reconstructions (see Fig. 4 for this, and see *Materials and Methods* for the computing hardware implementation).

## Optical Experiment

We prepared a series of 4 glass wafers with etched structures representing patterned layers from an actual IC design. A schematic cross-section of the sample is shown in Fig. 1*A*. The glass plates are held together and aligned in a custom-made holder. Immersion oil is added between the plates to minimize parasitic reflections and also tune the phase shift associated with the patterns. The pattern depth was measured to be 575 nm, yielding a phase shift of −0.33 rad for the particular oil used. Note that the phase shift is negative as the refractive index of the oil is lower than the refractive index of the glass. Details about the sample preparation and phase-shift measurements are given in *Materials and Methods*. The particular patterns etched on the sample are shown in Fig. 1*B*.

The experimental apparatus is detailed in Fig. 2. A collimated monochromatic plane wave from a continuous wave (CW) laser is incident on the sample, which is mounted on a 2-axis rotation stage. The sample is imaged through a demagnifying telescope to increase the field of view. The detector (an EM-CCD camera) is defocused from the image plane to simulate free-space propagation in an X-ray experiment where no imaging system can be used. Further details are given in *Materials and Methods*.

The strength of diffraction effects can be quantified with the Fresnel number

## Computation of the Approximant

As mentioned above, the task of the DNN is significantly facilitated if the raw measurements are preprocessed to give an approximation of the solution. We use a simple gradient descent method to generate the approximant of the sample refractive index distribution.

Light propagation through the object can be computed using the following split-step Fourier BPM. In this method each sample layer l is modeled as a 2D complex mask

From the measurements, we produce an approximation **2** written in operator form**4** and **5** is given after the conclusion in the *Derivation of the Gradient* section. In Fig. 3, we give an example of 1 experimental intensity measurement (**6**.

## DNN Architecture and Training

We use a DNN to map the approximant to the final reconstruction **2**. The synthetic measurements were subject to simulated shot noise and read noise equivalent to the noise levels found in the experiment. The shot noise was accounted for by converting the simulated measurement pixel intensities I expressed in average photon count per pixel per integration period of the detector to integer numbers of photons N following Poisson statistics. The actual optical power on the camera was measured with a power meter and converted to an average photon flux per detector pixel. Read noise following Gaussian statistics was subsequently added. The parameters (variance and average) of the noise were measured from a series of dark frames from the camera taken with the same gain (EM gain of 1) and integration time (2 ms) as the experimental measurements.

From each set of measurements, we produce a multilayer approximant using the gradient descent in Eq. **6**. The examples are split in a training set of 5,000 examples, a validation set of 450 examples, and a test set of 50 examples. Each set of measurements (1 example) comprises 22 views corresponding to different orientations of the sample. The DNN is then trained to map the approximant to the ground truth used for the simulation. Each layer of the sample is assigned to a different channel in the DNN. We use the negative Pearson correlation coefficient (

## Results

The method described in the previous sections was applied to the glass phantom shown in Fig. 1*B*. The synthetic measurements were subject to Poisson noise resulting from *E*–*H* for *I*–*L* for *M*–*P* (*Q*–*T* (*I*–*L*, the IC patterns (where the phase shift actually occurs) are typically reconstructed with 0 phase due to the rectified linear units (which project all negative values to 0) at the output layer of the DNN. We reassign the phase of the pattern to the nominal phase of −0.33 rad so that it can be visually compared to the ground truths in Fig. 4 *M*–*P*. An alternate approach leading to very similar results is to assign a 0 phase to the background.

The DNN reconstructions can be compared to those obtained using learning tomography (LT), a previously demonstrated optical tomography technique (60, 61) based on proximal optimization (FISTA) (63) with TV regularization (64). The role of the TV regularizer is to favor piecewise constant solutions while preserving sharp edges, which is especially well suited for IC patterns. LT was initially designed for holographic measurements and was modified here to work on intensity measurements by computing the gradient for the data fidelity term in Eq. **3**. The essential difference in the LT optimization is that a TV filter playing the role of a proximal operator is applied at each iteration on the current solution. The LT reconstructions for the experimental dataset are shown in Fig. 4 *A*–*D*. These particular reconstructions were obtained after 30 iterations of gradient descent, a step size of 0.05, a regularization parameter *Materials and Methods* for hardware details).

In Table 1, we summarize the values of the PCC, which we use to quantify the quality of the reconstructions. The values are given for reconstructions on the synthetic test set (50 examples) and also the reconstruction of the single experimental example. Because the reconstruction quality turns out to depend strongly on the particular layer, we display the value for each of the 4 layers separately. As may be expected, the values for the LT are higher (better reconstructions) than those for the approximant as LT was run for 30 iterations vs. 8 for the approximant and that the latter was not regularized. The DNN reconstructions appear to be the best according to the PCC metric, which shows that, even while starting from a poor approximation, the DNN was able to outperform LT. Note that a direct comparison between the performance of the DNN on the synthetic data and that on the experimental example is not fair because the ground truth is not known in the experiment. The ground truth used for the experimental example is an idealization from the design parameters used to fabricate the sample. We also indicate the values for reconstructions based on the regularized approximant (using the same regularization parameters as in the LT algorithm). In terms of PCC, there is no significant difference from the unregularized case.

The reconstructions based on regularized approximants are shown in Fig. 5. By comparing these images with the unregularized reconstructions shown in Fig. 4, and also by considering the value of the PCC in Table 1, we conclude that the regularization has little effect, especially on the experimental reconstructions. Moreover, the TV regularization may not operate as a favorable preconditioner for the DNN. The choice of the TV operator as a regularizer is arbitrary and only based on our assumption that the solution should be piecewise constant. In fact, because of the small angle range, the approximants for the different phantom layers are quite similar to each other and the regularization may cancel information that the DNN could use to discriminate between them. Layers 3 and 4 can be said to look visually better with the regularized approximant, but the situation is reversed for layers 1 and 2. Iterating more, i.e., using

In the regularized case *I*–*L*). For the regularized

So far, we have reconstructed the phase shift distribution

## Conclusion

In this paper, we have demonstrated through an emulated X-ray experiment that DNNs can be used to improve the reconstruction of IC layouts from tomographic intensity measurements acquired within an angle range of

One significant challenge for the method we demonstrated is to provide a proper training set for the DNN. In the case of ICs, the training set is simply given by the many real layouts that are available. For more generic objects, a problem needs first to be formulated to clearly define the class of object on which the training will be performed. This is, however, the case in general for problem solving that involves DNNs. What has been shown here is the compelling improvement that DNN can bring to a phase tomography problem when the class of object is known.

More generally, there is a trade-off between the specificity of the required priors and the “difficulty” of an inverse problem—measured as degree of ill-posedness, e.g., the ratio of largest to smallest eigenvalue in a linearized representation of the forward operator. The problem we addressed here is severely ill-posed due to the presence of strong scattering within the object and the limited range and number of angular measurements we allowed ourselves to collect. Therefore, the rather restricted nature of ICs as a class prior is justified; while, at the same time, our approach is addressing an indisputably important practical problem. Detailed determination of the relationship between the degree of ill-posedness and the complexity of the object class prior would be a worthwhile topic for future work.

### Derivation of the Gradient.

This derivation follows a path similar to the derivation given in ref. 61. We start from Eq. **3**:**9** and, by linearity of the derivation operation and denoting **6**:**15**, **4** and **5** that describe the forward model where we drop the index i to simplify the notation as the expression assumes the same form for each sample orientation:**16**,**17** as a recursive relationship for the optical field **23**–**25**:**24** and **25**,**35** which reads, for pure phase objects (**24**. Finally, according to Eq. **22**, we obtain layer l of

## Materials and Methods

The experimental apparatus is shown in Fig. 2*A*. The light source is a CW He-Ne laser at *x* and *y* axes. The sample middle plane is imaged using a demagnifying telescope (

The sample corresponds to a *C* that contain copper patterns that would induce a significant phase delay in the X-ray regime. The 4 glass plates corresponding to the IC layers were cut in a 500-μm-thick fused silica wafer and *B*) with a refractive index

The sample layers are fabricated on double-sided polished 150-mm-diameter and 500-μm-thick fused silica wafers. A 1-μm-thick positive tone resist (Megaposit SPR700) is spin coated at 3,500 rpm on both sides of the wafer and soft baked at 95 °C for 30 min in a convection oven. The backside was also coated for protection from the forthcoming wet etch. Scaled versions of IC designs in GDSII format are then patterned directly using a maskless aligner (MLA150; Heidelberg Instruments) with a 405-nm laser and developed using an alkaline developer (Shipley Microposit MF CD-26) for 45 s followed by a deionized (DI) water rinse and

The computation of the approximant is performed with the MATLAB software on an Intel i9-7900X processor running at 3.3 GHz. The DNN training and testing are performed under Keras with Tensorflow backend running on an NVIDIA Titan Xp graphics processing unit.

## Acknowledgments

This work was supported by the Intelligence Advanced Research Projects Activity (FA8650-17-C-9113).

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: agoy{at}goyman.com.

Author contributions: A.I.A. and G.B. designed research; A.G. performed research; G.R., S.L., and K.A. contributed new reagents/analytic tools; A.G. analyzed data; and A.G., G.R., and G.B. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

- Copyright © 2019 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- ↵
- ↵
- R. N. Bracewell

- ↵
- ↵
- ↵
- G. N. Hounsfield

- ↵
- J. Ambrose

- ↵
- D. J. de Rosier,
- A. Klug

- ↵
- D. J. de Rosier,
- A. Klug

- ↵
- ↵
- J. Radon

- ↵
- B. K. P. Horn

- ↵
- B. K. P. Horn

- ↵
- ↵
- D. W. Sweeney,
- C. M. Vest

- ↵
- C. M. Vest

- ↵
- J. C. Petruccelli,
- L. Tian,
- G. Barbastathis

- ↵
- ↵
- ↵
- A. N. Tikhonov

- ↵
- A. N. Tikhonov

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- D. J. Brady,
- A. Mrozack,
- K. MacCabe,
- P. Llull

- ↵
- M. Elad,
- M. Aharon

- ↵
- ↵
- ↵
- K. Gregor,
- Y. LeCun

- ↵
- M. Mardani et al.

- ↵
- ↵
- ↵
- S. Li,
- G. Barbastathis

- ↵
- Z. D. C. Kemp

- ↵
- T. Shimobaba,
- T. Kakue,
- T. Ito

- ↵
- Z. Ren,
- Z. Xu,
- E. Y. Lam

- ↵
- Y. Wu et al.

- ↵
- C. Dong,
- C. Loy,
- K. He,
- X. Tang

- ↵
- C. Dong,
- C. Loy,
- K. He,
- X. Tang

- ↵
- J. Johnson,
- A. Alahi,
- L. Fei-Fei

- ↵
- ↵
- M. Lyu,
- H. Wang,
- G. Li,
- G. Situ

- ↵
- S. Li,
- M. Deng,
- J. Lee,
- A. Sinha,
- G. Barbastathis

- ↵
- N. Borhani,
- E. Kakkava,
- C. Moser,
- D. Psaltis

- ↵
- Y. Li,
- Y. Xue,
- L. Tian

- ↵
- M. Lyu et al.

- ↵
- C. Chen,
- Q. Chen,
- J. Xu,
- V. Koltun

- ↵
- A. Goy,
- K. Arthur,
- S. Li,
- G. Barbastathis

- ↵
- M. Mardani et al.

- ↵
- J. Schlemper,
- J. Caballero,
- J. V. Hajnal,
- A. N. Price,
- D. Rueckert

- ↵
- ↵
- M. T. McCann,
- K. Hwan Jin,
- M. Unser

- ↵
- H. Gupta,
- K. Hwan Jin,
- H. Q. Nguyen,
- M. T. McCann,
- M. Unser

- ↵
- Y. Sun,
- Z. Xia,
- U. S. Kamilov

- ↵
- T. Nguyen,
- V. Bui,
- G. Nehmetallah

- ↵
- M. D. Feit,
- J. A. Fleck Jr

- ↵
- M. A. Mak

- ↵
- ↵
- U. S. Kamilov et al.

- ↵
- G. Huang,
- Z. Liu,
- K. Q. Weinberger

- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Engineering