## Abstract

Due to incompleteness of input data inherent to Limited Angle Tomography (LAT), specific additional constraints are usually employed to suppress image artifacts. In this work we demonstrate a new two-stage regularization strategy, named Generalized Total Variation Iterative Constraint (GTVIC), dedicated to semi-piecewise-constant objects. It has been successfully applied as a supplementary module for two different reconstruction algorithms: an X-ray type solver and a diffraction-wise solver. Numerical tests performed on a detailed phantom of a biological cell under conical illumination pattern show significant reduction of axial blurring in the reconstructed refractive index distribution after GTVIC is added. Analogous results were obtained with experimental data.

© 2016 Optical Society of America

## 1. Introduction

Quantitative Phase Imaging (QPI) techniques are the answer to a demand for nondestructive measurement methods designed to determine refractive index (RI) distribution of transparent and semi-transparent samples. Classical QPI methods for label-free RI investigation allow to obtain integrated data, where information about optical path difference is mixed with information about RI variance [1–3]. Recently these techniques have been enhanced to “2.5D” methods, where it is possible to retrieve object’s topography based on the integrated phase measurements [4,5]. However, these methods are applicable to objects that are either thin (with known object volume size) or have quasi-constant RI distribution (e.g. red blood cells). An important technique which meets the requirement of true 3D measurement is optical diffraction tomography (ODT) [6]. In typical full projection angle ODT a biological cell is inserted into a glass capillary filled with cell culture medium. The capillary is located in an immersion liquid, illuminated by a laser beam in an optical tomography setup, in which a series of optical projections is acquired for different angular positions of the capillary. Usually the angular range of rotation is 360° around the axis of symmetry of the capillary, which guarantees good quality reconstructions. Unfortunately, this approach requires relatively complicated and time-consuming process of placing cells into the capillary. Also, the process of acquiring projections is slow, which makes it impossible to investigate dynamic processes. What is more, the glass capillary introduces strong aberrations into the object beam wavefront [7–9].

Because of these problems, the limited angle tomography (LAT) is gaining much interest [10–14]. In LAT, biological microstructure is placed in a cell culture medium at a Petri dish or a glass plate which is then introduced into the optical tomography system built in a vertical configuration. Instead of rotating the object, in LAT the sample is illuminated from several directions, whereas the Petri dish is stationary. This approach has multiple advantages: (1) it minimizes vibrations in the system, (2) increases the measurement speed, since changing a position of laser illumination beam can be carried out much faster than object rotation, thus allowing analysis of biological dynamic processes, (3) it simplifies the preparation of a sample, since many biological cells are grown in Petri dishes. The biggest disadvantage is the fact, that projections from a limited angular range can be acquired only. Therefore, according to the Fourier Diffraction Theorem (FDT) [15], when object projections are recorded within a limited angular range, the Fourier spectrum contains areas, where no information about the measured object is present. When an inverse Fourier transform of such spectrum is calculated the results are corrupted due to two effects: the obtained reconstruction has anisotropic resolution (relatively high in *x-y* plane and much lower along *z* axis) and it is elongated in the *z* direction, where *z* is the system optical axis. In order to compensate these effects, special reconstruction algorithms are necessary.

Numerous iterative methods have been developed for limited angle tomography [16]. State-of-the-art methods include iterative constraint tomography techniques, where different types of constraints are applied in each iteration [6,16–19]. These approach give limited improvement of the reconstruction quality or are not suitable for optical tomography. Another approach utilizes Total Variation minimization (TVM), which iteratively minimizes the L1 norm of a gradient of a reconstruction. This technique is applied together with algebraic reconstruction methods, like SART (Simultaneous Algebraic Reconstruction Technique) or SIRT (Simultaneous Iterative Reconstructive Technique) [12,13,20,21]. Since TVM is a very strong regularizer, it can recover the original signal thus improving reconstruction quality to a great extent. It can be, however, used only with samples whose inner function can be approximated with a piecewise-constant function, which limits its applicability mainly to technological objects. Application of this method to investigation of biological microstructures results in corrupted reconstructions. Recently it was proved that utilization of *a priori* knowledge in the form of automatically calculated 3D mask with object’s external geometry to tomographic solver which assumes that light rays travel in straight lines through the object, enhances both convergence rate of the algorithm and quality of 3D reconstruction. The method is called Total Variation Iterative Constraint (TVIC) and can be applied to nonpiecewise-constant objects. It combines the advantages of TVM with the possibility of retrieving nonpiecewise-constant information [10,22]. In this paper we present a novel, comprehensive approach to reconstruction of 3D RI distribution of nonpiecewise-constant samples. The new method also uses 3D RI mask calculated automatically as an *a priori* information, but in contrast to the TVIC technique, it may work with an arbitrary tomographic solver and thus is applicable to very broad class of objects. It returns reconstructions with minimized deformation that otherwise would be present due to limited angular range of projections, and at the same time it allows for an arbitrary type of an inner function of a measured sample. Because of this versatility, it is called Generalized TVIC (GTVIC).

