Research Article

Physics-informed neural networks for diffraction

tomography
Amirhossein Saba ,* Carlo Gigli ,,† Ahmed B. Ayoub,,† and Demetri Psaltis
École Polytechnique Fédérale de Lausanne, Optics Laboratory, Lausanne, Switzerland

Abstract. We propose a physics-informed neural network (PINN) as the forward model for tomographic
reconstructions of biological samples. We demonstrate that by training this network with the Helmholtz
equation as a physical loss, we can predict the scattered field accurately. It will be shown that a pretrained
network can be fine-tuned for different samples and used for solving the scattering problem much faster than
other numerical solutions. We evaluate our methodology with numerical and experimental results. Our PINNs
can be generalized for any forward and inverse scattering problem.
Keywords: deep learning; physics-informed neural networks; scattering; three-dimensional imaging; optical diffraction
tomography.
Received Jul. 9, 2022; revised manuscript received Sep. 7, 2022; accepted for publication Oct. 27, 2022; published online Nov.
22, 2022.
© The Authors. Published by SPIE and CLP under a Creative Commons Attribution 4.0 International License. Distribution or
reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
[DOI: 10.1117/1.AP.4.6.066001]

1 Introduction experimental measurements of the projection as a loss function,

and update the estimation of the refractive index upon the min-
Optical diffraction tomography (ODT) is an imaging technique
imization of this loss function plus possibly any other prior
for extracting the three-dimensional (3D) refractive index distri-
knowledge, e.g., sparsity conditions. Importantly, such an iter-
bution of a sample, e.g., a biological cell using multiple two-
ative scheme requires an analytical/semi-analytical model to
dimensional (2D) images acquired at different illumination
angles. The refractive index of the sample provides useful mor- backpropagate the computed loss and update the estimation
phological information, making ODT an interesting approach for of the refractive index. This precludes the use of common
biological applications.1–3 The conventional method to reconstruct mesh-based numerical solvers such as finite difference and finite
the 3D refractive index from multiple projections was proposed element methods (FEMs). In Ref. 5, the authors use a linear
by Wolf in 1969.4 In Wolf’s method, the 3D Fourier domain of (single-scattering) forward model, in the approach proposed
the refractive index is filled with 2D Fourier transforms of the in Refs. 6–8, referred to as learning tomography. The forward
measured scattered fields. However, due to the finite number model is the beam propagation method (BPM) and, in Ref. 9,
of projections, the limited numerical aperture of the optical sys- the authors resort to the Lippmann–Schwinger equation. The
tem, and the single scattering approximation, missing frequencies forward models used in these iterative solutions either have in-
are causing some distortion, elongation, and underestimation of accuracies in the cases of multiple-scattering and high-contrast
the reconstructed refractive index.5 samples, or they are computationally demanding. As a result,
In the last several years, many different iterative methods presenting a fast, accurate, and differentiable forward model
have been proposed to reconstruct accurate refractive indices is necessary to be used in iterative ODT. Physics-informed
from ill-posed measurements.5–9 The main idea behind these neural networks (PINNs) can be a good candidate for solving
iterative approaches is to use a forward model that predicts forward scattering problems and for being used in iterative
the candidate projections for the current estimation of the refrac- tomographic reconstruction.
tive index in that iteration, compare this 2D prediction with the PINNs have recently gotten intense research attention for
solving different complex problems in physics.10,11 These net-
*Address all correspondence to Amirhossein Saba, amirhossein.sabashirvan@ works use physics laws as the loss function instead of the
epfl.ch data-driven loss functions. In conventional supervised deep
†
These authors contributed equally to this work. learning, a large dataset of labeled examples is used for the

Advanced Photonics 066001-1 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

training process: by comparing the known ground truth with the In this paper, we extend this idea for inhomogeneous and 3D
predictions from a deep multi-layer neural network, one can cases and present a MaxwellNet that is able to solve different
construct a loss function and tune the parameters of the network forward scattering problems, such as light scattering from bio-
to solve complex physical problems. Different examples of logical cells. In the first part of the work, we train MaxwellNet
these data-driven neural networks are proposed for optical ap- for 2D digital phantoms and show how this pretrained network
plications such as resolution enhancement,12 imaging through can be fine-tuned to predict light scattering from more complex
multi-mode fibers,13,14 phase retrieval,15 ODT,16 and digital and experimentally relevant samples, in our case, HCT-116
holography.17,18 In these networks, the knowledge acquired by cells. We benchmark the performance of MaxwellNet in solving
the network strongly depends on the statistical information pro- scattering problems for 2D and 3D objects. Next, we demon-
vided in the dataset, and training such a network requires access strate that such PINN can be efficiently used to invert the
to a large dataset. In contrast, PINNs directly minimize the scattering problem through an iterative scheme and improve
physical residual from the corresponding partial differential the results of conventional ODT. We first demonstrate the
equation (PDE) that governs the problem instead of extrapolat- reconstruction of the refractive index distribution from synthetic
ing physical laws after going through a large amount of exam- data and then we validate the technique with experimental mea-
ples. In the pioneering approach proposed by Lagaris et al.,19 the surements of scattering from polystyrene microspheres.
neural network maps independent variables, such as spatial and
time coordinates, to a surrogate solution of a PDE. By applying
the chain rule, e.g., through auto-differentiation integrated in
2 Methodology
many deep-learning packages, one can easily extract the deriv- The main idea of our work, shown in Fig. 1, consists of two
atives of the output fields with respect to the input coordinates blocks. The first, MaxwellNet, is a neural network that takes
and consequently construct a physics-based loss.20 The correct as an input the refractive index distribution nðrÞ and predicts
prediction can be therefore retrieved by minimizing the loss with the scattered field Us . Its structure is based on the U-Net archi-
respect to the network weights. This approach has been used to tecture,28 and the training is performed on a large dataset of dig-
solve nonlinear differential equations,21–24 to realize the forward ital phantoms using a physics-defined loss function. Then, this
model in the inverse design of optical components,25 and to network is used as a forward model in a second optimization
extract material parameters in near-field microscopy.26 task that compares the fields predicted by MaxwellNet for a can-
Having the independent variables of PDE as the input of the didate refractive index (RI) distribution with the ground truth
neural network limits the use of PINNs when fast inference is projections, e.g., computed numerically or evaluated experi-
required. For the example of optical scattering, the neural net- mentally, and updates nðrÞ up to convergence.
work should be trained for each refractive index distribution
separately. A different idea was proposed recently in Ref. 27 2.1 Forward Model: MaxwellNet
to solve Maxwell’s equations for micro lenses with different
permittivity distributions. The calculation of physical loss, in In this section, we describe the implementation of a PINN that
this case, is based on the finite difference scheme, and in con- predicts the scattered field for a known input RI distribution. For
trast to the previous approach that is trained for a single exam- the sake of simplicity, we first describe the method for the 2D
ple, this model proved to be well-suited for cases in which fast case, but we will show the extension to 3D in the following. In
inference is required. However, such a PINN was only demon- this case, MaxwellNet takes as an input the RI distribution as a
strated to work for homogeneous 2D samples. discrete array of shape N x × N z × 1 and we do expect an output

