Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Optical diffraction tomography of 3D microstructures using a low coherence source

Open Access Open Access

Abstract

Optical diffraction tomography (ODT) is a label-free technique for three dimensional imaging of micron-sized objects. Coherence and limited sampling of 3D Fourier space are often responsible for the appearance of artifacts. Here we present an ODT microscope that uses low temporal coherence light and spatial light modulators to retrieve reliable 3D maps of the refractive index. A common-path interferometer, based on a spatial light modulator, measures the complex fields transmitted by a sample. Measured fields, acquired while scanning the illumination direction using a digital micro-mirror device, are fed into a Rytov reconstruction algorithm to obtain refractive index maps whose accuracy is directly evaluated on microfabricated 3D test objects. Even for challenging shapes such as pyramids, bridges, and dumbbells, we obtain volumetric reconstructions that compare very well with electron microscopy images.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Optical tomography provides a label-free solution for three dimensional imaging of micron-sized objects with diffraction limit resolution. To produce a 3D map of the sample refractive index, in all the diverse approaches reported in the literature, a set of 2D images is acquired while some sort of scanning is performed. A simple and effective strategy is to illuminate the sample with a low coherence light source using a high numerical aperture condenser and collect a stack of intensity patterns by scanning the image plane along the optical axis [13]. In this case, the reconstruction algorithm analyzes a three dimensional intensity pattern whose fluctuations, with respect to a constant background, are assumed to be linearly dependent on the refractive index. Assuming that the response of the optical system to a point scatterer is known, then within the single scattering approximation, i.e., the Born approximation [4,5], the volumetric intensity produced by the entire sample could be expressed as the superposition of point contributions distributed over the sample’s volume and its refractive index map can be retrieved by a deconvolution algorithm.

A first method to solve the inverse problem consists in minimizing the difference between the experimental data and a dataset computed using a forward numerical model [69]. This optimization procedure converges at the cost of large computational time. Alternatively, using the beam-propagation method as a forward numerical model, it is possible to compute the derivative of the transmitted field with respect to the distribution of the refractive index and therefore minimize a cost function with a steepest-descent algorithm [10].

Apart from these methods, which are based on the minimization of a cost function, faster inversion algorithms are reported in the literature. For instance, quantitative phase maps of cells, acquired with a Mach-Zehnder heterodyne interferometer, can be analyzed using a Radon transform based algorithm to obtain the refractive index in 3D [11,12]. This approach, which neglects diffraction, is suited for thin objects with low contrast that are not far from the objective focal plane. In this respect, Fourier diffraction tomography gives a more accurate reconstruction [1320]. For thin and low contrast objects, multiple scattering can be neglected and therefore the scalar wave equation can be approximated to find a simple relation between the 2D measured scattered field and the 3D refractive index in the Fourier space [4]. Along this vein, one can assume slowly varying perturbations on the phase gradient [4,5,21,22], i.e., the Rytov approximation, to find a simple solution to the scalar wave equation. Using this framework, one can solve the inverse problem with no severe restriction on the observed objects size and contrast [23] and avoiding time-consuming direct-search algorithms.

In this work, we use Rytov reconstruction algorithm in combination with a LED light source whose tilt angle can be scanned. Compared to lasers, light sources with low temporal coherence result in a substantial improvement in terms of speckle background and image stability so that reconstruction artifacts are more easily avoided [24]. More specifically, previous work took advantage of the lower temporal coherence of LED to reduce the speckle noise [2,3,7,9,15,22]. To retrieve the transmitted complex field, which is not a straightforward task given the low temporal coherence of the source, we opted for using a common-path interferometer based on a Spatial Light Modulator (SLM). SLMs have proven to be useful tools for microscopy, especially with regards to phase contrast imaging [25] and quantitative phase imaging [2628]. In preexisting quantitative phase imaging based on a differential interference microscope [2931] the phase gradient is measured along a single direction. Differently, here the SLM duplicates and laterally translates one of the two copies of the field and, upon phase-shifting it, it is possible to retrieve the phase gradient [31] on both axes. This allows for a reliable measurement of the phase map passing through the computation of its Laplacian. Having access to the field phases, we obtain the complex fields that are fed to the tomographic reconstruction algorithm to compute the 3D refractive index maps of biological and synthetic micro-fabricated samples.

2. Methods

2.1 Experimental setup and field measurement

Our experimental setup is sketched in Fig. 1(a). We use a green LED (Thorlabs M530L4, $\lambda =530\pm 10$ nm) as light source. The lens L1 creates a magnified image of the LED on a digital micromirror device (DMD, Texas Instruments DLP3000), whose active surface is completely illuminated. Light impinging on the DMD is oriented in such a way that light reflecting from the pixels in the “on” state is collected by the lens L2, while light reflecting from pixels in the “off” state is blocked. A magnified image of the DMD display is reproduced on the back-focal plane of the condenser lens L3 ($\text {NA}=0.8$). The magnification is adjusted so that the DMD minor axis (3.4 mm) fits into the 1 inch diameter of the lens. On the sample plane, to each pixel in the “on” state it corresponds a collimated illumination beam with a given incidence. During the acquisition, to have enough light, we turn on a group of pixels whose beam’s average incidence angle is pre-calibrated as explained later in the text. The size of the “on” pixels sub-matrix controls the spatial coherence of each impinging beam since it acts as a diaphragm regulating the numerical aperture of the illumination.

 figure: Fig. 1.

Fig. 1. (a) Experimental setup. The sample is illuminated by a low temporal coherence source. The illumination beam’s angle and aperture are selected using a DMD whose image is reformed by the lenses L1 and L2 on the back-focal plane of the condenser L3. Light is polarized at 45$^\circ$ with respect to the SLM active axis by the first polarizer P. On the camera plane, the polarization component along the SLM active axis is laterally shifted by $\Delta$ as a consequence of a linear grating displayed on the SLM. A second polarization analyzer projects the two polarization components on a 45$^\circ$ axis so that the two interfere. (b) Given a shift direction ${\bf \Delta }$, the SLM changes the relative phase of the two fields so that the gradient of the phase can be computed as in Eq. (3). The phase gradients computed with ${\bf \Delta }=(\pm \Delta,0)$ are combined to obtain a cleaner gradient along $x$; the same is done for ${\bf \Delta }=(0,\pm \Delta )$. The two phase gradients $\partial _x\theta$ and $\partial _y\theta$ are differentiated numerically to compute the Laplacian that is integrated to finally obtain $\theta$. Scalebar is 2 µm.