The paper is organized as follows. In Section 2 we describe the general strategy of GTVIC method. Next the methodology of numerical tests and their results are presented in Sections 3 and 4 respectively. The effectiveness of the GTVIC strategy when applied to LA ODT of biological cells is shown in Section 5, while the final conclusions are given in Section 6.

## 2. Generalized total variation iterative constraint strategy

Figure 1 presents the flow diagram of GTVIC strategy. The method is a modular approach and consists of two stages. The principle behind this method is to generate a constraint in the first module and apply it in the second module, which is an arbitrary iterative tomographic solver. The GTVIC constraint itself, adaptively optimized during the first stage, has a form of a 3D mask, defining the spatial limits of the object volume. The sample should thus be distinguishable from the surrounding medium in terms of investigated phenomena. GTVIC method requires a set of object projections acquired for different illumination directions and a list of projection vectors that specifies these directions, as an input. The output is a 3D distribution of investigated values of the sample with minimized artifacts that otherwise would be present due to limited angular range of object projections.

The GTVIC mask geometry is based on a preliminary TV-regularized reconstruction performed by the Chambolle-Pock algorithm [23]. In principle, this algorithm minimizes the sum of two functionals, as shown in Eq. (1):

where*A*– system matrix, which depends on the optical system geometry,

*x*– reconstruction,

*B*– input data (sinogram),

*λ*– weighting factor,

*||·||*– TV norm, which is equal to

_{TV}*||grad(·)||*. The second functional is TVM. It minimizes the first norm of the gradient of a reconstruction. In the case of biological microobjects whose inner function is smooth, we obtain a reconstruction with erroneously retrieved inner structures, but with good approximation of the real 3D external geometry – free from deformation which normally is present due to limited angular range of projection. This result is then binarized, with a threshold level equal to 0.7*

_{1}*h*, where

*h*is the threshold level calculated automatically with Otsu’s method [24]. Next, in order to achieve a proper 3D mask with object’s external geometry, the binarized mask must be apodized. This is carried out through morphological dilation of the binary mask. This operation guarantees that every fragment of object volume is inside the 3D mask. The mask is then passed to the second GTVIC module, where it is used as a powerful prior.

The second module of GTVIC method can be an arbitrary iterative solver for 3D LAT, which is applied together with two constraints: (1) nonnegativity constraint, (2) external geometry constraint, which can be utilized thanks to the 3D binary mask. The input to the second module are: the same projections and list of projection vectors that were passed to the first module, as well as 3D mask, which is the result of the first module.

## 3. Methodology of numerical study

To demonstrate generality of the GTVIC technique, we test its performance in combination with two competing LAT solvers characterized by different modeling assumptions. The comparison between the two algorithms, with and without the support of GTVIC, constitutes main part of the numerical analysis presented in Sections 4.1 and 4.2. Sections 4.3 and 4.4 report additional investigation which includes only the more accurate diffraction-wise solver.

The first LAT reconstruction technique under consideration is the Generalized Data Replenishment Algorithm (GDRA) [22]. Similarly to solutions based on the inverse Radon transform and all algebraic reconstruction techniques [25,26], it works under assumption of straight-line propagation of projection rays. The principle behind this method is to iteratively replenish the known part of the sinogram with the original sinogram. Unlike traditional DRA, however, GDRA operates on 3D data directly, allowing for completely arbitrary projection geometries. The implementation is based on the SIRT3D solver available in ASTRA Tomography Toolbox [27]. In each iteration of GDRA the standard non-negativity constraint (NNC) is enforced by default, which prohibits object RI values to fall under RI of the surrounding medium. Employment of GTVIC introduces an additional constraint generated by the first GTVIC module, which is the 3D mask, defining spatial limits on the reconstructed object volume. The second iterative solver is the classical Fourier inverse algorithm [28–30], here called the Fourier Transform Reconstruction Algorithm (FTRA), operating within the Rytov weak scattering approximation. The technique is based on direct implementation of the Fourier Diffraction theorem [26]. Following other ODT researchers [6,28,29,31], we implement a simple variant of FTRA, which relies on the nearest neighbor interpolation of projection spectra inside the image Fourier space. In case of small and well-focused objects recorded with excessive resolution, which is typical situation in ODT imaging of living cells, the mapping artifacts have little contribution to the overall reconstruction error. Again, the algorithm is tested in two variants, with the non-negativity condition only and with additional support of the GTVIC mask.