Fig. 1 Schematic description of MaxwellNet, with U-Net architecture, and its application for tomo-
graphic reconstruction. The input is a refractive index distribution and the output is the envelope
of the scattered field. The output is modulated by the fast-oscillating term e jk 0 n 0 z to compute
the physics-informed loss for tuning the weights. For tomographic reconstruction, we minimize
a data-driven loss based on the difference between measured and calculated projections using
MaxwellNet. A regularization term can be added to improve the reconstruction.

Advanced Photonics 066001-2 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

with size N x × N z × 2, where the two channels correspond to scratch. This process, referred to in the following as fine-tuning,
the real and imaginary parts of the complex field. Among all is much faster than the original training of MaxwellNet. We will
the available architectures, the choice of U-Net appears favor- elaborate and discuss this interesting feature in Sec. 3.
able as we do expect to embed the latent features of the RI It should be mentioned that we train MaxwellNet based on
distributions in a lower dimensional space through consecutive the Helmholtz equation with a scalar field approximation, as
2D convolutions and then retrieve the complex electromagnetic described in Eq. (1). The scalar approximation allows us to have
field in the same spatial points through the decoding step. A a network with 2-channel output, representing the real and
similar architecture was also proven to provide good accuracy imaginary parts of the scalar field. We can also consider the
for the computation of the scattered field from micro lenses.27 full-vectorial Helmholtz equation where we need a larger net-
We implement the present network in TensorFlow 2.6.0. For work with 6-channel output to represent the real and imaginary
each step in the encoder, we use two Conv2D layers, each fol- parts of the three components of the field vector. However, the
lowed by batch-normalization and the elu activation function. A depolarization term can be neglected for samples with low re-
total number of five layers are adopted to encode the informa- fractive index gradients,31,32 allowing us to have a MaxwellNet
tion and each one is terminated with average pooling to reduce with fewer parameters and the scalar Helmholtz equation as the
the dimension. The maximum number of channels that we get in loss function.
the latent space is 512. On the decoder side, we used transposed
convolutional layers to the output with the size N x × N z × 2 (or
N x × N y × N z × 2 in the 3D case). It should be noted that we 2.2 Optical Diffraction Tomography Using MaxwellNet
also use residual skip connections from the encoder branch. In
common data-driven training, we would tune the weights of this Once MaxwellNet has been trained on a class of RI distribu-
network by minimizing the difference between predictions and tions, it can be used to rapidly backpropagate reconstruction er-
ground truth data computed with numerical solvers, in turn re- rors with an approach similar to learning tomography.6 Let us
quiring a large database of simulations and consequently assume that we measure L projections U m i , with i ¼ 0; …; L,
a massive computational cost. Here, we do not provide input- from an unknown RI distribution nðx; zÞ for different rotational
output pairs; instead we train the network by requiring that the angles. From these data, we can reconstruct a first inaccurate
Helmholtz equation is satisfied on the predicted field. To speed candidate nðx; zÞ through the Wolf’s transform using the
up the training and improve performances, we require the net- Rytov approximation. Then, we need to calculate the projec-
work to predict the slowly varying envelope of the scattered tions by MaxwellNet for different illumination angles. To imple-
field U senv being the scattered field obtained after demodulated ment illumination angle rotation, we can geometrically rotate
by the fast-oscillating component along the propagation direc- nðx; zÞ based on the corresponding illumination angle and cal-
tion Us ¼ U senv ejk0 n0 z . We define a physics-informed loss func- culate the scattered field for the rotated refractive index. By
tion to be minimized by updating the weights of the network: feeding MaxwellNet with ni ðx; zÞ ¼ Ri fnðx; zÞg, where Ri is
the image rotation operator that corresponds to the i-th projec-
X1 tion, we predict the complex scattered fields U si for the same
LPh ¼ k½∇2 þ k20 n2 ðrÞUs þ k20 ½n2 ðrÞ − n20 U i k2 ; (1) L angles. Consequently, we can construct a data-driven loss
r
N
function LD given by the difference kUsi − Um i k plus any addi-
2

