3D scattering microphantom sample to assess quantitative accuracy in tomographic phase microscopy techniques

In this paper we present a structurally-complex biomimetic scattering structure, fabricated with two-photon polymerization, and utilize this object in order to benchmark a computational imaging system. The phantom allows to tailor the scattering by modifying its degrees of freedom i.e. refractive index contrast and scattering layer dimensions and incorporates a 3D imaging quality test, representing a single cell within tissue. While the sample may be used with multiple 3D microscopy techniques, we demonstrate the impact of scattering on three tomographic phase microscopy (TPM) reconstruction methods. One of these methods assumes the sample to be weak-scattering, while the other two take multiple scattering into account. The study is performed at two wavelengths (visible and near-infrared), which serve as a scaling factor for the scattering phenomenon. We find that changing the wavelength from visible into near-infrared impacts the applicability of TPM reconstruction methods. As a result of reduced scattering in near-infrared region, the multiple-scattering-oriented techniques perform in fact worse than a method aimed for weak-scattering samples. This implies a necessity of selecting proper approach depending on sample’s scattering characteristics even in case of subtle changes in the object-light interaction.

One of the modern challenges in computational optics is to image scattering samples with high resolution 1 . This can be attributed to the fact that complex biological structures such as spheroids or organoids tend to be more relevant models than 2D cell cultures e.g. for drug discovery 2 . Also, most in vivo imaging techniques require the probing light to pass through the complex structure of a tissue which highly limits imaging depths due to multiple scattering. This demand stimulates the development of new methods 1,3-5 , however, it is difficult to select an appropriate one based on the scattering strength of the analyzed specimen. For this reason a versatile, repeatable and quantitative method for the evaluation of different imaging systems and algorithms is essential to determine their limits of applicability depending on the object's scattering properties. One possibility is to use calibrated microphantoms as imaging targets. Unfortunately, existing microphantoms are typically either weak-scattering (e.g. index-matched microspheres) or overly simplistic (e.g. index-mismatched microspheres) 6,7 compared to the types of heterogenously scattering multicellular samples that multiple-scattering methods are intended for. This is a critical limitation when characterizing computational imaging methods that utilize nonconvex solvers, where iterative convergence depends on the complexity of the energy landscape and directly associates with a sample's 3D complexity 8 .
In this work we present a 3D-printed microphantom with multiple-scattering refractive index (RI) distribution. To do so, we leverage recent developments in 3D printing via direct laser writing [9][10][11][12] . Among multiple available 3D printing techniques [13][14][15][16][17][18] , we chose a two-photon polymerization that enables 3D printing of microphantom samples with known geometry and calibrated RI. When compared to other implementations of direct laser writing, it allows (1) to control the RI with relatively high modulation range, (2) to adjust the RI contrast or scattering strength post-fabrication using different immersion liquid and (3) to handle and measure the microphantom in the same way as biological specimens. Next, we present the application of the phantom in the field of tomographic phase microscopy (TPM), a technique which has demonstrated impressive biological www.nature.com/scientificreports/ imaging results in prior works. However, it is important to note that all computational imaging methods can be evaluated with the proposed procedure. TPM is a quantitative, label-free imaging method that utilizes optical projections through a semi-transparent sample along various illumination angles to reconstruct the sample's 3D RI. This method has found several applications in biological imaging, where RI is directly related to the dry mass distribution at the cellular and subcellular levels. Refractive index and dry mass are known to be crucial factors in analyzing the current stage of cell cycle 19 , cell structure 20,21 , photobiochemical effects on cells 22 , influence of external factors on cellular parameters 23,24 and many others. Given the wealth of information provided by the dry mass at the single-cell level, there is significant demand to extend the capability to analyze dry mass to large multicellular clusters, thick tissue slices, or whole microorganisms. However, to reconstruct 3D RI, traditional TPM methods utilize critical assumptions in their computational reconstruction methodologies that rely on the sample being weakly scattering 25 . These assumptions limit samples to having thicknesses on the order of only tens of microns. For thick, complex samples, reconstruction frameworks that accommodate for multiple scattering must be utilized. To this end, numerous TPM approaches have been proposed in recent years that introduce new frameworks to accommodate multiple scattering [26][27][28][29][30][31] . Notably, these approaches utilize nonlinear and nonconvex solvers to iteratively solve for a sample's 3D RI. Though these methods have demonstrated impressive results in reconstructing RI in multiple-scattering samples, their quantitative accuracy has not yet been robustly characterized experimentally, and the presented results usually do not allow comparison of different methods in order to select proper approach for a given scattering level in a sample. The general strategy to experimentally evaluate quantitative accuracy in TPM methods is to reconstruct 3D RI in samples with known RI distributions 6,7,32 . A multiple-scattering TPM method that outputs accurate 3D RI reconstructions of a weakly scattering or overly simplistic microphantom cannot be expected to output similarly accurate RI reconstructions for more complex multiple-scattering samples, where the probability of converging to local minima is drastically higher. To robustly characterize the quantitative accuracy of multiple-scattering TPM methods, it is imperative to use microphantoms with known 3D RI that mimic the structural complexity of the types of samples that the TPM methods are intended to image. To the best of our knowledge, these types of gold-standard multiple-scattering microphantoms do not exist.

