Total Variation-Based Reconstruction and Phase Retrieval for Diffraction Tomography with an Arbitrarily Moving Object

We consider the imaging problem of the reconstruction of a three-dimensional object via optical diffraction tomography under the assumptions of the Born approximation. Our focus lies in the situation that a rigid object performs an irregular, time-dependent rotation under acoustical or optical forces. In this study, we compare reconstruction algorithms in case i) that two-dimensional images of the complex-valued wave are known, or ii) that only the intensity (absolute value) of these images can be measured, which is the case in many practical setups. The latter phase-retrieval problem can be solved by an all-at-once approach based utilizing a hybrid input-output scheme with TV regularization.


Introduction
We are interested in the tomographic reconstruction of a three-dimensional (3D) rigid object, which is illuminated from various directions.In optical diffraction tomography, the wavelength of the imaging wave, visible light with wavelength hundreds of nanometers, is in a similar order of magnitude as features of the µm-sized object we want to reconstruct.This setup substantially differs from the X-ray computerized tomography, where the wavelength is much shorter and therefore the light can be assumed to propagate along straight lines.x 3 " r M Figure 1: Experimental setup with measurement plane located at x 3 " r M , see [3].
Formally, in optical diffraction tomography we want to recover an object or, more specifically, its scattering potential f : R 3 Ñ R, which has compact support.At each time step t P r0, T s, the object undergoes a rotation described by a rotation matrix R t P SOp3q, such that its scattering potential becomes f pR t ¨q.Here, the rotation is not restricted to some direction, we only require that R t depends continuously differentiably on t, which is the case for objects in acoustical traps, cf.[1].We further assume that the rotation R t is known beforehand; for the detection of motion from the tomographic data we refer to [2].The object is illuminated by an incident, plane wave u inc pxq :" e ik 0 x 3 , x " px 1 , x 2 , x 3 q P R 3 , of wave number k 0 ą 0, which induces a scattered wave u t : R 3 Ñ C. We measure the resulting total field u tot t pxq " u t pxq `uinc pxq, x P R 3 , at a fixed measurement plane x 3 " r M , see Fig. 1.
Under the Born approximation, the scattered wave u t is a solution of the partial differential equation ´p∆ `k2 0 qu t " pu t `uinc q f pR t ¨q, see [4].We assume the Born approximation to be valid, which is the case for small, mildly scattering objects, see [5, § 3.3].Then the relationship between the measured field u t and the desired function f can be expressed via the Fourier diffraction theorem [6,7] by where Here the 3D Fourier transform F and the partial Fourier transform F 1,2 in the first two coordinates are given by f px 1 , x 3 q e ´i x 1 ¨k1 dx 1 for k P R 3 and k 1 P R 2 .The left-hand side of ( 1) is fully determined by the measurements u tot t p¨, ¨, r M q, and the right-hand side provides non-uniform samples of the Fourier transform Frf s, evaluated on the union of semispheres which contain the origin and are rotated by R t .Therefore, the reconstruction problem in diffraction tomography can be seen as a problem of inverting the 3D Fourier transform F given non-uniformly sampled data in the so-called k-space.If these samples have a positive Lebesgue measure, the scattering potential is uniquely defined by the measurements.
Theorem 1.1 (Unique inversion, [3,Thm. 3.2]) Let f P L 1 pR 3 q have a compact support in tx P R 3 : }x} 2 ď r M u, and let the Lebesgue measure of be positive.Then f is uniquely determined by u tot t p¨, ¨, r M q, t P r0, T s.
In § 2 of this paper, we consider algorithms for the tomographic reconstruction of f from u tot t p¨, ¨, r M q, t P r0, T s.Since in practice one can often measure only intensities |u tot t |, cf.[8], in § 3 we will incorporate phase retrieval methods to retrieve f .In contrast to [3], we focus here on the numerical studies for the case of an arbitrarily moving rotation axis.