where k0 is the wave number which is k0 ¼ 2π∕λ and tional regularizer, compute its gradient through auto-differentia-
λ ¼ 1.030 μm is the wavelength. nðrÞ is the RI distribution tion, update nðx; zÞ, and iterate up to convergence:
and n0 is the RI of the background medium. The summation
in Eq. (1) is done over the pixels of the computational domain X
L
1
and N is the number of pixels. To implement the Laplacian in LD ¼ kU s ðRi fngÞ − U m s
i k þ Regfn; U ðnÞg;
2
(3)
i¼1
L
Eq. (1), we follow the Yee grid finite difference scheme, com-
puting the derivative of the variables by 2D convolutions with a
5 × 5 kernel.29 Additionally, light scattering is by defining an ∂LD
open boundary problem. To satisfy the Sommerfeld radiation n → n − γD : (4)
condition and confine the problem in a finite space we use ∂n
a stretched-coordinate perfectly matched layer (PML)30 at the Also, in this case, we use an Adam optimizer for updating
edges of the simulation domain by introducing a complex the RI values. The regularizer in Eq. (3) consists of three
coordinate transformation [x → x þ ifðxÞ] when calculating parts: a total-variation (TV), a non-negativity, and physics-
the derivatives inside the PML region. We use the gradient of informed terms, Regfn;U sl ðnÞg ¼ λTV RTV ðnÞ þ λNN RNN ðnÞþ
the so-computed physical loss function to update the weights of λPh LPh ðn; Us Þ. The TV regularizer helps smooth the RI
the neural network, w through an Adam optimizer: reconstruction and the non-negativity regularizer adds the prior
information that nðx; zÞ should be larger than the background
∂LPh refractive index:
w → w − γ Ph : (2)
∂w
X qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
When we train MaxwellNet for a class of samples, it can ac- RTV ðnÞ ¼ j∇x nðrÞj2 þ j∇y nðrÞj2 þ j∇z nðrÞj2 ; (5a)
curately calculate the field for unseen samples from the same r
class. However, the key point to mention is that if we want
to use MaxwellNet for a different set of RI distributions, we X
can fix some of the weights and adjust only a part of the network RNN ðnÞ ¼ min ½nðrÞ − n0 ; 02 : (5b)
for the new dataset instead of re-starting the training from r

Advanced Photonics 066001-3 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

Importantly, we have to remark that MaxwellNet is trained It should be noted that once MaxwellNet is trained, the scat-
for a specific dataset and accurately predicts the scattered field tered field calculation is much faster than numerical techniques
for RI distributions that are not too far from this set. To take into such as FEM. We present a time comparison in Table 1. For
account this effect, we add the physics-informed loss to the reg- the test phantoms in Fig. 2, it took 17 ms for MaxwellNet in
ularizer. This further correction term helps to find RI values in comparison with 13 s for COMSOL, meaning three orders of
a way that MaxwellNet can predict the scattered field for them magnitude acceleration.
correctly. In contrast to TV and non-negativity constraints that Furthermore, performing a physics-based instead of direct
are used due to the ill-posedness of the ODT problem, the data-driven training holds promises for exploiting the advan-
physics-informed regularizer is necessary in our methodology tages of transfer learning.33 Maxwell equations are general
to ensure that the index distributions remain within the domain but having a neural network that predicts the scattered field
in which MaxwellNet has been trained. for any class of RI distribution in milliseconds with a negligible
The key advantages of using MaxwellNet with respect to physical loss is usually unfeasible. Most of the previous PINN
other forward models are: differently from BPM, it can accu- studies for solving partial differential equations are trained for
rately calculate field scattering, considering reflection, multiple- one example, and they will work for that specific example. In
scattering, or any other electromagnetic effects;5–8 once trained, our case, the U-Net architecture proved to be expressive enough
the field computations are performed in milliseconds, much to predict the field for a class of samples. However, if we use
faster than the Lipmann–Schwinger model; and finally, the MaxwellNet for inference on an RI distribution completely un-
data-driven error in Eq. (3) can be easily backpropagated differ- correlated with the training set, the accuracy drops. To evaluate
ently from commercially available full-vectorial solvers. We the MaxwellNet extrapolation capability, we considered the
discuss the reconstruction results and compare them with other model trained on phantoms samples in Fig. 2 and use it for in-
methods in Sec. 3.2. ference on HCT-116 cancer cells. The comparison between
MaxwellNet and COMSOL is shown in Fig. 2(c). The input
3 Results and Discussion of the network is a 2D slice of the experimentally measured
HCT-116 cell in the plane of the best focus. The discrepancy
3.1 MaxwellNet Results between MaxwellNet and COMSOL is due to the fact that
the former does not see examples of such RI distributions during
In this section, we evaluate the performance of MaxwellNet for the training. As a result, if we require accurate results for a new
the prediction of the scattered field from RI structures such as set of samples with different features, we have to re-train
biological cells. First, we check the performance on a 2D sample MaxwellNet for the new dataset, which would take a long time
assuming that the system is invariant along the y axis. The num- as shown in Table 1. However, it turns out that learning a physi-
ber of pixels for our model is N x ¼ N z ¼ 256 for both the cal law, as Maxwell equations, even though on a finite dataset, is
x and z directions, and their size is dx ¼ 100 nm. We also use better suited than data-driven training for transfer learning on
PML with a thickness of 1.6 μm at the edges of our computa- new batches. Indeed, we can use the pretrained MaxwellNet
tional domain. We create a dataset of digital cell phantoms and on digital phantoms and fine-tune some parts of the network
divide it into the training and testing sets. MaxwellNet has for HCT cells achieving good convergence in a few epochs.
∼5.9 × 106 parameters to train and we use the Adam optimizer In this example, we create a dataset of 136 RI distributions
with a learning rate of 1 × 10−4 and batch training. The details of HCT-116 cancer cells and divide them into the training
about the dataset and training and validation losses are discussed and validation sets. Some examples of the HCT-116 refractive
in Appendix B. We train and test MaxwellNet in TensorFlow index dataset are shown in Appendix B. A wide range of cells
2.6.0 on a desktop computer (Intel Core i7-9700K CPU, with different shapes are included in the dataset. We have single
3.6 GHz, 64 GB RAM, GPU GeForce RTX 2080Ti). cancer cells, such as shown in Fig. 2(c), examples of cells in the
In Figs. 2(a) and 2(b), we choose two random examples of mitosis process, or examples with multiple cells. In this case, we
the digital phantoms in the test set (which is not seen by the freeze the weights of the encoder part and fine-tune the decoder
network during the training). For each test case, in the second with the new dataset. We can see in Fig. 2(d) that after this cor-
and third rows, we present the prediction of the envelope of the rection step, the calculated field is much more accurate. As can
scattered field by the network, and we compare it with the result be seen in Table 1, the fine-tuning process is two orders of mag-
achieved by the FEM using COMSOL Multiphysics 5.4. We can nitude faster than a complete training from scratch.
see a very small difference between the results of MaxwellNet The 2D case is helpful for demonstrating the method and rap-
and COMSOL, which we attribute to discretization errors. idly evaluating the performances. Nevertheless, full 3D fields
There are different schemes of discretization in two methods are required for many practical applications. We can straightfor-
that can cause such differences. To quantitatively evaluate the wardly recast MaxwellNet in 3D using arrays of size N x × N y ×
performance of MaxwellNet, we define the relative error of N z × 1 as inputs of the network and requiring an N x × N y × N z ×
MaxwellNet with respect to COMSOL as 2 output, with the two channels corresponding to the real and
R imaginary parts of the envelope of the scattered field. In this
kU MaxwellNet ðrÞ − UCOMSOL ðrÞk2 dr case, the network consists of Conv3D, AveragePool3D, and
Error ¼ R ; (6)
kUCOMSOL ðrÞk2 dr Conv3DTranspose layers instead of 2D counterparts. As a
benchmark test, we created a dataset of 3D phantoms with
where UMaxwellNet and U COMSOL are the total fields calculated 200 examples (180 for training and 20 for testing). The compu-
with MaxwellNet and COMSOL. The integration is done ex- tational domain is defined with N x ¼ N y ¼ N z ¼ 64, dx ¼
cluding the PML regions. The calculated relative errors for test 100 nm, and PML thickness of 0.8 μm. To show the proof
case 1 and test case 2 in Fig. 2 are 4.1 × 10−2 and 4.6 × 10−2 , of concept of 3D MaxwellNet with limited computational re-
respectively. sources, we used a lower number of pixels per dimension with