Download Full Size | PDF

For each illumination angle, our goal is to measure the transmitted complex field. Given the low temporal coherence, to detect the field phase, we must opt for common-path interferometry. The phase of a field is measured by duplicating the beam, e.g., using a diffraction grating [14,16], and, before the two duplicates are recombined, one is spatially filtered by a pinhole [14,16] or phase shifted [25,26,32]. A practical option, which does not rely on accurate alignment, is to measure the phase gradient of the beam by adding a programmable phase retarder to a Differential Contrast Microscope [29,30]; this method has also been referred to as Gradient Light Interference Microscopy (GLIM) [31]. In analogy to GLIM, here we measure the phase gradients, but both the beam duplication and the phase shifting are obtained through a single SLM (Meadowlark Optics HSP256–0532). Using the analyzer P placed right after the tube lens (see Fig. 1(a)), light is polarized along a direction aligned at 45$^\circ$ with respect to the SLM active axis (in our case the vertical axis). On the CCD plane, a translation ${\bf \Delta }$ of the vertical component with respect to the horizontal component is obtained by applying a linear phase mask on the SLM that is placed in the Fourier plane:

$$\Phi(\textbf{r})=\frac{2\pi}{\lambda f}\boldsymbol{\Delta}\cdot\textbf{r}+\phi$$
where $\textbf{r}$ are the 2D coordinates of the SLM pixels, $f=100$ mm is the focal length of the lens L5 that Fourier transforms the field from the SLM plane to the camera plane, and finally $\phi$ is a constant phase whose role will be described shortly hereafter. Both polarizations are projected again along an axis rotated by 45$^\circ$ with respect to the vertical direction, therefore, the resulting recorded intensity is:
$$I_\text{camera}^{\phi}(\textbf{r})=I(\textbf{r})+I(\textbf{r}+\boldsymbol{\Delta})+2\sqrt{I(\textbf{r})I(\textbf{r}+\boldsymbol{\Delta})}\cos\left(\theta(\textbf{r}+\boldsymbol{\Delta})-\theta(\textbf{r})+\phi\right)$$
where the intensity of the transmitted field $E$ is $I$ and $\theta$ is its phase. When $\Delta$ is small, the phase difference in the cosine can be approximated with the phase gradient; this can be extracted using four images, acquired while varying $\phi$ among the values $[0,\pi /2,\pi,3\pi /2]$, as follows [29]:
$$\boldsymbol{\nabla}\theta\cdot\boldsymbol{\Delta}=\arctan\left( \frac{I_\text{camera}^{\pi/2}-I_\text{camera}^{3\pi/2}}{I_\text{camera}^{0}-I_\text{camera}^{\pi}}\right).$$

If the SLM phase mask is oriented along the $x$ axis then $\boldsymbol{\Delta }=(\Delta,0)$ and the gradient $\partial \theta /\partial x$ is obtained. Given the phase mask in Eq. (1), the translation $\Delta$ is measured for the LED central wavelength, with small imaging bias due to dispersion. To minimize artifacts, $\partial \theta /\partial x$ is computed by averaging the two gradients obtained with $\boldsymbol{\Delta }=(\Delta,0)$ and $\boldsymbol{\Delta }=(-\Delta,0)$. Similarly, $\partial \theta /\partial y$ is obtained using both values of $\boldsymbol{\Delta }=(0,\pm \Delta )$. The numerical gradients $\partial \theta /\partial x$ and $\partial \theta /\partial y$ are differentiated numerically to obtain the phase Laplacian $\nabla ^2\theta$ that can be in turn integrated as $\theta =F^{-1}\left [ F\left [\nabla ^2\theta \right ]/k^2 \right ]$, where $F[\cdots ]$ indicates the 2D Fourier operator and $\textbf{k}$ is the coordinate in the reciprocal space. From Eq. (2) we can also extract the field intensity by simply averaging $I_\text {camera}^{\phi }$ in Eq. (2) for $\phi$ and $\boldsymbol{\Delta }$ and approximating $I(\textbf{r}+\boldsymbol{\Delta })\approx I(\textbf{r})$. The entire procedure is summarized in Fig. 1(b).

To compute the scattered field, the incident field must be previously recorded. This can be done by measuring a background field with an empty field of view:

$$U_B=\left(\sqrt{I(\textbf{r})} e^{i\theta(\textbf{r})} - \sqrt{ I_\text{bg}(\textbf{r})}e^{i\theta_\text{bg}(\textbf{r})}\right)e^{{-}i\theta_\text{bg}(\textbf{r})+i\textbf{q}\cdot\textbf{r}}$$
where the subscript bg indicates the intensity and the phase acquired with an empty field of view and $\textbf{q}$ is the wavevector of the incident field projected on the $x-y$ plane. The term outside the parentheses in Eq. (4) is needed to improve the image quality: theoretically the background phase $\theta _\text {bg}$ should be simply a linear phase gradient due to a tilted illumination direction, but in practice undesired noise is always present. Such noise can be removed by subtracting $\theta _\text {bg}$ and adding the expected phase gradient $\textbf{q}\cdot \textbf{r}$ from the measured phase. The scattered field is processed to reconstruct the refractive index inverting the scalar wave equation in the Born approximation. In the Rytov approximation, the field is not decomposed as a sum of two terms $e^{i\textbf{q}\cdot \textbf{r}}+U_B$, i.e., incident plus scattered field, but it is decomposed as $e^{i\textbf{q}\cdot \textbf{r}}e^{\theta _R}$ that is the product of the incident field times a perturbation term $e^{\theta _R}$. $\theta _R$ can be obtained starting form the measured fields as:
$$\theta_R=\log\left( \frac{\sqrt{I} e^{i\theta}}{ \sqrt{ I_\text{bg}}e^{i\theta_\text{bg}}}\right)=i\left(\theta-\theta_\text{bg}\right)+\frac{1}{2}\log\left(\frac{I}{I_\text{bg}}\right).$$