The numerical tests in Sections 4.1-4.4 were carried out on one synthetic data set generated from a biological cell phantom (Fig. 2(a)). The phantom is surrounded by water and consists of cytoplasm with varying RI distribution (from n = 1.359 to n = 1.371), a core with varying RI distribution (from n = 1.367 to n = 1.376) with mean RI value higher than the mean RI value of cytoplasm, and a vacuole with a constant RI value equal to RI value of water (n = 1.332). The RI distribution of the phantom was based on published papers that investigate RI values of biological cells [32–34]. The boundaries of the phantom were smoothed by rescaling with a 3x3x3 cubic interpolation kernel, so that there is no sharp edge between the cell and immersion. The structure of the phantom is presented in Fig. 2(a). Figure 2(b) depicts geometry of the virtual acquisition setup, corresponding to the system configuration typically used in LAT [10,14].

The sample is placed on a Petri dish fixed in horizontal position and illuminated by a set of plane waves from above; the scattered light is collected by a microscopic objective placed below. For this study the conical illumination pattern was chosen, with uniformly distributed projection directions inclined at angle *α* = 45° with respect to the optical axis of the system.

The synthetic sinogram of the phantom cell was prepared using a Fourier-based forward projector operating on Rytov fields. The projector essentially repeats main steps of the FTRA solver in reversed order, with additional object-domain padding for minimization of nonphysical overlapping between neighboring projection replicas. Consequently, the inherent inaccuracies of the FTRA solver related to the usage of the Rytov fields get cancelled out and are not manifested in the phantom reconstructions.

The related quality gain, however, applies to the GDRA solver as well, since weaker diffraction effects in projections translate to weaker image artifacts due to the straight line propagation approximation. This assures fair comparison conditions between the two solvers.

It also puts increased emphasis on reconstruction errors related to the incompleteness of the input data, which makes it easier to assess effectiveness of the additional GTVIC constraint. The following parameters were used in simulation of projection data: wavelength in vacuum *λ _{0}* = 0.633µm, objective numerical aperture

*NA*= 1.3, phantom cube voxel size

*dx*=

*dy*=

*dz*= 0.18µm, phantom grid size: 200 voxels, phantom grid size after upsampling (for

*K*-space extension to region occupied by projection Ewald spheres): 480 voxels,

*K*-space grid size after additional symmetric padding of the phantom: 720 voxels, projection size after symmetric clipping: 200 pixels. The calculations were performed in single precision.

## 4. Numerical tests

#### 4.1 GTVIC with generalized data replenishment algorithm

In this part we examine the reconstruction quality gain provided by GTVIC when applied to the diffraction-blind GDRA solver. The resulting RI images (normalized by subtracting RI of water), before and after addition of GTVIC, are shown in Figs. 3(a) and 3(d). Figures 3(b) and 3(e) present corresponding progress curves of the main GDRA loop. There are two measures of reconstruction accuracy used here, both relying on the original phantom as a reference: (a) total RMS calculated inside a cube containing the object, which reflects global reconstruction error of the 3D refractive index distribution and (b) the universal quality index (QI) [35] evaluated for the *XZ* plane only within the region occupied by the phantom (ignoring the immersion liquid volume), which puts emphasis on the internal structures of the object and ignoring misplacement of its borders. QI is a combination of three factors: loss of correlation between images, luminance distortion and contrast distortion, and can be expressed with Eq. (2):

*x*and

*y*,${\sigma}_{x}^{2}$ is sample variance and

*x*,

*y*are images compared.

The global RMS curve was used to define a slope-saturation condition terminating the tomographic solver loop; exactly the same condition was used consequently in all numerical tests. Also, for each reconstruction the accuracy of a 3D image has been assessed quantitatively by means of global QI curves presented in Figs. 3(c) and 3(f). The curves describe three separate cross-section planes rotated around *X*, *Y* and *Z* axis, correspondingly, so that in each case the full image volume is scanned. The global accuracy of a 3D image is expressed by the average value of the three curves *QI _{av}*.