Advanced Photonics 066001-4 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

(a) (b)

(d)

Fig. 2 Results of MaxwellNet and its comparison with COMSOL. (a) and (b) Two test cases from
the digital phantom dataset and the prediction of the real and imaginary parts of the envelope of the
scattered fields using MaxwellNet, COMSOL, and their difference. (c) Scattered field predictions
from the network trained in (a) and (b) for the case of an experimentally measured RI of HCT-116
cancer cells and comparison with COMSOL. The difference between the two is no longer
negligible. (d) Comparison between MaxwellNet and COMSOL after fine-tuning the former for
a set of HCT-116 cells. MaxwellNet predictions reproduces much more accurate results after
fine-tuning.

Table 1 Computation time comparison.

Dataset 2D phantoms (training) 2D HCT-116 (fine-tuning) 3D phantoms (training)

Training details 2700 samples 122 samples 180 samples
5000 epochs 600 epochs 5000 epochs
MaxwellNet training/fine-tuning 30.5 h 0.18 h 15.5 h
MaxwellNet inference 17 ms 17 ms 44.9 ms
COMSOL 13 s 13 s 2472 s

Advanced Photonics 066001-5 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

Fig. 3 Results of 3D MaxwellNet and its comparison with COMSOL. The RI distribution is shown
in (a). The real part of the envelope of the scattered field calculated by 3D MaxwellNet is shown in
(b), calculated by COMSOL in (c), and their difference in (d). The imaginary part of the envelope of
the scattered field calculated by 3D MaxwellNet, COMSOL, and their difference are presented in
(e)–(g), respectively.

respect to the 2D case, keeping the pixel size, dx, the same to subsection, we demonstrate how this method can be applied for
have an accurate finite difference calculation. As a result, we improving ODT reconstruction fidelity.
have a limited computational domain size, which can be im-
proved using more powerful resources.
3.2 Tomographic Reconstruction Results
The 3D version of MaxwellNet has ∼1.72 × 107 parameters.
We use the Adam optimizer (learning rate ¼ 1 × 10−4 ) and a To show the ability of MaxwellNet to be used for different
batch size of 10. The results of the predicted field for an unseen imaging applications, we implement an optimization task with
example and its comparison with COMSOL are shown in Fig. 3. MaxwellNet as the forward model for ODT, as explained in
We can see that MaxwellNet performs as well as COMSOL Sec. 2.2. In this example, we consider one of the digital phan-
in field calculation. The quantitative error described in Eq. (6) toms in the test set of Fig. 2, and we use 2D MaxwellNet as
is 3.4 × 10−3 for the 3D example of Fig. 3. There are some the forward model to compute the 1D scattered field along the
marginal differences due to different discretization schemes. transverse direction x for N ¼ 81 different rotation angles θ.
However, we can see in Table 1 that MaxwellNet is about We restrict ourselves to the range θ ∈ ½−40 deg; 40 deg as
50,000 times faster than COMSOL in predicting fields (44.9 ms this is consistent with the typical conditions in common
versus 41.2 min). This result and the significant efficiency in tomographic setups. As is shown in Fig. 4(a), the Rytov
the computation time highlight the MaxwellNet potential for reconstruction obtained from these field projections is elongated
the calculation of the field in different applications. In the next along the z axis and underestimated due to missing frequencies.