Results
In this section we present the design of the 3D printed microphantom and its application in evaluation of three TPM reconstruction methods.
3D printed microphantom. The scattering microphantom we designed and 3D printed consists of a celllike target with internal test structures 33 embedded within a pseudo-random distribution of rods that switch their orientation across various layers (Fig. 1a). The width and height of each rod is equal to 0.5 µ m and 1.8 µ m respectively. The lateral distance between the rods in each layer is randomized between 0.7 µ m to 3 µ m and the layers are stacked vertically every 1.4 µ m. The resulting structure is transparent (over 99% transmittance for the www.nature.com/scientificreports/ extinction coefficient of 0.1 mm −1 34 ) and multiple scattering 35 . The final scattering region is a 60 µ m × 60 µ m × 40 µ m cube with a fill-factor of roughly 25% (fraction of volume occupied by the polymer). The cell-like target embedded within the layers of randomly distributed rods mimic a single biological cell encased within a scattering volume, which represents e.g. tissue. The cell-target comprises of substructures that enable assessment of quantitative accuracy in TPM imaging systems 36 . Information about 3D printing procedure are given in Sect. "3D scattering microphantom fabrication". The main features within the cell target include resolution test targets, nucleoli suspended in a nucleus, and a region of slow RI variation (see Fig. 1d). Notably, the resolution test target comprises of lines with increasing spatial frequency 37 up to 1667 lp/mm. By assessing the maximum spatial frequency of lines that can be distinguished within the cell test target, the imaging resolution of a TPM method of choice can be characterized and potentially compared to different methods of assessing resolution described by other research groups 38,39 . The pseudo-random distribution of rods that compose the scattering portion of the whole microphantom are suppressed within 0.5 µ m of the cell target and do not intersect with any of the test structures. A model of the phantom is available in Dataset 1 40 .
TPM systems evaluation. Three TPM reconstruction methods were implemented and evaluated through comparison of tomographic reconstructions of the proposed microphantom. The methods are: (1) Gerchberg-Paopulis with support constraint (GPSC) 41 , (2) multi-slice beam-propagation with electric-field measurements (MSBP-E) 42 and (3) multi-slice beam-propagation with intensity-only measurements (MSBP-I) 31 . In order to perform the comparison, the phantom has been measured with the TPM device and two datasets were composed from the complex-valued scattered field measurements of the microphantom being illuminated at different angles, using both 633 nm and 835nm wavelengths. For GPSC and MSBP-E, these datasets were used directly. For MSBP-I, only the amplitude components were used. The measurements and their corresponding 3D RI reconstructions are available in Dataset 1 40 . Figure 2 shows the reconstruction results for both wavelengths. We characterize lateral (x and y) resolution of the TPM reconstructions by visualizing the resolution test lines inside the microphantom. To quantitatively compare the three TPM reconstruction methods, horizontal and vertical cross-sectional plots across these resolution tests were generated by computing the average and standard-deviation of the pixel-values across rows or columns adjacent to the dashed white a-a and b-b lines, by ±4 pixels. These cross-sectional plots are shown below.