In our experiments, we used 38 different directions for sample illumination. These tilted illumination beams were produced by displaying sequentially on the DMD sub-matrices of 30$\times$30 pixels in the “on” state, centered on a square grid with a step of 50 pixels. When the beam is focused by the condenser lens, the sample is reached by a quasi plane wave with an incident angle (in air) of up to 52$^\circ$, corresponding to the maximum achievable angle of our condenser lens ($\text {NA}=0.79$). The spatial coherence is approximately $\ell _c\approx \lambda /fM_\text {L2}W=4$ µm, where $f=16$ mm is the focal length of L3 in Fig. 1(a), $M_\text {L2}=9$ is the magnification due to the lens L2 reimaging the DMD on the back-focal plane of L3, and $W=7.45$ µm $\times$ 30 pixels is the size of the pixel sub-matrix in the “on” state.

Starting from Eq. (2), a good measurement of the field intensity and phase is obtained when $I(\textbf{r})\approx I(\textbf{r}+\boldsymbol{\Delta })$ and $\theta (\textbf{r}+\boldsymbol{\Delta })-\theta (\textbf{r})\approx \boldsymbol{\nabla }\theta (\textbf{r})\cdot \boldsymbol{\Delta }$, which is the case if $\Delta$ is much smaller than the inverse of the maximum observable spatial frequency $\lambda /\text {NA}\approx 0.4$ $\mu$m. To observe the interference term of Eq. (2), $\Delta$ must be also smaller than the transverse spatial coherence length $\ell _c\approx 4$ $\mu$m. Considering the microscope magnification, we choose the value $\Delta =0.07$ µm for the SLM displacement, a quantity that is well below $\ell _c$ and $\lambda /\text {NA}$. A lower bound for $\Delta$ is given by noise that is inevitably present in $I_\text {camera}^{\phi }(\textbf{r})$. As $\Delta$ is lowered, the term $\theta (\textbf{r}+\boldsymbol{\Delta })-\theta (\textbf{r})\approx \boldsymbol{\nabla }\theta (\textbf{r})\cdot \boldsymbol{\Delta }$ approaches zero and thus becomes hardly detectable through Eq. (3).

The entire acquisition procedure takes about 2 minutes. This time lapse is mainly occupied by shutter time, which has to be large to compensate the low light intensity. Considering all the illumination directions, approximately 30 seconds are spent to display on the SLM the phase masks that give the four shifts of the global phase $[0,\pi /2,\pi,3\pi /2]$ and the four tilts $\boldsymbol{\Delta }=(\pm \Delta,0)$ and $\boldsymbol{\Delta }=(0,\pm \Delta )$.

2.2 Tomographic reconstruction

Before describing our reconstruction algorithm based on the Rytov approximation, we will first introduce the Fourier diffraction tomography method based on the Born approximation. We start from the scalar wave equation and consider a sample that is illuminated by a plane wave:

$$\nabla^2U_B+k_0^2U_B\approx{-}2k_0^2\delta me^{i\textbf{q}\cdot\textbf{r}}$$
where $k_0=2\pi n_0/\lambda$ with $n_0$ being the medium refractive index, $\textbf{r}$ is the coordinate in the 3D space, $\textbf{q}$ is the illumination wavevector and $\delta m=(n-n_0)/n_0$ is the relative refractive index contrast of the sample. The equation above is obtained by expressing the field as $U=e^{i\textbf{q}\cdot \textbf{r}}+U_B$ and neglecting the term $\delta m U_B$ meaning that we are working under the assumption of single scattering. Also, the difference of the squared of the refractive indices $(n^2-n_0^2)$ is approximated as $2n_0^2\delta m$. By taking the Fourier transform of Eq. (6), one can isolate $U_B$:
$$U_B(\textbf{k}) ={-}\frac{2 k_0^2}{k_0^2-k^2} \delta m(\textbf{k}-\textbf{q})$$
where $\textbf{k}$ is the 3D coordinate in the reciprocal space. Equation (7) relates $\delta m$ with the 3D scattered field that also contains unmeasurable near field contributions. To express the 3D map of the relative refractive index contrast $\delta m$ as a function of the 2D measured field one must compute the inverse Fourier transform at $z=0$:
$$U_B(k_x,k_y,z=0) = \frac{2\pi i k_0^2}{K_z} \delta m(k_x-q_x,k_y-q_y,K_z-q_z)$$
where $K_z=\sqrt {k_0^2-k_x^2-k_y^2}$ and $U_s(k_x,k_y,z=0)$ is the 2D Fourier transform of the observed scattered field. Equation (8) is at the base of Fourier diffraction theorem and can be used to fill the 3D Fourier space of $\delta m$ upon changing the incident wavevector $\textbf{q}$ [1320]. To avoid reconstruction artifacts in $\delta m$, complex interpolation algorithms [33] must be used to fit $U_B$ in the 3D reciprocal space. Alternatively, Fourier interpolation can be avoided by constructing $\delta m$ directly in the real space. Given the illumination wavevector $\textbf{q}$, the measured scattered field $U_B^{(q)}(x,y,z=0)$ can be back-propagated in the 3D space and then multiplied by $\exp (i\textbf{q}\cdot \textbf{r})$ so that, in the reciprocal space, the only Fourier components populated are the ones laying on a translated hemisphere as in Eq. (8). By adding up these terms for every $\textbf{q}$, the entire observable Fourier components are populated. However, since the hemispheres in the Fourier space partially overlap, a normalization is necessary:
$$\delta m(\textbf{r})= N(\textbf{r}) \otimes \left( \frac{-1}{2\pi i k_0}\sum_\textbf{q} U_B^{(q)}(\textbf{r})e^{{-}i\textbf{q}\cdot\textbf{r}}\right )$$
where the term $N$ is the normalization filter. In Fourier transform based tomography, real-space back-projection is sometimes used to obtain an image phantom that must be filtered to compute the correct image [34]. Here, the normalization filter $N$ is built in the reciprocal space as follows: we start from a 3D array filled with zeros and for each $\textbf{q}$ the value 1 is added to the voxels through which the hemisphere $k_z-q_z=(k_0^2 - (k_x-q_x)^2- (k_y-q_y)^2)^{1/2}$ passes. Finally, for the voxels with value different from zero, the inverse is taken. A summary of the steps described above is shown in Fig. 2.

 figure: Fig. 2.