Advanced Photonics 066001-6 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

(a) Rytov Reconstruction Ground-truth (b) RI profile

x
z

5µm

Fig. 4 Tomographic reconstruction of RI using MaxwellNet. (a) The RI reconstruction was

achieved by Rytov, MaxwellNet, and the ground truth. (b) 1D RI profile at z ¼ 0 (plane of best
focus), for Rytov (green), MaxwellNet (blue), and the ground truth (orange).

We then minimize the loss function [Eq. (3)] to improve the MaxwellNet, which is slightly pixelated. This issue happens
RI reconstruction choosing as regularizer parameters λTV ¼ because the beam propagation method, as the forward model
3.1 × 10−7 , λNN ¼ 1 × 10−1 , λPh ¼ 5 × 10−2 , and the Adam op- of learning tomography, is a smooth forward model with respect
timizer with an initial learning rate of 3 × 10−4 . We also sched- to the voxels of the refractive index distribution, which is not
uled the learning rate, halving it every 1000 epochs to speed up the case for a deep neural network such as MaxwellNet.
convergence. The resulting RI distribution after 3000 epochs is However, the reconstructions are quantitatively comparable.
shown in Fig. 4. It can be seen that the reconstructed RI is no If we assume the reconstruction error of εðnrecon ; ntruth Þ ¼
longer underestimated nor elongated along the z axis. This is a knrecon − ntruth k22 ∕kntruth − n0 k22 , we get an error of 0.613 for
significant improvement in comparison with the Rytov predic- Rytov, an error of 0.146 for learning tomography, and an error
tion. The missing details in the reconstructed RI, which can be of 0.116 for MaxwellNet reconstructions, as shown in Fig. 5. In
better visible in the 1D cutline in Fig. 4(b), can be due to the terms of computation time with the desktop specifications we
missing information in the 1D fields that the optimization of RI mentioned earlier, we used 3000 epochs for iterative optimiza-
could not retrieve. tion with MaxwellNet, each epoch taking 570 ms and 600
Next, we try a 3D digital phantom from the test set and we epochs for learning tomography, each epoch taking 710 ms,
use the 3D MaxwellNet as the forward model in our tomo- which means a four-fold factor for MaxwellNet in the compu-
graphic reconstruction method. Since generating synthetic data tation time.
with COMSOL is time-consuming for multiple angles, we We also evaluated our methodology experimentally. We
create synthetic scattered fields from the phantom with the mentioned earlier that MaxwellNet takes care of reflection as
Lippmann–Schwinger equation. 9 Later, we will show an exper- a forward model, and therefore, our reconstruction technique
imental example where we illuminate the sample with a circular can be used for samples with high contrast. In our experimental
illumination pattern with an angle ≈10 deg. As a result, in this analysis, we try a polystyrene microsphere immersed in water,
numerical example, we rotate the sample for 181 angles (includ- where we expect to have a ∼0.25 refractive index contrast.
ing 1 normal incidence), equivalently to an illumination rotation Polystyrene microspheres (Polybead Polystyrene 2.0 Micron)
with a fixed illumination angle of 10 deg. We keep the exper- are immersed in water and placed between two #1 glass
imental conditions, λ ¼ 1.030 μm and n0 ¼ 1.33. Then, we use coverslips. We have an off-axis holographic setup where we
these synthetic measurements for our optimization task along use a yttrium-doped fiber laser (Amplitude Laser Satsuma)
with TV, non-negativity, and physics-informed regularization. with λ ¼ 1.030 μm and we change the illumination angle with
The reconstruction is achieved after 6000 epochs with λTV ¼ two Galvo mirrors. Using a delay path, the optical lengths of
1.2 × 10−8 , λNN ¼ 2 × 101 , and λPh started with 5 × 10−1 and the reference and signal arms are matched. We measure holo-
divided by two in every 500 epochs. The reconstructions are grams for 181 illumination angles and extract the phase and
shown in Fig. 5 in YX, YZ, and XZ planes. The first row shows amplitude of the complex scattered fields using Fourier holog-
the Rytov reconstruction where we can see a significant raphy. More details about the experimental setup are discussed
underestimation and elongation along the z axis that is due in Appendix C. Then, we use the extracted scattered fields for
to the small illumination angle (10 deg). The details in the different projections for our optimization task to reconstruct the
reconstruction achieved using MaxwellNet are slightly blurred 3D RI distribution of the sample. The experimental projections
in comparison with the ground truth as a result of low resolution are 2D complex fields that are imaged in the center of the sample
with λ ¼ 1.030 μm. using a microscope objective lens, and we can propagate them
Additionally, we performed learning tomography6 for the in the background medium to calculate the scattered field in
synthetic measurements using 181 projections. The 3D tomo- any other plane, perpendicular to z axis, after the sample. This
graphic reconstruction using learning tomography is shown 2D field can be compared with the output of MaxwellNet
in the third row of Fig. 5. In comparison with MaxwellNet, in that plane, as described in Eq. (3). Additionally, the exper-
learning tomography has some elongated artifacts, which can imental projections are based on illumination rotation and we
be due to the fact that reflection is neglected in its forward interpolate them to achieve the equivalent sample rotation pro-
model. However, the reconstruction with learning tomogra- jections. We iteratively optimize the loss function in Eq. (3) for
phy is smoother in comparison with the reconstruction of 2000 epochs where we use the regularization parameters of