Discussion
The known complex geometry and RI distribution of the developed microphantom allows to show that changing the wavelength from visible into near-infrared impacts the applicability of TPM reconstruction methods 43,44 . The magnified views of the resolution test region after reconstructing the microphantom using the GPSC algorithm with 633 nm wavelength light reveals significant grainy artifacts which occlude the line features within the microphantom's test region. We note that these artifacts are heavily decreased when reconstructing with 835 nm wavelength light. This suggests that, at 633 nm wavelength, the microphantom is too highly-scattering for GPSC to be applicable. The decreased GPSC reconstruction artifacts for 835 nm matches conventional knowledge that longer wavelengths of light are more resistant to scattering than shorter wavelengths. In order to confirm that this observation is due to reduced scattering and not different noise characteristics between the two light sources, we analyzed the standard deviation of phase-noise in an object-free region within the input datasets. For 633 nm wavelength light, a phase-noise standard deviation of σ = 0.10 radians was observed, while for 835 nm wavelength light, a phase-noise standard deviation of σ = 0.08 radians was observed. Given such a small variation between these noise characteristics, we conclude that the major factor in the GPSC reconstruction quality that affects its ability to visualize the microphantom's test lines is the scattering strength of the microphantom at the two different wavelengths.
Notably, both MSBP-E and MSBP-I use total-variation (TV) regularization in order to stabilize the convergence of the nonconvex iterative solver in the presence of noise 45 . Especially in the case of using 633 nm wavelength light, TV regularization results in 3D RI reconstructions with less noise compared to the 3D reconstructions computed via GPSC, which does not use TV regularization. This can be directly visualized in the 2D cross-sections in Fig. 2, and is confirmed by the bounds of standard deviation shown in the 1D cross-sectional plots. However, the drawback of regularization is that it has a blurring effect on high-resolution features. Because the microphantom is less scattering in 835 nm wavelength light, only GPSC managed to reconstruct the high spatial-frequency test lines within the microphantom. Typically, the strength of TV regularization is manually tuned to fit experimental factors and balance the tradeoff between achieving iterative stability versus high imaging resolution.
In terms of the average RI, the knowledge about the ground truth RI distribution of the phantom makes it possible to quantitatively compare reconstruction results with GPSC, MSBP-E, and MSBP-I. We see that all methods successfully capture the bulk characteristics of the microphantom. As described above, the 3D reconstruction via GPSC exhibits grainy artifacts when using 633 nm wavelength light, likely due to the microphantom being multiple-scattering at that wavelength. Furthermore, MSBP-I outputs slightly overestimated RI values, and also suffers from low-frequency spatial artifacts (which have been observed in other intensity-only phase-imaging techniques 46,47 ). Other works have shown that MSBP-I demonstrates higher accuracy when using partiallycoherent illumination, which drastically reduces coherent noise 31 . Future work may include repeating this analysis across a larger range of TPM reconstruction techniques with more complex scattering microphantoms.
With the presented results we show that changing the illumination wavelength affects the scattering nature of the microphantom. Specifically, though the microphantom is multiple scattering with 633 nm wavelength light, it is weak scattering with 835 nm wavelength light. This naturally indicates that the optimum choice for illumination wavelength must balance between resolution ( NA for single projection) and scattering strength. As has been www.nature.com/scientificreports/ shown, there are cases when instead of applying multiple-scattering methods, it is advantageous to increase the illumination wavelength (thus decreasing the scattering strength of the sample) and apply a method based on the Fourier Diffraction Theorem that does not utilize the total-variation constraint. More fundamentally however, we showed that 3D printed microphantoms enable quantitative assessment of 3D RI reconstruction accuracy across various TPM methodologies. This capability is important when choosing a TPM method optimized for specific classes of samples and imaging conditions. Future works will focus on developing methods to quantify scattering strengths of phantoms and relate these quantities with the scattering properties of different tissue types. If successful, this would enable us to design and fabricate (using the presented methods) 3D microphantom structures to mimic a wide range of biological specimens ranging from multicellular clusters to bulk tissues and small organisms. Another possible direction is to tune scattering parameters in the phantom to characterize imaging performance for various techniques used to image into scattering tissue, such as optical coherence tomography or confocal reflectance microscopy. Finally, we envision that one can exploit the flexibility of two-photon polymerization to fabricate microphantoms on different substrates (e.g. at the end-face of the optical fibre for sample rotation tomography), using biocompatible www.nature.com/scientificreports/ resins (to combine test targets with the living cells in the single measurement volume 48 ) or modify the resin with functional particles 18 (e.g. to accommodate systems that also measure absorption [49][50][51][52][53] ).