Fig. 2. Scheme of the reconstruction algorithm of Eqs. (9),(11). For each illumination wavevector $\textbf{q}$, the input field is back-propagated in the 3D space. In the Born approximation, the input field is the scattered field while, in the Rytov approximation, it is constructed as $[\theta _Re^{i\textbf{q}\cdot \textbf{r}}]_{z=0}$ where $\theta _R(z=0)$ is computed as in Eq. (5). The Fourier transform of the back-propagated field stack has non-zero components on a hemisphere with radius $k_0$. Each stack is multiplied by $e^{-i\textbf{q}\cdot \textbf{r}}$ that is equivalent to a shift in the Fourier space. The Fourier space is filled by summing over all $\textbf{q}$. A convolution with a normalization function ensures that all the Fourier components are weighted properly.

Download Full Size | PDF

The Born approximation is accurate only for low contrast samples ($\delta m\ll k_0^{-1}L^{-1}$ for an object with size $L$). A more robust approximation is the one given by Rytov [5], in which the field is decomposed as $U=e^{i\textbf{q}\cdot \textbf{r}}e^{\theta _R}$. The phase $\theta _R$ can be expressed as a series whose first order term can be found by solving a differential equation whose form is close to the scalar wave equation (Eq. (20) in [4]):

$$\left(\nabla^2+k_0^2\right)\left[e^{i\textbf{q}\cdot\textbf{r}}\theta_R\right]\approx2k_0^2\delta m.$$

Equation (10) is similar to Eq. (6) with $\delta m e^{i\textbf{q}\cdot \textbf{r}}$ taking the place of $\delta m$ and $U_B$ being replaced by the fictitious field $\hat {U}=\theta _Re^{i\textbf{q}\cdot \textbf{r}}$. Following the same steps that lead to Eq. (9), we can conclude that:

$$\delta m(\textbf{r})= N(\textbf{r}) \otimes \left( \frac{1}{2\pi i k_0}\sum_\textbf{q} \hat{U}(\textbf{r})e^{{-}i\textbf{q}\cdot\textbf{r}}\right)$$
where $\hat {U}(\textbf{r})$ is computed by back-propagating in the 3D space the 2D field $[\theta _Re^{i\textbf{q}\cdot \textbf{r}}]_{z=0}$ which is obtained experimentally as in Eq. (5).

3. Results

To validate our technique, and to identify the range of experimental conditions in which our method is capable of producing quantitatively reliable refractive index tomograms, we imaged spherical particles of different sizes and at different refractive index contrasts with respect to the medium in which they are immersed. Slices of the 3D reconstructions along the $x-y$ and $x-z$ planes passing through the center of each sphere are shown together with plots of the refractive index contrast along $x$ in Fig. 3. In all reconstructions, a priori information on the samples (i.e., a strictly positive or a strictly negative refractive index contrast) is used by zeroing all voxels that violate these constraints. These artefacts in the reconstructions are due to incomplete sampling of the Fourier space, as shown in Fig. 2 (see also Appendix A for a more extensive discussion).

 figure: Fig. 3.

Fig. 3. 3D reconstructions of the refractive index contrast with respect to the medium ($\delta n=n-n_0$) for four different microspheres. Slices of the $x-y$ and the $x-z$ planes are shown together with a plot along the axis parallel to the $x$ direction passing through the microsphere’s center. The reconstruction of the 2 µm silica bead shown in (a) is based on the Born approximation, while the other tomograms (b)-(e) are reconstructed using Rytov’s approximation. For the microspheres with a refractive index higher than the medium (a),(b),(d),(e) the condition of positive contrast is imposed by setting to zero the voxels for which $\delta n<0$; the opposite is done for the silica bead in glycerol (c) where the unphysical values $\delta n>0$ are set to zero. Dashed lines indicate the expected size of the spherical bead. All scale bars are 2 µm

Download Full Size | PDF

A comparison between the reconstructed refractive index maps of a 2 µm silica microsphere immersed in water (refractive index 1.334 at 530 nm) using Born and Rytov approximations is shown in Figs. 3(a),(b). Silica beads’ refractive index is often smaller than the value for fused silica (1.457 at 530 nm) as a consequence of the fabrication process that leaves pores in the structure and can vary between 1.37 and 1.46 [35,36]. The reconstruction based on the Born approximation yields a tomogram with non-uniform values of the refractive index within the bead’s volume, and the particle’s shape appears irregular and elongated in the $z$ direction. Conversely, the Rytov algorithm provides a more regular reconstruction with a refractive index of $1.43\pm 0.01$; small fluctuations around this value are present in the volume occupied by the microsphere. The tomograms produced by the two approximations differ significantly as a consequence of the relatively high refractive index contrast ($\delta n\simeq 0.1$) of the silica microbead in water, confirming that the Rytov approach provides a much more faithful reconstruction when the single-scattering hypothesis starts to fail. We note that the measured value of the refractive index is consistent with the expected range of values for (porous) silica microspheres.