Advanced Photonics 066001-7 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

Fig. 5 Tomographic RI reconstruction of 3D sample using MaxwellNet. The RI reconstruction

is achieved by Rytov, MaxwellNet, learning tomography, and the ground truth in different rows
at the Y X , Y Z , and X Z planes in the center of the sample. A z-stack demonstration of the
reconstruction is shown in Video 1 (Video 1, MP4, 1.4 MB [URL: https://doi.org/10.1117/1.AP
.4.6.066001.s1]).

λTV ¼ 3.8 × 10−9 , λNN ¼ 5 × 101 , and λPh started with 1.5 × 10−1 4 Conclusion
and divided by two in every 500 epochs. The reconstruction is
shown in Fig. 6 using Rytov, MaxwellNet, and learning tomog- In summary, we proposed a PINN that rapidly calculates the
raphy. It can be seen that the underestimation and z axis elon- scattered field from inhomogeneous RI distributions such as
gation in the Rytov reconstruction are remarkably improved. biological cells. Our network is trained by minimizing a loss
The reconstruction using learning tomography in Fig. 6 has function based on Maxwell equations. We showed that the net-
artifacts due to the high refractive index contrast of the poly- work can be trained for a set of samples and could predict the
styrene bead and reflections that cannot be considered in the scattered field for unseen examples that are in the same class. As
beam propagation method. our PINN is not a data-driven neural network, it can be trained

Advanced Photonics 066001-8 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

Fig. 6 Tomographic RI reconstruction of a polystyrene microsphere immersed in water. The

projections are measured with off-axis holography for different angles. The RI reconstructions
achieved by Rytov, MaxwellNet, and learning tomography are presented at the Y X , Y Z , and
X Z planes in the center of the sample, respectively.

for different examples in different conditions. Even though the variables can be calculated using the chain rule. In this
network is not efficiently extrapolating to classes that are sta- implementation, the weights of the network can be trained to
tistically very different from the training dataset, we showed that minimize the loss function for a single refractive index distribu-
by freezing the encoder weights and fine-tuning the decoder tion, nðrÞ in Eq. (1). In our approach, the PINN gets the refrac-
branch, one can get a new predictive model in a few minutes. tive index, nðrÞ, on a uniform grid as the input and finds the field
We believe that this can be further used for changing the wave- on the same grid which minimizes the loss function for that re-
length, boundary condition, or other physical parameters. fractive index. The output of the network is the 3D array of the
We used our PINN as a forward model in an optimization scattered field envelope, and we use a finite difference scheme
loop to retrieve the RI distribution from the scattered fields to calculate the derivative of the field with respect to the co-
achieved by illuminating the sample from different directions, ordinates:
known as ODT. This example shows the ability of MaxwellNet
to be used as an accurate forward model in optimization loops for
∂Us Us ðði þ 1ÞΔx; jΔy; kΔzÞ − Us ðði − 1ÞΔx; jΔy; kΔzÞ
inverse design or inverse scattering problems. ¼ ;
∂x 2Δx
5 Appendix A: Calculation of Physics- (7)
Informed Loss in which ði; j; kÞ are the pixel indices, and Δx, Δy, and Δz are
During the training of MaxwellNet, we calculate at each epoch the pixel sizes along the x-, y-, and z-axes. This way, we can
s
the loss function in Eq. (1) for the network output. To evaluate calculate ∂U s
∂x by convolving U with a kernel of ½−1∕2, 0; 1∕2
the Helmholtz equation residual, we should numerically com- along the x axis. When computing electromagnetic fields, since
2 s 2 s 2 s
pute the term ∂∂xU2 þ ∂∂yU2 þ ∂∂zU2 . In the previous PINN papers the curl of the electric field gives the magnetic field and vice
for solving PDEs,19–21,23–26 the inputs of the network are the spa- versa, a smart technique to improve accuracy is to use two stag-
tial coordinates x, y, z, and the derivatives with respect to these gered grids for discretizing fields, commonly referred to as the

Advanced Photonics 066001-9 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

Yee scheme.29 In practice, this can be easily implemented size of 10 for 5000 epochs. This training took 30.5 h and after
through two shifted convolutional kernels for the two grids, 5000 epochs, no significant decrease in the validation loss could
½−1∕2,1∕2,0 and ½0; −1∕2,1∕2. be observed. The training and validation curves of the physical
To minimize the discretization error, one can use a smaller loss are shown in Fig. 7(a). This figure shows that MaxwellNet
pixel size, Δx, or higher-order approximations. Here, we use performs very well for out-of-sample cases.
the fourth-order finite difference scheme34 in which convolu- We discussed in Sec. 3 using MaxwellNet that was trained
tional kernels of ½0; þ1∕24; −9∕8; þ9∕8; −1∕24 and ½þ1∕24; for cell phantoms to predict the scattered field for real cells. A
−9∕8; þ9∕8; −1∕24, 0 are used for the calculation of the dataset of HCT-116 cancer cells is used for this purpose. The 3D
derivatives in Eq. (1). refractive index of these cells is reconstructed using the Rytov
approximation, with projections achieved with an experimental
6 Appendix B: Training and Fine-tuning of setup utilizing a spatial light modulator, as described in Ref. 8.
Then, a 2D slice of the refractive index is chosen in the plane of
MaxwellNet best focus. A total number of 8 cells are used, and we rotated
As mentioned in Sec. 3, we create a dataset of digital cell phan- and shifted these cells to create a dataset of 136 inhomogeneous
toms to train and validate MaxwellNet. The dataset for 2D cells whose refractive index range is (1.33, 1.41). We use 122 of
MaxwellNet includes 3000 phantoms with elliptical shapes these images for training and 14 for validation. Some examples
oriented in different directions. The size of these phantoms is of the HCT-116 refractive index dataset are shown in Fig. 7(c).
in the range from 5 to 10 μm, their refractive index varies in We freeze the encoder of MaxwellNet and fine-tune its decoder
the range of (1.38, 1.45), and the background refractive index for this new dataset. The training and validation losses are
is n0 ¼ 1.33. Two examples of these phantoms are shown in shown in Fig. 7(b).
Fig. 2. We divide this dataset into 2700 phantoms for training For 3D MaxwellNet, a dataset of 200 phantoms is created.
and 300 phantoms for testing. We use batch training with a batch These 3D phantoms have a spherical shape with some details