The results reveal that the GTVIC-supported reconstruction exhibits accelerated convergence, which is an expected consequence of additional *a priori* information provided by the RI mask. The main GDRA loop acting with the non-negativity condition alone reaches the *RMS* saturation condition after 96 iterations (computation time: 77 minutes), whereas GDRA with GTVIC saturates after 9 iterations (computation time: 13 minutes). Additionally, the reconstruction quality is visibly improved with GTVIC, which is reflected by 2 percentage points increase in the global averaged QI measure. This comes in hand with reconstruction errors being uniformly distributed across the object volume, as follows from curves in Figs. 3(c) and 3(f).

#### 4.2 GTVIC with Fourier transform reconstruction algorithm

Figure 4 depicts numerical reconstructions obtained with the diffraction-wise FTRA solver. In Fig. 4(b) one can notice that this time the *QI(xz)* curve exhibits a local maximum before RMS reaches its saturation limit. This is due to a specific image distortion introduced by FTRA, which is visible in Fig. 4(b) as excessive elevation of RI values around the central horizontal plane. The effect is more pronounced in the first case (FTRA without GTVIC) due to larger number of FTRA iterations (27 iterations, computation time: 1.5 minute). In parallel to the increasing distortion, however, the positions of the top and bottom borders of the cell are being brought closer to their correct positions, which delays saturation of the global RMS measure. If the optimal number of iterations is exceeded, the RMS also exhibits a local extremum (coincides with the saturation point).

In the FTRA + GTVIC case, in which the number of iterations is reduced to only 5 (computation time: 7 minutes), the RI distribution distortion is significantly weakened. The quality improvement has two aspects: (a) sharpening the top and bottom boundaries of the object (trivial consequence of the GTVIC mask) and (b) enhancement of the refractive index structures close to these borders (nontrivial consequence of GTVIC), as a result of accelerated restoration of the object frequencies not covered directly by projection data. On the other hand, further examination of the RMS curve in Fig. 4(e) reveals a sharp local extremum when the saturation point (*n* = 5) is exceeded. This poses a problem in practical application of FTRA + GTVIC as the optimal iteration count may depend on the object and is not known in advance. A reliable formulation of the stop condition in experimental applications is yet to be developed. This issue is not present in the case of GDRA + GTVIC. Also, the slight RI distortion remaining in the *XZ* plane in Fig. 4(d), but not present in Fig. 3(d), suggests that FTRA is slightly less effective in restoration of the missing projection data than GDRA. This can be justified by fundamental differences in architecture between the two algorithms. Despite this weakness, it follows from charts in Fig. 3 and Fig. 4 that FTRA is always slightly more accurate than GDRA ( + 0.73% in averaged QI with GTVIC), as expected from a solver accounting for diffraction. Also, FTRA running on CPU (Intel® Core i7-3770, 3.40GHz × 8) is about 25 times faster than GDRA running on a graphic card (NVIDIA GeForce GTX 660Ti, 3GB, CUDA Driver v6.0). With the GTVIC strategy in work the overall computation time is reduced by factor 3.8 when GDRA is replaced by FTRA.

Further comparison of reconstructions calculated with GDRA and FTRA with and without GTVIC support is presented in Fig. 5, where 1-dimensional cross-sections of reconstruction from Figs. 3 and 4 are shown.

#### 4.3 Reconstruction from phase-only projections

It should be noted that GDRA, as any X-ray type tomographic solver (bound by the straight-line propagation approximation), takes real-valued data as input, while the Fourier-based FTRA operates on complex data. In the case of Optical Diffraction Tomography, in which the refractive index distribution of transparent objects is reconstructed, the real-valued projections encode unwrapped phase of the optical field, which is interpreted by the model as the integrated phase (along the projection direction). The complex-valued projections, on the other hand, represent both amplitude and phase of the optical field combined. When the Rytov approximation is used, both parts need to be split apart to enable proper phase unwrapping before the complex phase (Rytov field) is calculated.

Apart from these two cases, examined in previous sections, one can also consider a “middle way”, which is to pass phase-only data to the diffraction-wise solver. In case of a transition from an X-ray type solver this would mean a change in spectral representation of the projection data – from plane surfaces to spherical surfaces in the *K*-space. This could be a useful option, as the amplitude patterns are usually more sensitive to noise and other disturbances during real-world data acquisition. It is then of interest to examine how much quality improvement one should expect from introducing the amplitude patterns in ideal conditions.