Figures 3(c)-(e) also shows more tests of the Rytov algorithm with samples of varying microspheres’ size and refractive index contrast. In particular, we analysed the case of negative refractive index contrast (panel c) using 3.7 µm silica beads in glycerol (refractive index 1.474 at 530 nm) and the case of a larger sphere size of 7 µm (panel d) silica in water. Notably, we obtained tomograms with refractive index reasonably constant within particles’ volume, with values of $1.394\pm 0.006$ and $1.48\pm 0.03$, respectively. The first value lays within the expected range while the second is slightly larger then expected. Differently, the reconstruction shown in Fig. 3(e) of the 5 µm polystyrene sphere (refractive index 1.599 at 530 nm) in water shows a non-uniform distribution of the refractive index together with contrast values that are lower than the expected $\delta n=0.265$. This indicates that a moderately high refractive index contrast is challenging also for Rytov’s method. We note that tomographic reconstructions of colloidal beads are difficult to be found in the literature as usually these techniques are tested against biological samples whose shape is not known in advance. To the best of our knowledge, the only exception is a work in which polystyrene 3 µm diameter beads immersed in microscope oil ($\delta n \approx ~0.08$) are reconstructed [37]. As label-free tomography has shown to be a useful tool for 3D imaging of biological samples, here our technique is tested on a microglia cell adhering to the substrate of an Ibidi µ-Slide multi-well plate. The cells were kept at 37$^\circ$ Celsius via a custom built thermostat mounted on the microscope stage. A 3D rendering of the cell refractive index and the $x-y$ slice of the plane on the substrate-medium interface are shown in Fig. 4(d). As compared to spheres, here the contrast is very low ($\delta n$ in the range $0.025-0.05$, with $n_\text {medium}=1.334$) and the cell shows a thick body with thin and wide (north in the image) or long (south) protrusions laying on the substrate. Even in this case the Rytov algorithm yields a good reconstruction with refractive index values between $\approx$ 1.36 and $\approx$ 1.39, as expected for most cells [11,38].

 figure: Fig. 4.

Fig. 4. 3D reconstructions of SU8 microstructures (refractive index 1.603) fabricated with two-photon polymerization: pyramid in oil (a), arch in oil (b) and dumbbell in water (c). For each structure 2D sections of the reconstruction are show more quantitatively the structure size. Each tomographic reconstruction is paired with a corresponding SEM image for comparison. (d) 3D reconstruction of the refractive index map of a microglia cell. Inset shows a 2D slice of the reconstruction at the coverglass interface. All color bars indicate the refractive index contrast $\delta n$. The cell refractive index contrast with respect to the medium lays between 0.025 and 0.05 as expected for most cells [38].

Download Full Size | PDF

We showed that our optical tomography based on the Rytov approximation is capable to extract reliable tomograms of samples that, in the case of micro-beads, have a high contrast and a simple spherical shape, and, in the case of a cell, have a low contrast and a complex shape. Objects with known, non-trivial shape and contrast $\delta n\approx ~0.1$ represent a challenging goal, open to the possibility of checking whether the reconstructed shape reproduces accurately the real one. To this aim, we imaged three SU-8 microstructures (refractive index 1.603 at 530 nm), fabricated via 2-photon polymerization [39,40] as test objects. Figure 4 shows the three reconstructions, a dumbbell in water (a), a pyramid in oil (b), and a square arch in oil (c), along with the corresponding electron microscopy images (d),(e),(f). Scanning Electron Microscopy (SEM) images of microstructures sputtered with a chrome layer were acquired using a FEI Quanta 400 microscope. Imaging was performed on samples tilted by $22^\circ$ with respect to the vertical direction, using an extracting voltage of 20 kV. In the case of the dumbbell (Fig. 4(c)), despite the refractive index contrast is very high ($\delta n=0.269$), it is possible to reconstruct the two spheres as well as to resolve the thin bar that connects them. We threshold the 3D map to remove background noise and compare the structure size, which we assume to occupy the voxels above threshold, with the one measured on the electron microscopy images. For the sphere we measured an inplane diameter of 2.2 µm and a size approximately 2.4 along the optical axis, which compares well with the value of 2.5 µm measured with the SEM, while for the bar thickness we measure 0.66 µm that is also close to the value of 0.63 µm obtained with the SEM. A fair agreement between the tomographic and the SEM images is achieved but the high refractive index contrast respect to the surrounding medium makes the reconstruction quite challenging. Reconstruction of SU8 structures in oil (refractive index 1.518 at 530 nm, $\delta n=0.085$), having a much lower contrast, are shown in Figs. 4(a),(b).

Figure 4(a) shows a 3D reconstruction of a pyramid. The vertical slice of the reconstruction, passing through the apex and the midpoint of the base edges, shows that the angle (60 degrees by design) is compatible with our observations. Similarly horizontal slices show clearly the square base of the pyramid. The arch is a more challenging structure (Fig. 4(b)), as a consequence of the horizontal elevated crossbar and the void space below it. We measured a leg thickness of 2.9 µm (2.1 µm on SEM), the total length of the crossbar 7.0 µm (6.8 on SEM), the crossbar section has a vertical and horizontal size of respectively 4.1 and 1.7 µm (4.0 and 1.7 µm on SEM).

4. Conclusions

We developed a new tomographic microscopy technique that uses an LED source and a common-path interferomenter. Unlike most of tomographic techniques using low temporal coherence light [2,3,7,9,15], we used an algorithm based on Rytov approximation that provides a fast sample reconstruction while being not limited in terms of sample contrast. The algorithm requires the knowledge of the complex transmitted field that is not easily accessible if a low coherence source is used. We addressed the problem of measuring the phase of a field with low temporal coherence by acquiring multiple images, each corresponding to a controlled phase modulation applied by an SLM placed in the image Fourier plane. The advantage of our method lies in the use of a common-path interferometer that, unlike heterodyne interferometers, does not suffer from phase fluctuations and does not require any optical path compensation in the reference arm, which in the case of low temporal coherence light must be extremely precise. We characterized our technique and tested it on both biological and synthetic samples, and for the second ones we compared the obtained 3D reconstructions with SEM images. Results assessed the robustness of the Rytov approximation and demonstrated that we can get a reliable tomogram of small structures like a dumbell (a submicron thick bar connecting two $\approx$ 2 µm diameter spheres) with high refractive index contrast ($\approx$ 0.27) with respect to the water medium. In the case of the pyramids and the arch, the bigger size of approximately 10 µm forced us to reduce the contrast to $\approx$ 0.08, by using microscope oil as the medium around the structures, to get faithful tomographic reconstructions. Our method provides reliable reconstructions of 3D refractive index maps of both biological and microfabricated structures. These features make it particularly useful for 3D imaging of complex structures that integrate biological components within artificial microstructures [41].