Materials and methods
We present below the (1) methodology with which we 3D-print multiple scattering microphantoms; (2) the optical design of the TPM imaging systems that we use to experimentally collect scattering electric-field measurements of the microphantom; and (3) short theoretical descriptions of three tomographic algorithms that were used to reconstruct 3D RI from the measured data.
3D scattering microphantom fabrication. The phantom is fabricated using two-photon laser lithography, in which a focused laser beam is scanned within liquid resin. The resin within the laser's focal volume is locally polymerized. Adjusting the scanning trajectory and the exposure time of the laser beam enables simultaneous control over the 3D printed geometry (accuracy at the order of 100 nm) and RI (accuracy at the order of 5 × 10 −4 , maximal RI = 0.03 within the structure) in three dimensions. We used Photonic Professional GT (Nanoscribe GmbH) equipped with a 1.3 NA 100× microscope objective and piezo scanning stage. The phantom is fabricated in the IP-Dip resin (Nanoscribe GmbH) on top of a #1.5H coverslip (dip-in configuration 54 ). After fabrication the structure was developed in PGMEA (Propylene glycol monomethyl ether acetate; 12 min), followed by isopropyl alcohol (10 min) and then blow-dried. The full methodology for fabrication and validation of the features can be found in our previous work 33 .
To conduct our TPM imaging experiments, the microphantom was immersed in Zeiss Immersol 518F oil (RI 632 nm = 1.5123), which provides similar RI contrast as in the case of cells immersed in culture medium. By using immersion oils with varying RI, it is possible to adjust the scattering properties of the microphantom post-fabrication.
Measurement system. In this work, an optical system, as shown in Fig. 3a) was used in order to study the scattering phantom. The system is a Mach-Zehnder-based TPM microscope 55 , working in a limited-angle configuration with stationary sample and illumination rotated with a galvo mirror (Thorlabs GVS212/M) 56 . The research was performed with two wavelengths and thus there were two modified versions of the presented tomographic microscope. First version, TPM 633 works with wavelength = 633 nm and the second, TPM 835 with = 835 nm. The input beam (S in Fig.3a) is delivered with an optical fiber, collimated and then split into object and reference arm. In the TPM 633 system the light source was a volume Bragg grating laser (Necsel NovaTru Chroma 633 SLM), S 633 , providing a single longitudinal mode and offering long coherence length. The TPM 835 system utilized a swept source (Superlum Broadsweeper BS-840-2-HP, = 800-870 nm), S 835 set at = 835 nm. Due to difference in coherence length, an additional delay module was placed in the reference beam for TPM 835 measurements. The beam-splitting cubes in this work were either coated for 400-700 nm or 700-1100 nm depending on the wavelength used. The focal length of the tube lens TL1 was EFL 633 = 150 mm and EFL 835 = 200 mm respectively. Both microscope objectives (MO1 and MO2) in Fig. 3 were 100× NA 1.  www.nature.com/scientificreports/ system was a CMOS sensor in both cases, with 3.45 µm pixel size (JAI BM500GE) in case of CAM 633 and 5.5µm pixel size (Basler acA2040-180km) in case of CAM 835 . The minimum magnification, which is imposed by pixel size and wavelength in order to assure correct hologram recording for each projection 57 is M 633min = − 44.2 and M 835min = − 53.5 , which is satisfied in both cases. A sample hologram is presented in Fig. 3b). Both systems were set to illuminate the sample with a circular scanning scenario (see Fig. 3d) at zenith angle θ = 47 • and provided 180 projections spaced at ϕ = 2 • . A sample of the phase and amplitude provided by the system at 835nm is presented in Fig. 3c and e.
Reconstruction algorithms. There were three different TPM reconstruction methods for computing 3D RI of the microphantom used for comparison. The data for each method was acquired through the multi-angle scattering measurements captured as described in Sect. "Measurement system". We provide a short description of these methods below. Complete descriptions of these methods is given in respective references.
Gerchberg-Papoulis with object support. To provide a baseline standard to compare to TPM reconstruction algorithms utilizing multiple-scattering models, we first reconstruct the microphantom's 3D RI using a weakscattering TPM method. We specifically use the Gerchberg-Papoulis algorithm enhanced with additional finite object support constraint (GPSC) 41 . Complex-valued electric-fields measured by our TPM systems are used as inputs. The procedure is performed in two steps. First, an initial tomographic 3D RI distribution is reconstructed from the electric-field measurements with strong total-variation regularization 45 . This is performed through the Chambolle-Pock 58 optimization method and implemented with ASTRA tomography toolbox 59 . The result undergoes binarization and a finite object support is generated. Secondly, a classic Gerchberg-Papoulis algorithm is used which is an iterative version of Direct Inversion method (also known as the Wolf transform) 60 . This iterative procedure is based on the Fourier Diffraction Theorem 61 and utilizes first-order scattering approximation. Here, the reconstruction and its Fourier transform are calculated alternately and constraints are applied: nonnegativity and finite object support in the signal domain, and replenishment of original projections in the frequency domain.