(a) Training (b) Fine-tuning (c)

Fig. 7 Training and fine-tuning of MaxwellNet. (a) Training (blue) and validation (orange) loss of
MaxwellNet for the digital cell phantoms dataset. (b) Fine-tuning of the pretrained MaxwellNet for a
dataset of HCT-116 cells for 1000 epochs. (c) Examples of the HCT-116 dataset.

Fig. 8 Experimental setup for multiple illumination angle off-axis holography. HW: half-wave plate;
P: polarizer; BS: beam splitter; L: lens; Obj: microscope objective; and M: mirror.

Advanced Photonics 066001-10 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

inside them and the range of their diameter is 1.8 to 2.4 μm. We 12. Y. Rivenson et al., “Deep learning microscopy,” Optica 4(11),
randomly choose 180 phantoms for training and 20 phantoms 1437–1443 (2017).
for testing. We train 3D MaxwellNet with the training dataset 13. N. Borhani et al., “Learning to see through multimode fibers,”
with a batch size of 10. The example of Figs. 3 and 5 is one of Optica 5(8), 960–966 (2018).
14. B. Rahmani et al., “Multimode optical fiber transmission with
the phantoms in the testing dataset. a deep learning network,” Light Sci. Appl. 7(1), 69 (2018).
15. Y. Rivenson et al., “Phase recovery and holographic image
7 Appendix C: Experimental Setup for ODT reconstruction using deep learning in neural networks,” Light
For ODT, we require complex scattered fields from multiple il- Sci. Appl. 7(2), 17141 (2018).
16. J. Lim, A. B. Ayoub, and D. Psaltis, “Three-dimensional tomog-
lumination angles. The off-axis holographic setup to accomplish raphy of red blood cells using deep learning,” Adv. Photonics 2(2),
that is shown in Fig. 8. It relies on an ytterbium-doped fiber laser 026001 (2020).
at λ ¼ 1.030 μm whose power is controlled with a half-wave 17. Z. Ren, Z. Xu, and E. Y. Lam, “End-to-end deep learning frame-
plate (HW) and a polarizing beam splitter. The optical beam is work for digital holographic reconstruction,” Adv. Photonics 1(1),
divided into the signal and reference arms using a beam splitter 016004 (2019).
(BS1). In the signal arm, we use two galvo mirrors, GM-V and 18. D. Pirone et al., “Speeding up reconstruction of 3D tomograms in
GM-H, to control the illumination angle in the vertical and holographic flow cytometry via deep learning,” Lab Chip 22(4),
horizontal directions. Using two 4F systems (L1–L4), we image 793–804 (2022).
these galvo mirrors on the sample plane, so the position of the 19. I. E. Lagaris, A. Likas, and D. I. Fotiadis, “Artificial neural net-
works for solving ordinary and partial differential equations,”
beam remains fixed while changing the illumination angle. This
IEEE Trans. Neural Networks 9(5), 987–1000 (1998).
way, we can illuminate the sample with a condensed plane wave. 20. L. Lu et al., “DeepXDE: a deep learning library for solving differ-
The sample is then imaged on the camera (Andor sCMOS Neo ential equations,” SIAM Rev. 63(1), 208–228 (2021).
5.5) using another 4F system consisting of a 60× water dipping 21. M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed
objective (Obj1) and a tube lens L5. The signal and reference neural networks: a deep learning framework for solving forward
arms are then combined with another beam splitter, BS2, to and inverse problems involving nonlinear partial differential equa-
create the off-axis hologram on the camera. A motorized delay tions,” J. Comput. Phys. 378, 686–707 (2019).
line controls the optical path of the reference arm to match the 22. S. M. H. Hashemi and D. Psaltis, “Deep-learning PDEs with
optical path of the signal arm. unlabeled data and hardwiring physics laws,” arXiv:1904.06578
(2019).
23. Z. Mao, A. D. Jagtap, and G. E. Karniadakis, “Physics-informed
Acknowledgments neural networks for high-speed flows,” Comput. Methods Appl.
This project was funded by the Swiss National Science Mech. Eng. 360, 112789 (2020).
Foundation (SNSF) under funding number 514481. Additionally, 24. X. Jin et al., “NSFnets (Navier-Stokes flow nets): physics-
informed neural networks for the incompressible Navier-Stokes
A. Saba would like to thank Dr. Babak Rahmani and Dr. Joowon equations,” J. Comput. Phys. 426, 109951 (2021).
Lim for their useful comments about the implementation of 25. Y. Chen et al., “Physics-informed neural networks for inverse
MaxwellNet and tomographic optimization in TensorFlow. The problems in nano-optics and metamaterials,” Opt. Express 28(8),
authors have no competing interests to declare. 11618–11633 (2020).
26. Y. Chen and L. Dal Negro, “Physics-informed neural networks for
References imaging and parameter retrieval of photonic nanostructures from
near-field data,” APL Photonics 7(1), 010802 (2022).
1. W. Choi et al., “Tomographic phase microscopy,” Nat. Methods 27. J. Lim and D. Psaltis, “Maxwellnet: physics-driven deep neural
4(9), 717–719 (2007). network training based on maxwell’s equations,” APL Photonics
2. Y. Sung et al., “Optical diffraction tomography for high resolution 7(1), 011301 (2022).
live cell imaging,” Opt. Express 17(1), 266–277 (2009). 28. O. Ronneberger, P. Fischer, and T. Brox, “U-Net: convolutional net-
3. D. Jin et al., “Tomographic phase microscopy: principles and works for biomedical image segmentation,” Lect. Notes Comput.
applications in bioimaging,” J. Opt. Soc. Am. B 34(5), B64–B77 Sci. 9351, 234–241 (2015).
(2017). 29. K. Yee, “Numerical solution of initial boundary value problems
4. E. Wolf, “Three-dimensional structure determination of semi- involving Maxwell’s equations in isotropic media,” IEEE Trans.
transparent objects from holographic data,” Opt. Commun. 1(4), Antennas Propag. 14(3), 302–307 (1966).
153–156 (1969). 30. W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium
5. J. Lim et al., “Comparative study of iterative reconstruction algo- from modified Maxwell’s equations with stretched coordinates,”
rithms for missing cone problems in optical diffraction tomogra- Microwave Opt. Technol. Lett. 7(13), 599–604 (1994).
phy,” Opt. Express 23(13), 16933–16948 (2015). 31. A. Ishimaru, Wave Propagation and Scattering in Random Media,
6. U. S. Kamilov et al., “Learning approach to optical tomography,” Vol. 2, Academic Press, New York (1978).
Optica 2(6), 517–522 (2015). 32. A. Saba et al., “Polarization-sensitive optical diffraction tomogra-
7. S. Chowdhury et al., “High-resolution 3D refractive index micros- phy,” Optica 8(3), 402–408 (2021).
copy of multiple-scattering samples from intensity images,” 33. C. Tan et al., “A survey on deep transfer learning,” Lect. Notes
Optica 6(9), 1211–1219 (2019). Comput. Sci. 11141, 270–279 (2018).
8. J. Lim et al., “High-fidelity optical diffraction tomography of 34. A. Fathy et al., “A fourth order difference scheme for the Maxwell
multiple scattering samples,” Light Sci. Appl. 8(1), 1 (2019). equations on Yee grid,” J. Hyperbol. Differ. Equ. 5(3), 613–642
9. T.-A. Pham et al., “Three-dimensional optical diffraction tomog- (2008).
raphy with Lippmann–Schwinger model,” IEEE Trans. Comput.
Imaging 6, 727–738 (2020).
10. G. E. Karniadakis et al., “Physics-informed machine learning,” Amirhossein Saba is a PhD student in photonics at École Polytechnique
Nat. Rev. Phys. 3(6), 422–440 (2021). Fédérale de Lausanne (EPFL), Lausanne, Switzerland. He received his
11. S. Cai et al., “Physics-informed neural networks (PINNs) for fluid BS degree in electrical engineering with a minor in physics from Sharif
mechanics: a review,” Acta Mech. Sin. 37, 1727–1738 (2021). University of Technology, Tehran, Iran, in 2018. His current research