Appendix A: Fourier sampling and sign constraints

Here we want to discuss the effect of the limited sampling of the refractive index map in the Fourier space on the refractive index reconstruction in the real space. Each illumination direction allows us to sample the refractive index in the Fourier space on a hemisphere as shown in Fig. 2. By increasing the number of incident angles, the reciprocal space is more densely sampled but the regions above and below the red area represented in Fig. 2 will always remain uncovered. To cover the missing space one should be able to rotate the specimen while keeping it relatively stable in $x,y,z$ [18,21], but this is impracticable in most cases.

In Fig. 5(a) we plot as solid blue line the refractive index profile along $x$ of the same 2 $\mu$m silica microsphere shown in Fig. 3(b) along with the $x-y$ and $x-z$ slices. Note that there is a region around the microsphere where the refractive index is lower than the medium one. In Fig. 3(b) we used the a priori knowledge on the sample and set to zero the voxels having negative refractive index contrast. An alternative way to impose this positivity condition is described in [23]: i) set to zero the negative refractive index voxels; ii) take the inverse Fourier transform (now the sampled regions in the Fourier space have different values and the unsampled regions are no longer zero); iii) replace the accessible components of the Fourier space with the measured data; iv) take the inverse Fourier transform. This procedure is iterated until the algorithm converges. The obtained profile is plotted by the orange line in Fig. 5(a). Most of the voxels have a positive contrast but the uniformity inside the region occupied by the microsphere is lost as can be seen from the $x-y$ and the $x-z$ slices in Fig. 5(c). For this reason we impose sign constraints using the simpler procedure described in the text.

 figure: Fig. 5.

Fig. 5. (a) Reconstructions of the 2 $\mu$m silica bead in water shown in Fig. 3(b). The blue solid line plots the refractive index contrast before the voxels with negative contrast are set to 0. The $x-y$ and the $x-z$ slices of the same reconstruction are shown in (b). The orange line plots an alternative way to impose the a priori knowledge on the sample as explained in the appendix. The $x-y$ and the $x-z$ slices for the same reconstruction are shown in (c). (d),(e) Simulated reconstructions of a silica bead in water for different numbers of incident angles. (d) Profiles along the $x$ axis. (e) Profiles along the optical axis $z$.

Download Full Size | PDF

Finally, we investigated how changing the number of the incident beams affects the reconstructed image. For this purpose, using Mie theory [42], we computed the field scattered by a 2 $\mu$m silica bead immersed in water, assuming $\delta n=0.1$. The incident beam is modeled as a monochromatic plane wave with $\lambda =0.53$. The incident wavevectors are distributed uniformly on spherical cap and the maximum incident angle is the same one of our condenser, corresponding to a NA of 0.8. Figure 5(d) plots the reconstructed refractive index contrast along the $x$ axis while Fig. 5(e) plots the same quantity along the optical axis $z$. As the number of incident angles increases the refractive index value profile inside the particle becomes flatter, the image sharpens, and the area with negative refractive index appears. In particular, the reconstruction sharpens along the optical axis $z$. Conversely, with a low number of incident angles, the refractive index contrast does not drop to zero, even at a distance of $4$ $\mu$m from the bead center.

Funding

H2020 European Research Council (834615); Regione Lazio (B86C18000970002).

Acknowledgments

We thank Maria Rosito (Center for Life Nano- & Neuro-Science, Fondazione Istituto Italiano di Tecnologia (IIT), Rome, Italy) for providing the microglia cells examined in this work and Marco Albano (Scanning Electron Microscope Laboratory of the Earth Sciences Department, Sapienza University of Rome) for acquiring the SEM images shown in this work.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. T. Kim, R. Zhou, M. Mir, S. D. Babacan, P. S. Carney, L. L. Goddard, and G. Popescu, “White-light diffraction tomography of unlabelled live cells,” Nat. Photonics 8(3), 256–263 (2014). [CrossRef]  

2. J. M. Soto, J. A. Rodrigo, and T. Alieva, “Label-free quantitative 3D tomographic imaging for partially coherent light microscopy,” Opt. Express 25(14), 15699–15712 (2017). [CrossRef]  

3. J. A. Rodrigo, J. M. Soto, and T. Alieva, “Fast label-free microscopy technique for 3D dynamic quantitative imaging of living cells,” Biomed. Opt. Express 8(12), 5507–5517 (2017). [CrossRef]  

4. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Elsevier, 2013).

5. B. Chen and J. J. Stamnes, “Validity of diffraction tomography based on the first Born and the first Rytov approximations,” Appl. Opt. 37(14), 2996–3006 (1998). [CrossRef]  

6. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Learning approach to optical tomography,” Optica 2(6), 517–522 (2015). [CrossRef]  

7. L. Tian and L. Waller, “3D intensity and phase imaging from light field measurements in an LED array microscope,” Optica 2(2), 104–111 (2015). [CrossRef]  

8. T. Zhang, C. Godavarthi, P. C. Chaumet, G. Maire, H. Giovannini, A. Talneau, M. Allain, K. Belkebir, and A. Sentenac, “Far-field diffraction microscopy at λ/10 resolution,” Optica 3(6), 609–612 (2016). [CrossRef]  

9. S. Chowdhury, M. Chen, R. Eckert, D. Ren, F. Wu, N. Repina, and L. Waller, “High-resolution 3D refractive index microscopy of multiple-scattering samples from intensity images,” Optica 6(9), 1211–1219 (2019). [CrossRef]  

10. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comput. Imaging 2(1), 59–70 (2016). [CrossRef]  

11. W. Choi, C. Fang-Yen, K. Badizadegan, S. Oh, N. Lue, R. R. Dasari, and M. S. Feld, “Tomographic phase microscopy,” Nat. Methods 4(9), 717–719 (2007). [CrossRef]  

12. F. Merola, P. Memmolo, L. Miccio, R. Savoia, M. Mugnano, A. Fontana, G. D’ippolito, A. Sardo, A. Iolascon, A. Gambale, and P. Ferraro, “Tomographic flow cytometry by digital holography,” Light: Sci. Appl. 6(4), e16241 (2017). [CrossRef]  