Multi-slice beam-propagation with electric-field measurements.
Our main method to model multiple-scattering is the multi-slice beam-propagation method (MSBP) 62 , which has recently shown promising results for biological imaging 42,63 . In our first implementation of MSBP, we use the same exact electric-field dataset utilized by GPSC from above. An initial guess of the sample's 3D RI is selected to start off the iterative procedure. Afterwards, the MSBP method is used to simulate scattering measurements resulting from plane waves propagating through the sample's initial estimated 3D RI. The scattered fields resulting from this simulation are compared with those obtained experimentally with our TPM systems. The error computed between the simulated and experimental measurements is back-propagated through each layer of the 3D sample estimate to incrementally modify the RI value of each voxel. Continued iterations repeating these steps eventually result in the 3D sample estimate converging to a stable steady-state solution. We implemented the MSBP with electric-field measurements (MSBP-E) through the Learning Tomography algorithm (LT) 42 . The LT procedure is an iterative optimization algorithm with additional weak total-variation (TV) regularization applied in each iteration to ensure convergence. We found that the method works best when an initial-guess is chosen as a starting point for the iterative process. In this paper, we use the Direct Inversion method to provide the initial guess.

Multi-slice beam-propagation with intensity-only measurements.
Recent works have demonstrated that the gradient-update step within the MSBP method can be reformulated to reconstruct 3D RI from only non-interferometric intensity measurements 31 . The key advantages of this method include the use of a non-interferometric imaging system, which are resistant to mechanical instabilities that often limit long-term use of dual-arm interferometers without realignment. Furthermore, the light source can be partially coherent, to avoid coherent speckle artifacts in the measurements, while retaining sufficient coherence necessary for RI reconstruction. For the purposes of demonstrating 3D RI reconstruction using this intensity-only variant of MSBP (which we refer to here as MSBP-I), we simply use the amplitude component of the electric-field measurements used for the GPSC and MSBP-E reconstructions, described above. Similarly to MSBP-E, total-variation regularization is applied at every iteration. The starting point for MSBP-I is a matrix of zeros.
Stopping criterion. All described TPM reconstruction methods are iterative procedures that use the same stopping criterion to automatically terminate the computations. This criterion is a modification of a method presented earlier 64 , and is described in Alg. 1 below. The general intuition behind this procedure is to terminate the iterative computation process when the dynamics of the change between 3D sample estimates outputted from consecutive iterations drops below a certain saturation level ǫ . In order to be less dependent on outliers, the median value of the dynamics from the last 10 iterations is calculated. The value for ǫ is chosen empirically for each algorithm. For GPSC ǫ = 0.02 , for LT and MS ǫ = 0.01 . Values of ǫ were empirically chosen to balance between incomplete convergence and reconstruction speed.