Advanced Photonics 066001-11 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
 Saba et al.: Physics-informed neural networks for diffraction tomography

interests include computational imaging, polarization-sensitive, and non- de Lausanne, Lausanne, Switzerland. His current research interests in-
linear imaging techniques. clude optical imaging and three-dimensional refractive index reconstruc-
tions for biological samples.
Carlo Gigli received his MS degree in nanotechnologies for ICTs from
Politecnico di Torino and Université Paris Diderot in 2017, and his PhD Demetri Psaltis received his BSc, MSc, and PhD degrees from
in physics at Université de Paris in 2021. During this period, his research Carnegie-Mellon University, Pittsburgh, Pennsylvania. He is a professor
interests included the design, fabrication, and characterization of dielec- of optics and director of the Optics Laboratory at the École Polytechnique
tric resonators and metasurfaces for nonlinear optics. Since September Fédérale de Lausanne (EPFL), Lausanne, Switzerland. In 1980, he
2021, he has been working as post-doc in the Optics Laboratory at EPFL, joined the faculty at the California Institute of Technology, Pasadena,
where his main activities focus on optical computing, physics informed California. He moved to EPFL in 2006. His research interests include
neural network, and nonlinear tomography. imaging, holography, biophotonics, nonlinear optics, and optofluidics.
He has authored or coauthored over 400 publications in these areas.
Ahmed B. Ayoub received his BS degree in electrical engineering from He is a fellow of the Optical Society of America, the European Optical
Alexandria University in Egypt in 2013. He received his MS degree in Society, and SPIE. He was the recipient of the International Commission
physics from the American University, Cairo, Egypt, in 2017. He is pursu- of Optics Prize, the Humboldt Award, the Leith Medal, and the Gabor
ing his PhD in electrical engineering at the École Polytechnique Fédérale Prize.

Advanced Photonics 066001-12 Nov∕Dec 2022 • Vol. 4(6)

Downloaded From: https://www.spiedigitallibrary.org/journals/Advanced-Photonics on 03 Jul 2023
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use