13. 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]  

14. Y. Kim, H. Shim, K. Kim, H. Park, J. H. Heo, J. Yoon, C. Choi, S. Jang, and Y. Park, “Common-path diffraction optical tomography for investigation of three-dimensional structures and dynamics of biological cells,” Opt. Express 22(9), 10398–10407 (2014). [CrossRef]  

15. R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang, “Diffraction tomography with Fourier ptychography,” Optica 3(8), 827–835 (2016). [CrossRef]  

16. J. Jung, K. Kim, J. Yoon, and Y. Park, “Hyperspectral optical diffraction tomography,” Opt. Express 24(3), 2006–2012 (2016). [CrossRef]  

17. K. Lee, K. Kim, G. Kim, S. Shin, and Y. Park, “Time-multiplexed structured illumination using a DMD for optical diffraction tomography,” Opt. Lett. 42(5), 999–1002 (2017). [CrossRef]  

18. B. Simon, M. Debailleul, M. Houkal, C. Ecoffet, J. Bailleul, J. Lambert, A. Spangenberg, H. Liu, O. Soppera, and O. Haeberlé, “Tomographic diffractive microscopy with isotropic resolution,” Optica 4(4), 460–463 (2017). [CrossRef]  

19. S. Shin, D. Kim, K. Kim, and Y. Park, “Super-resolution three-dimensional fluorescence and optical diffraction tomography of live cells using structured illumination generated by a digital micromirror device,” Sci. Rep. 8(1), 9183 (2018). [CrossRef]  

20. D. Kim, N. Oh, K. Kim, S. Lee, C.-G. Pack, J.-H. Park, and Y. Park, “Label-free high-resolution 3-D imaging of gold nanoparticles inside live cells using optical diffraction tomography,” Methods 136, 160–167 (2018). [CrossRef]  

21. M. Ziemczonok, A. Kuś, P. Wasylczyk, and M. Kujawińska, “3D-printed biological cell phantom for testing 3D quantitative phase imaging systems,” Sci. Rep. 9(1), 18872 (2019). [CrossRef]  

22. K. Lee, S. Shin, Z. Yaqoob, P. T. So, and Y. Park, “Low-coherent optical diffraction tomography by angle-scanning illumination,” J. Biophotonics 12(5), e201800289 (2019). [CrossRef]  

23. 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]  

24. P. Ferraro, A. Wax, and Z. Zalevsky, Coherent Light Microscopy: Imaging and Quantitative Phase Analysis, vol. 46 (Springer Science & Business Media, 2011).

25. C. Maurer, A. Jesacher, S. Bernet, and M. Ritsch-Marte, “What spatial light modulators can do for optical microscopy,” Laser Photonics Rev. 5(1), 81–101 (2011). [CrossRef]  

26. G. Popescu, L. P. Deflores, J. C. Vaughan, K. Badizadegan, H. Iwai, R. R. Dasari, and M. S. Feld, “Fourier phase microscopy for investigation of biological structures and dynamics,” Opt. Lett. 29(21), 2503–2505 (2004). [CrossRef]  

27. C. Wang, N. Chen, and W. Heidrich, “Towards self-calibrated lens metrology by differentiable refractive deflectometry,” Opt. Express 29(19), 30284–30295 (2021). [CrossRef]  

28. N. Chen, C. Wang, and W. Heidrich, “Compact computational holographic display (invited article),” Front. Photonics 3, 1 (2022). [CrossRef]  

29. M. I. Shribak, J. LaFountain, D. S. Biggs, and S. Inoue, “Orientation-independent differential interference contrast microscopy and its combination with an orientation-independent polarization system,” J. Biomed. Opt. 13(1), 014011 (2008). [CrossRef]  

30. M. Shribak, K. G. Larkin, and D. Biggs, “Mapping optical path length and image enhancement using quantitative orientation-independent differential interference contrast microscopy,” J. Biomed. Opt. 22(1), 016006 (2017). [CrossRef]  

31. T. H. Nguyen, M. E. Kandel, M. Rubessa, M. B. Wheeler, and G. Popescu, “Gradient light interference microscopy for 3D imaging of unlabeled specimens,” Nat. Commun. 8(1), 210 (2017). [CrossRef]  

32. Z. Wang, L. Millet, M. Mir, H. Ding, S. Unarunotai, J. Rogers, M. U. Gillette, and G. Popescu, “Spatial light interference microscopy (SLIM),” Opt. Express 19(2), 1016–1026 (2011). [CrossRef]  

33. S. Lanzavecchia and P. Luigi Bellon, “Fast computation of 3D radon transform via a direct Fourier method,” Bioinforma. 14(2), 212–216 (1998). [CrossRef]  

34. M. Figl, “CT reconstruction,” in Applied Medical Image Processing: A Basic Course, p. 319 (CRC Press, 2011).

35. F. C. Cheong, K. Xiao, D. J. Pine, and D. G. Grier, “Holographic characterization of individual colloidal spheres’ porosities,” Soft Matter 7(15), 6816–6819 (2011). [CrossRef]  

36. F. Garcia-Santamaria, H. Miguez, M. Ibisate, F. Meseguer, and C. Lopez, “Refractive index properties of calcined silica submicrometer spheres,” Langmuir 18(5), 1942–1944 (2002). [CrossRef]  

37. 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]  

38. P. Y. Liu, L. K. Chin, W. Ser, H. F. Chen, C. M. Hsieh, C. H. Lee, K. B. Sung, T. C. Ayi, P. H. Yap, B. Liedberg, K. Wang, T. Bourouinaj, and Y. Leprince-Wang, “Cell refractive index for cell biology and disease diagnosis: past, present and future,” Lab Chip 16(4), 634–644 (2016). [CrossRef]  

39. S. Maruo, O. Nakamura, and S. Kawata, “Three-dimensional microfabrication with two-photon-absorbed photopolymerization,” Opt. Lett. 22(2), 132–134 (1997). [CrossRef]  