Reconstruction methods
We consider three different reconstruction methods for the situation that we know the complex-valued field u t p¨, ¨, r M q.First, the filtered backpropagation is based on a truncation and discretization of the inverse 3D Fourier transform, where Y is the set where the data Ff is available in k-space.We state the following filtered backpropagation formula for a moving rotation axis, the special case for constant axis is well established [9].Theorem 2.1 (Filtered Backpropagation, [6,Thm. 4.1]) Let the assumptions of the Fourier diffraction theorem be satisfied.Let α P C 1 pr0, T s, Rq and n P C 1 pr0, T s, S 2 q.For t P r0, T s, we denote the rotation around the axis nptq with the angle αptq by R J nptq,αptq :" where nptq " pn 1 , n 2 , n 3 q J , c :" cos αptq and s :" sin αptq.Denote by u t the scattered wave according to the rotated object f ˝Rnptq,αptq .Next, let be the map that covers the accessible domain Y " ΦpUq in k-space.Then, for all x P R 3 , where Card denotes the cardinality of a set and |∇Φ| is the absolute value of the Jacobian determinant of Φ.
If nptq and αptq are known analytically in (3), then one can compute |∇Φ|, but the Banach indicatrix or Crofton symbol CardpΦ ´1pΦpk 1 , k 2 , tqqq is hard to determine analytically.However, since the backpropagation tends to produce artifacts due to the limited coverage Y in k-space, we consider alternative approaches.The second approach consists in finding arg min where F denotes a discretization of F, a nonuniform discrete Fourier transform (NDFT) [10], f are uniform samples of f and g is the uniformly sampled data on the left-hand side of (1), which is determined by the measurements.The resulting normal equation can be solved with a conjugate gradient (CG) method, making use of fast implementations [11] of the NDFT and its adjoint.Such methods have also been applied in X-ray imaging [12], magnetic resonance imaging [13], and spherical tomography [14,15].
Since the CG reconstruction is sensitive regarding to larger measurement errors in order to be applied in phase retrieval, we further add a total variation (TV) regularization term and incorporate the non-negativeness of the object, see [3, § 5.3].For the discretized object pf q :" f px q, x :" 1 L P R 3 , P I, on a uniform Cartesian grid, where I Ă Z 3 denotes the employed index set, the total variation is given by TVpf q :" see e.g.[16], where the discrete gradient grad may be computed using finite forward differences.To incorporate the non-negativity constraint, we exploit the characteristic function χ ě0 pf q :" # 0, if f ě 0 for all P I, `8, otherwise.
Instead of minimizing (4), we propose to invert the NDFT by computing arg min f χ ě0 pf q `1 2 }F pf q ´g} 2 2 `λ TVpf q. (5) The minimization problem (5) can be solved using a primal-dual method [17] with backtracking [18].We compare the different reconstruction methods numerically for the rotation axis nptq " `sin `1 2 π sin t ˘, cos `1 2 π sin t ˘, 0 ˘J and angle αptq " t, t P r0, 2πs, for which the backpropagation formula is given in [6,Ex. 4.8].Our forward data u tot t is created based on the Born approximation.In Fig. 2, we see that the filtered backpropagation falls behind the other methods.Adding 5 % Gaussian noise to the data u t `uinc , we see in Fig. 3 that the TV-regularized inverse is superior to the others.Here, we perform 20 iterations for the CG and 200 iterations for the primal-dual TV-regularization.The forward model to generate the data also uses Born approximation, but a different algorithm based on convolution.

Phase retrieval
In many practical applications of diffraction tomography, only the intensities ˇˇu tot t px 1 , x 2 , r M q ˇˇ, px 1 , x 2 q P R 2 , t P r0, T s, are available since the measurement of the phase is technically very challenging and expensive.The missing phase information of the complex-valued function u tot t thus has to be recovered in advance.Especially because the inversion of the given magnitudes equipped with an arbitrary phase would yield a solution of f , phase retrieval in diffraction tomography is highly ill posed.Such behaviour is also well known for classical phase retrieval problems, where the ambiguities can be characterized analytically [19,20,21].Therefore, we require a priori knowledge to obtain a meaningful reconstruction.For our numerical experiments, we exploit that f is nonnegative with bounded support, and that f often has a cartoon-like appearance, e.g., for images of biological cells.Unfortunately, our phase retrieval problem in diffraction tomography fits neither into the classical setting for instance studied by Gerchberg, Saxton, and Fienup [22,23] nor into the modern approaches of PhaseLift [24] and its tensor-free variant [25], since the connection between u tot t and f is not linear but affine.For phase retrieval in diffraction tomography, we propose an all-at-one approach utilizing the hybrid input-output (HIO) scheme, see [3, § 6], which is based on Fienup's HIO method [23].In a nutshell, the main idea behind the iterative HIO approach consists of computing the forward model, projecting to the known magnitudes in the measurement domain, solving the known-phase reconstruction problem with one of the above-mentioned methods and incorporating the additional a priori information on f .Instead of strictly enforcing the a priori constraints, the HIO algorithm tries to correct the last input to push the next output in direction of a feasible object.In more details, we apply the HIO method in Algorithm 1, where the operator D models the measurement process that maps the discretizations of f to u tot t p¨, ¨, r M q, t P r0, T s, the sign function sgnpzq :" is applied element by element, and d denotes the element-wise multiplication.
The crucial point of the proposed method is that when computing D ´1 we have to solve the known-phase reconstruction in each step, which is performed via an iterative algorithm by itself.Since both iterative methods show a slow convergence, we made several numerical improvements to speed up computation.For instance, we restart the inner loop of the primal-dual algorithm with the parameters and dual variables obtained in the previous outer step and execute only a small number of inner iterations.In this manner, both iterative algorithms may be incorporated in each other.Furthermore, an initial reconstruction using CG as inner algorithm serves as starting point for the more accurate but slower HIO with primal-dual iteration.It is well known that the performance of HIO may be improved by utilizing a total variation denoising between the iterates [26,27].If the total variation regularization is used during the backward transform, this TV denoising is automatically incorporated into the HIO method.
Numerical simulations indicate that phase retrieval in our diffraction tomography setup is indeed possible utilizing the proposed HIO algorithm with the primal-dual TV regularization, see Fig. 4 and 5.Here we use the same rotation axes (6) as before.For the HIO with CG, we performed J HIO " 10 outer and 5 inner (CG) iterations, where more iterations yield a noisier result.For the HIO with primal-dual, we performed J HIO " 200 outer and 5 inner iterations.We note that, in contrast to the reconstructions of § 2, we incorporate additional a priori information via the support constraint r s " 40.A Matlab implementation of the algorithms is available at https: //github.com/michaelquellmalz/FourierODT.

Conclusion
We have compared reconstruction algorithms in diffraction tomography with a rigid object that undergoes a time-dependent rotation according to a moving rotation axis.All three algorithms work well for exact data with an advantage for the iterative methods (CG and primal-dual).In the presence of noise in the measured data, the TV-regularized primal dual shows superior results.If only the intensity data |u tot | is known, the allat-once HIO algorithm combined with the TV regularization is still able to recover the object accurately.