Figure 6 depicts reconstruction result calculated by FTRA + GTVIC after the amplitude parts of the complex sinogram have been equalized (set to unity). Interestingly, when Fig. 6 is compared with Figs. 3(d) and 4(d), it appears that the amplitude patterns encode most of the object information carried by the scattered field. The change of the solver alone adds 0.29% to the averaged Quality Index, while inclusion of the amplitude patterns provides additional 0.44%. The second improvement step, between Figs. 4(d) and 6, is visible as slight enhancement in sharpness of the cytoplasm details across the hole phantom volume, combined with cancellation of the border artifact (resembling a cell membrane) in the horizontal cross-section plane.

#### 4.4 Reconstruction from few-view data

The numerical analysis described so far was based on a single sinogram composed of 180 projections in a conical illumination pattern. In this part we examine reconstruction quality by FTRA( + GTVIC) as a function of projection count in a projection subset taken out from the original, full sinogram. The results will point to optimal acquisition schemes for quick measurements in a real-world setup.

We start from a sinogram acquired in full analogy to the previous one (Section 3) but this time with total number of 360 projections. Then, we choose subsets with increasingly smaller sizes, keeping the angular pattern of illumination azimuthally as uniform as possible.

The plots depicted in Fig. 7 present the final quality of images obtained from FTRA with and without GTVIC support, expressed by the global QI values averaged over 3x21 test plane orientations, as in the lower plots in Figs. 3(c) and 3(f) and Figs. 4(c) and 4(f). Figure 7(a) corresponds to an idealized situation, in which it was assumed that the refractive index mask is always known precisely (calculated from all 360 projections). Despite this uniformity condition, the quality curve displays erratic jumps, even between the neighboring values of projection count, both with and without the support of GTVIC. A closer inspection of the iteration progress curves in the non-GTVIC case (with iteration counts above 20) revealed that the RMS-based stop condition could not be blamed for the observed variations, that is the lower points in Fig. 7(a) correspond to reconstructions, which actually saturated at lower quality levels. It has also been checked that the quality variations persist even if separate sinograms were prepared for individual projection counts, with illumination directions distributed on a cone with the same uniformity as in the full-size sinogram (360 views). Therefore, the observed quality variations can only be attributed to inherent sensitivity of the solver to very slight changes in the passed input data. It is in fact a basic property well recognized in broader field of iterative/adaptive inverse problem solvers [36]. Unless the algorithm is specifically designed to search for a global optimum, it will gravitate towards the closest of many local solutions.

Figure 7(b) depicts results obtained in a realistic scenario, where the refractive index mask is calculated from the limited projection set. In order to reduce the computational load, only 9 values of projection count were chosen, each of them in two alternative configurations. The secondary (dotted) curves represent projection subsets with identical illumination geometries, but shifted azimuthally by 1° with respect to the full bundle of 360 projection directions. We again observe significant variations of the quality levels, which can only be attributed to very subtle differences in the input data contents.