40. S. Bianchi, G. Vizsnyiczai, S. Ferretti, C. Maggi, and R. Di Leonardo, “An optical reaction micro-turbine,” Nat. Commun. 9(1), 4476 (2018). [CrossRef]  

41. G. Vizsnyiczai, A. Búzás, B. L. Aekbote, T. Fekete, I. Grexa, P. Ormos, and L. Kelemen, “Multiview microscopy of single cells through microstructure-based indirect optical manipulation,” Biomed. Opt. Express 11(2), 945–962 (2020). [CrossRef]  

42. F. Saglimbeni, S. Bianchi, G. Bolognesi, G. Paradossi, and R. Di Leonardo, “Optical characterization of an individual polymer-shelled microbubble structure via digital holography,” Soft Matter 8(34), 8822–8825 (2012). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (5)

Fig. 1.
Fig. 1. (a) Experimental setup. The sample is illuminated by a low temporal coherence source. The illumination beam’s angle and aperture are selected using a DMD whose image is reformed by the lenses L1 and L2 on the back-focal plane of the condenser L3. Light is polarized at 45$^\circ$ with respect to the SLM active axis by the first polarizer P. On the camera plane, the polarization component along the SLM active axis is laterally shifted by $\Delta$ as a consequence of a linear grating displayed on the SLM. A second polarization analyzer projects the two polarization components on a 45$^\circ$ axis so that the two interfere. (b) Given a shift direction ${\bf \Delta }$, the SLM changes the relative phase of the two fields so that the gradient of the phase can be computed as in Eq. (3). The phase gradients computed with ${\bf \Delta }=(\pm \Delta,0)$ are combined to obtain a cleaner gradient along $x$; the same is done for ${\bf \Delta }=(0,\pm \Delta )$. The two phase gradients $\partial _x\theta$ and $\partial _y\theta$ are differentiated numerically to compute the Laplacian that is integrated to finally obtain $\theta$. Scalebar is 2 µm.
Fig. 2.
Fig. 2. Scheme of the reconstruction algorithm of Eqs. (9),(11). For each illumination wavevector $\textbf{q}$, the input field is back-propagated in the 3D space. In the Born approximation, the input field is the scattered field while, in the Rytov approximation, it is constructed as $[\theta _Re^{i\textbf{q}\cdot \textbf{r}}]_{z=0}$ where $\theta _R(z=0)$ is computed as in Eq. (5). The Fourier transform of the back-propagated field stack has non-zero components on a hemisphere with radius $k_0$. Each stack is multiplied by $e^{-i\textbf{q}\cdot \textbf{r}}$ that is equivalent to a shift in the Fourier space. The Fourier space is filled by summing over all $\textbf{q}$. A convolution with a normalization function ensures that all the Fourier components are weighted properly.
Fig. 3.
Fig. 3. 3D reconstructions of the refractive index contrast with respect to the medium ($\delta n=n-n_0$) for four different microspheres. Slices of the $x-y$ and the $x-z$ planes are shown together with a plot along the axis parallel to the $x$ direction passing through the microsphere’s center. The reconstruction of the 2 µm silica bead shown in (a) is based on the Born approximation, while the other tomograms (b)-(e) are reconstructed using Rytov’s approximation. For the microspheres with a refractive index higher than the medium (a),(b),(d),(e) the condition of positive contrast is imposed by setting to zero the voxels for which $\delta n<0$; the opposite is done for the silica bead in glycerol (c) where the unphysical values $\delta n>0$ are set to zero. Dashed lines indicate the expected size of the spherical bead. All scale bars are 2 µm
Fig. 4.
Fig. 4. 3D reconstructions of SU8 microstructures (refractive index 1.603) fabricated with two-photon polymerization: pyramid in oil (a), arch in oil (b) and dumbbell in water (c). For each structure 2D sections of the reconstruction are show more quantitatively the structure size. Each tomographic reconstruction is paired with a corresponding SEM image for comparison. (d) 3D reconstruction of the refractive index map of a microglia cell. Inset shows a 2D slice of the reconstruction at the coverglass interface. All color bars indicate the refractive index contrast $\delta n$. The cell refractive index contrast with respect to the medium lays between 0.025 and 0.05 as expected for most cells [38].
Fig. 5.
Fig. 5. (a) Reconstructions of the 2 $\mu$m silica bead in water shown in Fig. 3(b). The blue solid line plots the refractive index contrast before the voxels with negative contrast are set to 0. The $x-y$ and the $x-z$ slices of the same reconstruction are shown in (b). The orange line plots an alternative way to impose the a priori knowledge on the sample as explained in the appendix. The $x-y$ and the $x-z$ slices for the same reconstruction are shown in (c). (d),(e) Simulated reconstructions of a silica bead in water for different numbers of incident angles. (d) Profiles along the $x$ axis. (e) Profiles along the optical axis $z$.

Equations (11)

Equations on this page are rendered with MathJax. Learn more.

Φ ( r ) = 2 π λ f Δ r + ϕ
I camera ϕ ( r ) = I ( r ) + I ( r + Δ ) + 2 I ( r ) I ( r + Δ ) cos ( θ ( r + Δ ) θ ( r ) + ϕ )
θ Δ = arctan ( I camera π / 2 I camera 3 π / 2 I camera 0 I camera π ) .
U B = ( I ( r ) e i θ ( r ) I bg ( r ) e i θ bg ( r ) ) e i θ bg ( r ) + i q r
θ R = log ( I e i θ I bg e i θ bg ) = i ( θ θ bg ) + 1 2 log ( I I bg ) .
2 U B + k 0 2 U B 2 k 0 2 δ m e i q r
U B ( k ) = 2 k 0 2 k 0 2 k 2 δ m ( k q )
U B ( k x , k y , z = 0 ) = 2 π i k 0 2 K z δ m ( k x q x , k y q y , K z q z )
δ m ( r ) = N ( r ) ( 1 2 π i k 0 q U B ( q ) ( r ) e i q r )
( 2 + k 0 2 ) [ e i q r θ R ] 2 k 0 2 δ m .
δ m ( r ) = N ( r ) ( 1 2 π i k 0 q U ^ ( r ) e i q r )
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.