It also follows from Fig. 7(b) that best quality reconstructions with FTRA are possible for projection counts no larger than 60, which confirms known capabilities of this solver in limited angle and few-view configurations [13]. The GTVIC-enhanced variant of the algorithm seems to be slightly more restrictive, with about 120 projections required for results lying in the best quality range. This is partly due to lower relative depth of quality variations owed to the presence of the GTVIC mask (as follows clearly from Fig. 7(a), but also due to adaptive nature of the mask itself, which tends to fit the object tighter and more accurately as the projection count increases.

## 5. Experimental verification

In this part we examine the effectiveness of the GTVIC strategy when applied to limited angle ODT imaging of fixed and living cells. Since the superior performance of the Fourier inversion algorithm (FTRA) over the X-ray algorithm (GDRA) has been already demonstrated for this type of objects in Section 4, here only the FTRA solver is considered.

The object under study is a single, fixed C2C12 muscle cell immersed in water, measured in a fixed-object configuration, described in details in [10]. The measurement system was a digital holography-based tomographic phase microscope in Mach-Zehnder interferometer configuration. In this setup, the angular scanning of the sample (in order to obtain tomographic projections) is performed opto-electronically, using a spatial light modulator. Such system provides vibration-free scanning with a prospect for active phase correction. The data set consists of 360 projections acquired in a conical pattern (Fig. 2(b)) with the mean inclination angle *α* = 48°. Each projection direction has been independently calculated from orientation and frequency of the spatial carrier extracted from the corresponding interferogram. The model takes into account RI changes along the optical path between the sample and the objective lens, as well as magnification and nonlinear characteristics of the optical system [21]. The FTRA( + GTVIC) reconstructions were performed with the following parameters: wavelength in vacuum *λ _{0}* = 0.633µm, projection resolution 236x236, projection sample size

*dx*=

*dy*= 0.18µm, objective numerical aperture

*NA*= 1.25,

*K*-space grid size after zeropadding 354, image resolution 354x354x354, image sample size

*dx*=

*dy*=

*dz*= 0.12µm.

Figure 8 depicts a pair of spectral and spatial cross-sections of the images reconstructed with pure FTRA (including the nonnegativity constraint) and with FTRA additionally constrained by the refractive index GTVIC mask. The input image for mask generation, obtained by the Chambolle-Pock algorithm fed by 180 unwrapped-phase projections, was binarized and apodized in exactly the same way as for the synthetic data in Section 4. The iteration counts were set to 50 for FTRA and 5 for FTRA + GTVIC. The Chambolle-Pock algorithm for GTVIC mask optimization was run with 200 iterations. The total reconstruction times were 4 minutes for pure FTRA and 11 minutes for FTRA + GTVIC.

Comparison between Figs. 8(e) and 8(i), corresponding to a vertical cross-section plane, reveals that the FTRA + GTVIC procedure successfully restored all the structural details returned by FTRA (with no blur typical for edge-preserving reconstructions), while at the same time avoiding the high-frequency noise generated during 50 FTRA iterations (see Figs. 8(c) and 8(g)). Additionally, the cell borders are well defined with GTVIC, which alters the reconstruction in two ways: (1) the border geometry is more clear, which in this case helps to identify the cavities at the bottom side of the cell, and (2) the RI levels in the vicinity of the border are slightly elevated, which we attribute to the cancellation of the axial distortion characteristic for LAT images. This comes in pair with slight dumping of the RI levels in the central volume, which was also observed before on synthetic data (see Fig. 4 and Fig. 5).

Figures 8(f) and 8(j) depict a plane crossing the C2C12 cell near its base. In this view we observe a circular cavity, which is much better recovered with GTVIC. Here we can also observe very clearly the effect of refractive index elevation, expected in vicinity of the recovered bottom boundary surface.

## 6. Conclusion

The newly proposed Generalized Total Variation Iterative Constraint (GTVIC) is a meta-procedure designed to enhance Limited Angle Tomography (LAT) reconstructions by means of TV regularization when the object is not piecewise constant in the most part and would normally suffer from strong TV-related artifacts. In this study the effectiveness of GTVIC has been examined on the ground of Optical Diffraction Tomography of living biological cells. For the best demonstration of accuracy and versatility of the method, it has been tested in parallel with two different reconstruction engines: the non-diffractive (GDRA) one and the second operating within the first-order weak scattering approximation (FTRA). When applied to a synthetic data set, both algorithms required additional 6 minutes to optimize the GTVIC object mask. In the case of GDRA + GTVIC this represented 45% of the total computation time and provided 6-fold speedup relative to GDRA alone, thanks to accelerated convergence enforced by GTVIC. In the case of FTRA, which is by itself much faster than GDRA, the first GTVIC stage took 88% of total computation time, making the FTRA + GTVIC option 4.5 times slower compared to FTRA alone. The overall results suggest that the GTVIC strategy offers significant quality improvement in LAT reconstructions for a reasonable price in terms of computational effort. It can be conveniently combined with any algorithm for tomographic reconstruction that allows for arbitrary constraints within the image domain. However, to make the technique fully automated and reliable, some of its aspects need further development. These include: (a) optimal choice of the GTVIC mask parameters, which may depend on object characteristics, (b) optimal choice of number of iterations for the Chambolle-Pock algorithm during the GTVIC mask preparation stage and (c) a well-defined stop condition for the main solver loop, especially needed for FTRA + GTVIC, which tends to diverge quickly after the optimal number of iterations is reached.

## Acknowledgments

The research leading to the described results were carried out within the program TEAM/2011-7/7 of Foundation for Polish Science, co-financed from the European Funds of Regional Development. The authors would like to acknowledge the support from the statutory funds and the grants of the Dean of Mechatronics Faculty, Warsaw University of Technology. We would like to thank Kyoohyun Kim and his research team from Korea Advanced Institutes of Science and Technology, South Korea, for their kind support in implementation of our Fourier inversion solver, as well as to Folkert Bleichrodt and Tristan van Leeuwen from Centrum Wiskunde & Informatica, Netherlands, for providing the code for the Chambolle–Pock algorithm.

## References and links

**1. **T. E. Gureyev and K. A. Nugent, “Rapid quantitative phase imaging using the transport of intensity equation,” Opt. Commun. **133**(1-6), 339–346 (1997). [CrossRef]

**2. **C. J. Mann, P. R. Bingham, V. C. Paquit, and K. W. Tobin, “Quantitative phase imaging by three-wavelength digital holography,” Opt. Express **16**(13), 9753–9764 (2008). [CrossRef] [PubMed]

**3. **B. Rappaz, B. Breton, E. Shaffer, and G. Turcatti, “Digital holographic microscopy: a quantitative label-free microscopy technique for phenotypic screening,” Comb. Chem. High Throughput Screen. **17**(1), 80–88 (2014). [CrossRef] [PubMed]

**4. **G. Popescu, “Quantitative phase imaging of nanoscale cell structure and dynamics,” Methods Cell Biol. **90**, 87–115 (2008). [CrossRef] [PubMed]

**5. **H. Ding, Z. Wang, X. Liang, S. A. Boppart, K. Tangella, and G. Popescu, “Measuring the scattering parameters of tissues from quantitative phase imaging of thin slices,” Opt. Lett. **36**(12), 2281–2283 (2011). [CrossRef] [PubMed]

**6. **K. Kim, H. Yoon, M. Diez-Silva, M. Dao, R. R. Dasari, and Y. Park, “High-resolution three-dimensional imaging of red blood cells parasitized by Plasmodium falciparum and in situ hemozoin crystals using optical diffraction tomography,” J. Biomed. Opt. **19**(1), 011005 (2013). [CrossRef] [PubMed]

**7. **J. Kostencka, T. Kozacki, A. Kuś, M. Dudek, M. Kujawińska, and B. Kemper, “Holographic method for capillary induced aberration compensation for 3D tomographic measurements of living cells,” Proc. SPIE **8792**, 879204 (2013). [CrossRef]

**8. **J. Kostencka, T. Kozacki, A. Kuś, and M. Kujawińska, “Accurate approach to capillary-supported optical diffraction tomography,” Opt. Express **23**(6), 7908–7923 (2015). [CrossRef] [PubMed]

**9. **M. Kujawińska, W. Krauze, A. Kus, J. Kostencka, T. Kozacki, B. Kemper, and M. Dudek, “Problems and solutions in 3-D analysis of phase biological objects by optical diffraction tomography,” Int. J. Optomechatronics **8**(4), 357–372 (2014). [CrossRef]

**10. **A. Kus, W. Krauze, and M. Kujawinska, “Active limited-angle tomographic phase microscope,” J. Biomed. Opt. **20**(11), 111216 (2015). [CrossRef] [PubMed]

**11. **A. Kuś, W. Krauze, and M. Kujawińska, “Limited-angle, holographic tomography with optically controlled projection generation,” Proc. SPIE **9330**, 933007 (2015). [CrossRef]

**12. **A. Kuś, W. Krauze, M. Kujawińska, and M. Filipiak, “Limited-angle hybrid diffraction tomography for biological samples,” Proc. SPIE **9132**, 91320O (2014).

**13. **E. Y. Sidky, C.-M. Kao, and X. Pan, “Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,” J. Opt. Soc. Am. **25**, 1772–1782 (2009).

**14. **Y. Cotte, F. Toy, P. Jourdain, N. Pavillon, D. Boss, P. Magistretti, P. Marquet, and C. Depeursinge, “Marker-free phase nanoscopy,” Nat. Photonics **7**(2), 113–117 (2013). [CrossRef]

**15. **S. Vertu, I. Yamada, J.-J. Delaunay, and O. Haeberlé, “Tomographic observation of transparent objects under coherent illumination and reconstruction by filtered backprojection and Fourier diffraction theorem,” Proc. SPIE **6861**, 686103 (2008). [CrossRef]

**16. **T. Sato, S. J. Norton, M. Linzer, O. Ikeda, and M. Hirama, “Tomographic image reconstruction from limited projections using iterative revisions in image and transform spaces,” Appl. Opt. **20**(3), 395–399 (1981). [CrossRef] [PubMed]

**17. **D. Wang, H. Qiao, X. Song, Y. Fan, and D. Li, “Fluorescence molecular tomography using a two-step three-dimensional shape-based reconstruction with graphics processing unit acceleration,” Appl. Opt. **51**(36), 8731–8744 (2012). [CrossRef] [PubMed]

**18. **T. Humphries, A. Saad, A. Celler, G. Hamarneh, T. Möller, and M. Trummer, “Segmentation-based regularization of dynamic SPECT reconstruction,” in Nuclear Science Symposium Conference Record, (IEEE, 2009), 2849–2852.

**19. **D. Hyde, E. L. Miller, D. H. Brooks, and V. Ntziachristos, “Data specific spatially varying regularization for multimodal fluorescence molecular tomography,” IEEE Trans. Med. Imaging **29**(2), 365–374 (2010). [CrossRef] [PubMed]

**20. **X. Jin, L. Li, Z. Chen, L. Zhang, and Y. Xing, “Anisotropic total variation for limited-angle CT reconstruction,” in Nuclear Science Symposium Conference Record, (IEEE, 2010), 2232–2238. [CrossRef]

**21. **W. Krauze, A. Kuś, and M. Kujawinska, “Limited-angle hybrid optical diffraction tomography system with total-variation-minimization-based reconstruction,” Opt. Eng. **54**(5), 054104 (2015). [CrossRef]

**22. **W. Krauze, P. Makowski, and M. Kujawińska, “Total variation iterative constraint algorithm for limited-angle tomographic reconstruction of non-piecewise-constant structures,” Proc. SPIE **9526**, 95260Y (2015).

**23. **A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” J. Math. Imaging Vis. **40**(1), 120–145 (2011). [CrossRef]

**24. **N. Otsu, “A threshold selection method from gray-level histograms,” Automatica **11**, 23–27 (1975).

**25. **A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the art algorithm,” Ultrason. Imaging **6**(1), 81–94 (1984). [CrossRef] [PubMed]

**26. **A. C. Kak and M. Slaney, *Principles of Computerized Tomographic Imaging* (Society of Industrial and Applied Mathematics, 2001).

**27. **W. van Aarle, W. J. Palenstijn, J. De Beenhouwer, T. Altantzis, S. Bals, K. J. Batenburg, and J. Sijbers, “The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography,” Ultramicroscopy **157**, 35–47 (2015). [CrossRef] [PubMed]

**28. **Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express **17**(1), 266–277 (2009). [CrossRef] [PubMed]

**29. **K. Kim, K. S. Kim, H. Park, J. C. Ye, and Y. Park, “Real-time visualization of 3-D dynamic microscopic objects using optical diffraction tomography,” Opt. Express **21**(26), 32269–32278 (2013). [CrossRef] [PubMed]

**30. **B. P. Medoff, W. R. Brody, M. Nassi, and A. Macovski, “Iterative convolution backprojection algorithms for image reconstruction from limited data,” J. Opt. Soc. Am. **73**(11), 1493–1500 (1983). [CrossRef]

**31. **Y. Kim, H. Shim, K. Kim, H. Park, S. Jang, and Y. Park, “Profiling individual human red blood cells using common-path diffraction optical tomography,” Sci. Rep. **4**, 6659 (2014). [CrossRef] [PubMed]

**32. **B. Kemper, S. Kosmeier, P. Langehanenberg, G. von Bally, I. Bredebusch, W. Domschke, and J. Schnekenburger, “Integral refractive index determination of living suspension cells by multifocus digital holographic phase contrast microscopy,” J. Biomed. Opt. **12**(5), 054009 (2007). [CrossRef] [PubMed]

**33. **B. Kemper, L. Schmidt, S. Przibilla, C. Rommel, A. Vollmer, S. Ketelhut, J. Schnekenburger, and G. von Bally, “Influence of sample preparation and identification of subcelluar structures in quantitative holographic phase contrast microscopy,” Proc. SPIE **7715**, 771504 (2010). [CrossRef]

**34. **C. Fang-Yen, W. Choi, Y. Sung, C. J. Holbrow, R. R. Dasari, and M. S. Feld, “Video-rate tomographic phase microscopy,” J. Biomed. Opt. **16**(1), 011005 (2011). [CrossRef] [PubMed]

**35. **Z. Wang and A. C. Bovik, “A universal image quality index,” IEEE Signal Process. Lett. **9**(3), 81–84 (2002). [CrossRef]

**36. **P. C. Hansen, *Discrete Inverse Problems: Insight and Algorithms* (Society of Industrial and Applied Mathematics, 2010).