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

Numerical refractive index correction for the stitching procedure in tomographic quantitative phase imaging

Open Access Open Access

Abstract

Tomographic quantitative phase imaging (QPI) lacks an absolute refractive index value baseline, which poses a problem when large dense objects extending over multiple fields of view are measured volume by volume and stitched together. Some of the measurements lack the natural baseline value that is provided by the mounting medium with a known refractive index. In this work, we discuss the problem of the refractive index (RI) baseline of individual reconstructed volumes that are deprived of access to mounting medium due to the extent of the object. The solution of this problem is provided by establishing the RI offsets based on the overlapping regions. We have proven that the process of finding the offset RI values may be justifiably reduced to the analogous procedure in the 2D baseline correction (2D-BC). Finally, we proposed the enhancement of the state-of-the-art 2D-BC procedure previously introduced in the context of 2D QPI. The processing is validated at the examples of a synthetic dataset and a liver organoid.

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

1. Introduction

Tomographic quantitative phase imaging (QPI) is a 3D label-free imaging technique commonly used for analysis of biological specimens [15]. The technique retrieves spatial distribution of the refractive index (RI). When an analyzed sample is larger than the field of view (FoV) of the imaging device, it is necessary to capture multiple measurements of fragments of the object, reconstruct them independently and finally stitch them together to create a single large tomographic reconstruction [68]. Before the stitching procedure can be applied it is required to set a common coordinate system which is performed based on the spatial overlap of the adjacent measurements. The overlapping regions allow for precise positioning of the measurements by mutually aligning features captured in neighboring volumes. There are already developed procedures, such as the open-source BigStitcher [9], for determining the coordinate system of an arbitrary set of 3D images. However such procedures are not designed to match the baseline RI values in the stitched volume. Such a problem is unique to QPI, which fundamentally provides relative measurements.

Although the problem of baseline RI value originates from the basics of wave propagation and is present in the theory of optical diffraction tomography [10] it has not been accentuated in the literature in the context of 3D QPI stitching up to this point, most likely due to each reconstructed volume containing the mounting medium regions that pose a natural baseline RI value. In such cases where the mounting medium is present in each FoV the relative RI values are changed to absolute ones by setting the correct baseline level – an offset value is added such that the background regions are valued at an a priori known RI of the mounting medium. This is a standard treatment of the reconstructed volumes and other authors rightfully did not emphasise it. The problem arises when large samples are measured fragment-by-fragment and the mounting medium is not present in all FoVs. Hence, it is impossible to independently establish the correct offset value, which brings significant errors to the stitched volume. This is the case for objects like dense tissue slices, organoids or embryos.

A similar problem has been described in 2D QPI and addressed with appropriate numerical processing methods [11,12]. Those techniques utilized the overlapping parts of the measurements that are anyway required for fine spatial positioning to unify the baseline values in certain phase images that lack sample-free regions within their FoVs.

According to our knowledge the baseline correction (BC) is a novel consideration in the context of 3D QPI stitching. We begin this work with a theoretical discussion on why 3D-BC is a necessary procedure in the absence of background region within some tomographic measurements within the captured dataset. In the course of the discussion it is shown that the 3D-BC is essentially a problem of finding appropriate additive RI values (offsets) for each volume that lacks the background region. By drawing inspiration from the 2D-BC we propose to establish the offsets based on the overlapping regions. We conclude the theoretical analysis of the 3D-BC with the novel observation, that the process of finding the offset RI values may be justifiably reduced to the analogous procedure in the 2D-BC, hence blurring the line between 2D-BC and 3D-BC.

Based on this observation we introduce a novel BC method and compare it on the synthetic datasets with the previous SotA BC method that was initially introduced exclusively in the context of 2D-BC. We conclude the work by presenting baseline-corrected stitching of an experimental dataset – a measurement of a liver organoid.

2. Baseline in 3D QPI measurement: discussion and problem formulation

Typically in tomographic QPI a sequence of off-axis projections of an object are captured in order to reconstruct the objects’ 3D RI distribution. Any QPI modality may be utilized to perform such volumetric reconstruction, as exemplified by Fourier ptychography [13], lensless holographic imaging [14] and holographic tomography [15].

Most tomographic approaches comprise of multiple 2D QPI measurements (taken at different illumination angles or sample angular positions) coupled with a suitable tomographic reconstruction algorithm [10,1621]. Other 3D QPI approaches, that process intensity measurements [6,2225], indirectly leverage the connection between phase delay and 3D RI distribution by computation of multiple propagation steps. Hence, they may only model the phase differences that are present within a single FoV, regardless of the presence of the mounting medium.

With the a priori knowledge of the RI of the mounting medium the absolute RI of the object may be calculated, as long as the background region (characterized by the RI of the mounting medium) is present within the FoV.

Without the background region the phase differences captured for all illumination angles would have to be baseline-corrected as shown in the works regarding 2D QPI BC [11,12]. In general it can be stated that the phase shift is coupled with object’s RI by the following equation:

$$\begin{aligned}\mathrm{\Phi}(x, y) &= \int_L \Delta n(x, y, z) \mathrm{d}l,\\ \mathrm{therefore,}&\\ \mathrm{\Phi}(x, y) + C &= \int_L \left( \Delta n(x, y, z) + C' \right) \mathrm{d}l, \end{aligned}$$
with: $\Phi$ – phase map, $x, y, z$ - spatial coordinates, $\Delta n$ – variation of RI within the FoV, $\mathrm {d}l$ – distance along the beam path $L$, and $C$, $C'$ - constants. Such equation holds not only for measured phase projection used by conventional reconstruction methods, but also for phase shifts modeled by propagation in the intensity-based reconstruction algorithms.

With this equation it is clear that a change in phase projections by a constant offset value results in another offset value in the RI reconstruction. This also means that in order to perform this correction it is not required to correct all captured projections, because the problem may be simplified to the BC of the final reconstructed volume, and hence the BC is adequate for any 3D RI reconstruction method.

3. Methods

In this section we provide an overview of our data acquisition scenario with the overlapping regions (Sec. 3.1) and introduce a novel iterative averaged BC method that utilizes these overlaps (Sec. 3.2). The method is compared to the SotA BC (see Sec. 3.3) optimization-based method (opt-BC) that has originally been introduced as a 2D-BC [12], but may be applied universally as both 2D and 3D-BC method considering the argumentation in Sec. 3.2. Let us again emphasise that only the BC based on the overlapping measurements are within the scope of this work. In cases where the background is present in all measurements it is preferred to perform BC individually in each FoV by using the background regions directly.

3.1 Tomographic measurement setup

In this work we assume a measurement scenario, where a large sample is measured with an arbitrary optical tomographic setup. The specimen is larger than the FoV of the device, thus a motorized XY stage is used to measure it fragment by fragment. Each set of projections is then reconstructed independently and as a result one tomographic reconstruction is generated for each position of the stage. Reconstructions of neighbouring fragments have a common part, called an overlap margin. For the schematic view of described configuration see Fig. 1. We assume, that the overlap percentage is the same for all reconstructions, as it greatly simplifies the BC procedure. It is possible to generalize the implementation of the methods presented here to take into account irregular overlap, however, as we show in Sec. 4, the current implementation is resistant to minor overlap irregularities that are present in experimental practice.

 figure: Fig. 1.

Fig. 1. Schematic view of the spatial positioning of the tomographic measurements. The measurements have been captured such that the volumes are located on a regular grid with an overlap, which is assumed for the presented baseline correction method. The averaged RI maps visualized below are utilized by the BC procedure. They are calculated by taking the mean along Z axis of each volume.

Download Full Size | PDF

3.2 Iterative averaged baseline correction

Since the BC in 3D QPI requires finding an additive constant offset term, there is no need to consider detailed 3D RI volume. Thus, the preprocessing step is introduced to reduce the 3D problem to the 2D one. This is done through averaging of the 3D RI along the Z axis (perpendicular to the stitching plane, see Fig. 1). This operation vastly decreases computational load, reduces impact of local defects in the reconstruction and enables considering any 2D-BC algorithm as universal BC approach. As such it is an important contribution of this work. This step is performed once for each tomographic reconstruction.

Below we provide the details of the novel BC method introduced in the context of 3D QPI. We call the method iterative averaged baseline correction (IABC). Given the 2D RI maps the iterative part of the method is carried out: in each pass of the algorithm over all 2D frames the mean difference between the given image’s margins and all its overlapping counterparts is calculated, averaged and by the end of the iteration added to the current offset of a given frame (offset initialized as zero). Therefore each frame is corrected based on its neighbourhood as shown in Fig. 2(a)). By averaging the contribution of the neighboring frames we avoid the problem of error accumulation [12]. The addition of such averaged difference for all images is a step towards minimizing the overall RI differences between the overlapping regions.

 figure: Fig. 2.

Fig. 2. a) Scheme showing the number of neighboring frames relative to the location of a given frame from the boundary, b) Visualization of the relative placement of the images with notations used in Eq. (3).

Download Full Size | PDF

The images which represent the averaged 2D RI maps are arranged in a regular grid of $k$ rows and $l$ columns, according to the acquisition grid assigned to the motorized stage (see Fig. 2(b)). For all used variables the subscript indicates the location within the grid and the superscript indicates the current iteration number. Let $I_{i, j}$ be the image from $i^{\text {th}}$ row and $j^{\text {th}}$ column. Additionally let us designate the modified version of $I_{i, j}$ in $n^{\text {th}}$ iteration as $\widehat {I^n_{i, j}} = I_{i, j} + P^n_{i, j}$, where $P^n$ represents a set of individual offset values for all corrected images at $n^{\text {th}}$ iteration step. The side margins of $\widehat {I^n_{i, j}}$ are designated as $B^n_{i, j}, T^n_{i, j}, L^n_{i, j} \text { and } R^n_{i, j}$, respectively from bottom, top, left or right side of the image. Those margins are clipped from full images according to nominal overlap percentage. With this convention the RI correction step for each image is described as

$$\begin{aligned} \Delta P^n_{i, j} &= \frac{1}{2} \langle E_{B^n_{i, j}} + E_{T^n_{i, j}} + E_{L^n_{i, j}} + E_{R^n_{i, j}} \rangle\\ P^{n + 1}_{i, j} &= P^{n}_{i, j} + \Delta P^n_{i, j} \end{aligned}$$
with
$$\begin{aligned} E_{B^n_{i, j}} &= \langle T^n_{i+1, j} - B^n_{i, j} \rangle\\ E_{T^n_{i, j}} &= \langle B^n_{i-1, j} - T^n_{i, j} \rangle\\ E_{L^n_{i, j}} &= \langle R^n_{i, j-1} - L^n_{i, j} \rangle\\ E_{R^n_{i, j}} &= \langle L^n_{i, j+1} - R^n_{i, j} \rangle \end{aligned}$$
where $\left\langle \Box \right\rangle$ indicates the average value. $P^n_{i, j}$ is initialized as a set of zeros at $n=0$ and updated with $\Delta P^n_{i, j}$. $\Delta P^n_{i, j}$ is comprised of error terms $E$ that are the average difference between the overlapping margins. We also define the RI error at boundaries as zero, for example $E_{T^n_{0, j}} = 0$ (see Fig. 2(a)). The factor $\frac {1}{2}$ originates from halving the compensation step in order to equally distribute the discrepancies between neighbours, considering that for example $E_{B^n_{i, j}} = -E_{T^n_{i+1, j}}$. This also ensures that the algorithm converges – with a factor higher than $\frac {1}{2}$ it would be possible for $\Delta P^n_{i, j}$ to become unstable throughout the iterations. With each pass the differences between overlapping regions take a step towards minimum in terms of the cost function defined as:
$$f^n(P) = \sum_{i=1}^{k} \sum_{j=1}^{l} \left| \Delta P^n_{i, j} \right|$$
It is due to the fact that each iteration $n$ shifts individual $\Delta P^{n+1}_{i, j}$ towards zero. Although the cost function as a whole is not required to be computed throughout the iterations, tracking its value is a suitable tool for the stopping criterion. Error terms $E$ in Eq. (3) are direct estimates of the baseline error to be corrected. Given that the acquisition setup is well adjusted the small shift between the overlapping regions does not pose a problem. The discrepancy due to misaligned image features cancels out – positive error from one margin cancels out negative error from the second one. As long as the misalignment does not introduce a significant amount of new features unique to only one of the overlapping margins the difference between the images’ margins is an accurate estimation of the baseline error.

It should be noted that the offset values additively influence the error terms $E_{B^n_{i, j}}$, $E_{T^n_{i, j}}$, $E_{L^n_{i, j}}$ and $E_{R^n_{i, j}}$. For example:

$$E_{ab} = \langle \left( I_a + P_a \right) - \left( I_b + P_b \right) \rangle{=} \langle I_a - I_b \rangle{+} P_a - P_b$$
with: $I_a$ and $I_b$ – arbitrary images, $P_a$ and $P_b$ – arbitrary offset values. Hence, the terms may be evaluated only once and updated in each step without explicit recalculation according to Eq. (3). This may be done simply by adding and subtracting respective offset values $P^n$ to the error terms.

The procedure may be easily extended by utilizing the diagonally overlapping regions or by performing image registration of the overlapping margins in order to separately align each overlapping region. The former would utilize all available information, whereas the latter would ensure the alignment of the compared structures present within the margins. Both extensions may be integrated by applying appropriate weights in the calculation of $\Delta P_{i, j}$, such as the nominal overlap percentage in case of diagonal overlaps or relative increase or decrease of the overlapping region due to image alignment procedure.

Those extensions are however beyond the scope of this work.

3.3 Performance analysis in the simulated measurement

The performance analysis requires the ground truth which is not available experimentally. Unfortunately there are no known large phantoms that span over multiple FoVs. For this reason the performance analysis of the proposed method is assessed using the simulated data. Here we will compare the new IABC method against opt-BC, which we have introduced in our previous work [12] in the context of 2D QPI. opt-BC is an L2 minimization of the discrepancy between all overlapping margins using the Broyden-Fletcher-Goldfarb-Shanno optimization procedure with default parameters suggested by the SciPy package [26]. The comparison is made in terms of the relevant error metrics described below and in terms of the execution time.

The simulated scenario consists of a tomographic reconstruction of a large 3D sample comprised of spherical objects suspended in a medium. In the first step, the average 2D RI map is generated from the reconstruction. The resulting 2D map is then split into a $6 \times 6$ grid of images spatially located on a regular grid, with 10% overlap. Each image is $259 \times 295$ pixels. Those images constitute the ground truth to which the set of corrupted-and-corrected images are compared.

The maps are corrupted with four types of errors: (1) systematic error term (same for all images), (2) an individually randomized linear term, (3) an individually randomized offset to simulate the essential problem discussed in this work and (4) a individually randomized low-frequency error term. The justification for adding the first two error terms is that such corrections were required for the successful correction of the experimental data shown in Sec. 4. Error (4) is not corrected and the reason for implementing it is so that the overlapping parts of the images are not perfectly identical – this is to ensure that even if errors (1) and (2) are perfectly estimated there will be some error left to test the robustness of error (3) estimation in the presence of other errors that are always present in practice. The errors (1), (2) and (4) were both modelled as random 2D Legendre polynomials.

Because the results of the experiment depend on the relative magnitudes of modeled signal and disturbances we disclose below the magnitudes of those terms given that the RI contrast of a single sphere is normalized, $\Delta \text {RI} = 1$. For errors (1) and (4) the Legendre coefficients were uniformly randomized in range $-0.5$ to $0.5$ for up to the second degree. For error (2) we used the same range, but limited the degree to one, that is two orthogonal linear terms. The error (3) has been randomized uniformly in range from $0$ to $100$.

The errors (1)-(3) were corrected in the processing pipeline proposed in our previous work [12]. For the correction of error (1) we used the filtered averaged multiple exposure method that averages all available 2D images and estimates the systematic error as the low-pass-filtered version of such an average image. To correct error (2) we used the gradient-based plane fitting (grad-PF) that calculates the gradient of the image and estimates the slope of the linear trend by taking the median value of individual partial derivatives. Finally with the correction of offsets errors (3) the proposed IABC method may be compared with the previously developed opt-BC.

The experiment is divided into two cases. In the first case the images are processed like described above. This scenario represents the typical usage case. Second one is the distorted case, where after correcting the errors (1) and (2) an additional random linear trend is added to a single image within the grid. This case represents a scenario where the linear trend removal failed and the BC method needs to perform in such disturbed condition. The aim of this scenario is to test the robustness of the BC methods.

The discrepancies between the ground truth and the results after corrections are expressed using the following error metrics: $A$, is the root mean squared error (RMSE) of the averaged refractive index with respect to the mounting medium within the objects’ regions $\Delta \mathit {RI}_{\text {avg}}$, divided by its maximum value. Such an error summarizes the local RI error to be expected at any given point of the stitched FoV.

$$A = \frac{\mathrm{RMSE} \left(\Delta \mathit{RI}_{\text{avg}} - \mathbb{H}\left\{\Delta \widetilde{\mathit{RI}}_{\text{avg}}\right\} \right)} {\max(\Delta \mathit{RI}_{\text{avg}})}$$
with $\mathbb {H}\left \{ \Box \right \}$ indicating the operator of a given set of correction methods and $\widetilde {\mathit {RI}}_{\text {avg}}$ being the disturbed version of the ground truth $\mathit {RI}_{\text {avg}}$. The second metric, $B$, is the absolute value of the relative cumulative $\Delta \mathit {RI}_{\text {avg}}$ error.
$$B = \left| \frac{\sum \Delta \mathit{RI}_{\text{avg}} - \sum \mathbb{H}\left\{\Delta \widetilde{\mathit{RI}}_{\text{avg}}\right\}}{\sum \Delta \mathit{RI}_{\text{avg}}} \right|$$
Such a metric estimates an overall RI error to be expected within the stitched FoV. Such an error further propagates to the biomedically relevant metrics.

The randomized simulated experiment is repeated 50 times, for both regular and distorted cases. Metrics A and B are calculated based on the simulated ground truth after 20, 100 and 200 iterations, for each of 50 repetitions. In Table 1 we show an average and standard deviation within the 50 repetitions in terms of both scores, given in percentages. Additionally the average execution time of each approach is given. Note that the execution time depends on the grid size and in the case of opt-BC also on the number of pixels of the overlapping margins. Due to Eq. (5) in IABC the number of pixels influences only the execution time of the first iteration, and in practice is not relevant due to the speed of pixel-wise matrix subtraction and in the context of multiple iterations.

Tables Icon

Table 1. Comparison of the performance of opt-BC and IABC methods in terms of A and B performance metrics and execution time. The methods were applied to 50 randomized simulated datasets, hence the results are presented as an average with a standard deviation.

The results show that the IABC method performs on par with the opt-BC method given enough iterations. This is true for the regular case, where the typical errors are successfully corrected prior to BC correction. In the distorted case IABC performed less robustly. Notably, the IABC method computes over 20 times faster, which is a great advantage in cases with a significantly large number of individual measurements, such as the 320 measurements discussed in Sec. 4. The rapid convergence is visualized in Fig. 3, where the error expressed by cost function (Eq. 4) is shown in the logarithmic scale. The plot visualizes the convergence in the regular case. It was not possible to show the convergence of the distorted case on the same plot, as they overlapped very closely. This is because the cost function does not have access to the ground truth.

 figure: Fig. 3.

Fig. 3. Plot visualizing the minimization of the error expressed by the cost function (Eq. (4)) throughout the iterations of the IABC method in the simulated examples (regular case). The error is normalized – divided by the number of simulated tiles and the overlap percentage. As the plot is an average of 50 experiments, notice the narrow 95% confidence interval.

Download Full Size | PDF

4. Experimental investigation

In order to show the significance of proposed processing the volumetric stitching of a liver organoid is presented. The organoid measured in this work spans 425 µm along its most elongated direction by 250 µm along the orthogonal direction. Spatial sampling of the measurement is 0.1216 µm. The object has been immersed in the mounting medium with $\text {RI@0.633nm} = 1.3350$. Although the captured measurement consisted of 320 individual FoVs, the sample occupies a grid of $15 \times 13$ FoVs, while the rest are object-free. Reconstructed volumes are $530 \times 530 \times 800$ px each, with $58\%$ overlap between neighboring frames. The measurements formed a single layer, so the volumes overlapped side by side, but were not stacked in the axial (Z) direction. The phase projections were captured using optical diffraction tomography system based on Mach-Zehnder interferometer developed at Warsaw University of Technology as a technology demonstrator [15]. It is equipped with $1.3NA$ and $x48$ objective, operating under 0.633 µm illumination wavelength.

The measurements have been reconstructed with the direct inversion using Fourier diffraction theorem based algorithm with Rytov approximation [15]. This resulted in errors in some parts of the object due to phase unwrapping procedure, that is required by the Rytov approximation and often fails in case of highly scattering samples. Hence some regions of the organoid lack clarity due to the inconsistencies in the underlying data.

4.1 Processing path

Initially we attempted the correction with only the BC step. This has not provided satisfactory results, both for opt-BC and IABC methods. The problem was that the obtained correction offset map contained a significant amount of nonuniform offsets in the object-free regions (see successful example at Fig. 4(b)). This led us to assume that the volumes contain errors common in 2D QPI, namely the systematic and linear errors. To address this assumption, we preceded the BC with: (a) averaged multiple exposure method for systematic error subtraction and (b) gradient-based trend removal for linear error correction [12]. Those corrections were performed on the 2D RI maps, before passing them to BC, but estimated errors other then offsets were not transferred to the volumetric data.

 figure: Fig. 4.

Fig. 4. Visualization of the processing steps. Note that the overlapping regions have been cropped out for the visual clarity. a) Raw 2D RI maps averaged along the Z axis, tiled in their nominal positions (before stitching), b) offset terms estimated within 200 iterations of IABC, with range of values marked with red lines on the colorbar (offsets colorbar range fixed to match the dynamic range of the RI colorbar), c) the 2D RI maps after baseline correction. d)-e) Stitched volume cross-section (central) after stitching with BigStitcher obtained with volumes d) without any preprocessing, e) after proposed baseline correction. For the animated version of the comparison visit [29] and [30].

Download Full Size | PDF

Note that the baseline correction is performed prior to stitching, hence the 2D RI maps in Fig. 4(a))-c) are shown tiled in their nominal spatial position (as read from the motorized stage controller) with the overlapping regions cropped out.

To visually compare the results we have stitched both raw volumes (Fig. 4(d)) and volumes corrected with 3D-BC method (Fig. 4(e)). The final stitching procedure has been performed using the BigStitcher [9] software. Because BigStitcher does not operate on floating point images, the individual RI volumes have been rescaled and rounded to the unsigned 16-bit integer values before the stitching and rescaled back to RI values afterwards. BigStitcher has been used to:

  • (1) adjust the spatial position of the volumes in order to improve the alignment of the overlapping features of the object in all three dimensions, first using the pairwise phase correlation image registration technique, which provides shifts for each overlapping pair of images [9,27] followed by global optimization based on the pairwise shifts to establish a single shift for each volume [28],
  • (2) smoothly fuse individual volumes into one single volume using a per-pixel weighted average that minimizes boundary artefacts [9] (Fig. 4(d))-e)).

4.2 Results summary

The obtained results show the validity of the approach proposed in section 2. – the correction offsets range from $-0.0081$ to $0.0044$, therefore spanning $0.0126$ (see red marks in Fig. 4(b)), which is a substantial correction considering that the contrast of the measured sample is approximately $0.03$, hence the raw values are off by $\approx 40\%$. Within the objects 3D footprint the RI RMSE reached 0.0029. Comparative visual assessment between raw (Fig. 4(a)) and corrected (Fig. 4(c)) RI maps additionally indicates the significance of the results, by providing superior continuity between individual measurements of the object– the continuity is expected in the object such as the organoid. The corrections established on the 2D RI maps transfer well to the full volume, as shown in comparison of Fig. 4(d)) and e). We additionally prepared supplemenary materials in the form of the following animations:

  • • GIF switching between raw and corrected stitched volumes [29]
  • • Comparison of cross-sections through the volumes [30]

Overall, we argue that the proposed procedure granted correction that is significant to the 3D QPI measurement stitching and is a first step towards bridging the gap between single FoV measurements with background regions and multi-FoV measurements without ubiquitous access to the background regions.

4.3 Source code and data availability

The source code for the 3D-BC procedure is available here [31]. The raw reconstructed volumes used for this work, together with their spatial placement configuration file are available [32].

5. Conclusions

In this work we discussed the implications of 3D QPI relying on the relative measurements in the context of fragment by fragment acquisition. We argue that measurements of certain objects that do not provide an opportunity to capture object-free region within a single FoV require baseline correction based on the overlapping regions. Furthermore, we reduce the problem of such correction to the 2D case, such that the state-of-the-art 2D baseline correction algorithms can be utilized.

We also present the novel algorithm for 2D baseline correction, together with the quantitative analysis based on simulated dataset. The results show that the performance of our approach is on a par with the SotA in cases where the aberrations of the individual measurements are well corrected. However, it is subpar to the previous solution in the distorted cases. Notably, the novel solution converges much faster, making it vastly competitive in terms of computation time.

Finally, we present the experimental comparison of the tomographic reconstructions of raw measurements and baseline-corrected measurements. The results show the validity of the claims made throughout the work, namely that the baseline unification is a necessary preprocessing procedure for the stitching of 3D RI distributions of large and dense objects and that the problem may be efficiently resolved by reducing the correction to the case similar to the baseline correction of 2D QPI measurements.

Funding

H2020 Industrial Leadership (101016726); Fundacja na rzecz Nauki Polskiej (TEAM-TECH/2016/1/4).

Acknowledgments

CRediT statement: Piotr Stępień: Conceptualization, Formal analysis, Methodology, Software, Visualization, Writing – original draft. Michał Ziemczonok: Data curation, Investigation, Writing – review & editing. Małgorzata Kujawińska: Funding acquisition, Supervision, Formal analysis, Writing – review & editing. Maria Baczewska: Investigation, Resources. Luca Valenti, Alessandro Cherubini, Elia Casirati: Resources. Wojciech Krauze: Supervision, Methodology, Project administration, Writing – review & editing.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

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. K. Lee, K. Kim, J. Jung, J. Heo, S. Cho, S. Lee, G. Chang, Y. Jo, H. Park, and Y. Park, “Quantitative phase imaging techniques for the study of cell pathophysiology: from principles to applications,” Sensors 13(4), 4170–4191 (2013). [CrossRef]  

2. Y. Park, C. Depeursinge, and G. Popescu, “Quantitative phase imaging in biomedicine,” Nat. Photonics 12(10), 578–589 (2018). [CrossRef]  

3. Tomocube Inc., “Tomocube,” http://www.tomocube.com/ (2020). [Online; accessed 24-May-2022].

4. Nanolive SA, “Nanolive,” https://nanolive.ch (2020). [Online; accessed 24-May-2022].

5. V. Balasubramani, A. Kuś, H.-Y. Tu, C.-J. Cheng, M. Baczewska, W. Krauze, and M. Kujawińska, “Holographic tomography: techniques and biomedical applications,” Appl. Opt. 60(10), B65–B80 (2021). [CrossRef]  

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

7. A. J. Lee, H. Hugonnet, W. Park, and Y. Park, “Three-dimensional label-free imaging and quantification of migrating cells during wound healing,” Biomed. Opt. Express 11(12), 6812–6824 (2020). [CrossRef]  

8. H. Hugonnet, Y. W. Kim, M. Lee, S. Shin, R. H. Hruban, S.-M. Hong, and Y. Park, “Multiscale label-free volumetric holographic histopathology of thick-tissue slides with subcellular resolution,” Adv. Photonics 3(02), 026004 (2021). [CrossRef]  

9. D. Hörl, F. Rojas Rusak, F. Preusser, P. Tillberg, N. Randel, R. K. Chhetri, A. Cardona, P. J. Keller, H. Harz, H. Leonhardt, M. Treier, and S. Preibisch, “Bigstitcher: reconstructing high-resolution image datasets of cleared and expanded samples,” Nat. Methods 16(9), 870–874 (2019). [CrossRef]  

10. A. C. Kak, M. Slaney, and G. Wang, “Principles of computerized tomographic imaging,” (2002).

11. P. Stępień, D. Korbuszewski, and M. Kujawińska, “Digital holographic microscopy with extended field of view using tool for generic image stitching,” ETRI J. 41(1), 73–83 (2019). [CrossRef]  

12. P. Stępień, W. Krauze, and M. Kujawińska, “Preprocessing methods for quantitative phase image stitching,” Biomed. Opt. Express 13(1), 1–13 (2022). [CrossRef]  

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

14. S. O. Isikman, W. Bishara, S. Mavandadi, W. Y. Frank, S. Feng, R. Lau, and A. Ozcan, “Lens-free optical tomographic microscope with a large imaging volume on a chip,” Proc. Natl. Acad. Sci. 108(18), 7296–7301 (2011). [CrossRef]  

15. A. Kuś, W. Krauze, P. L. Makowski, and M. Kujawińska, “Holographic tomography: hardware and software solutions for 3d quantitative biomedical imaging,” ETRI J. 41(1), 61–72 (2019). [CrossRef]  

16. V. Lauer, “New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope,” J. Microsc. 205(2), 165–176 (2002). [CrossRef]  

17. S. J. LaRoque, E. Y. Sidky, and X. Pan, “Accurate image reconstruction from few-view and limited-angle data in diffraction tomography,” J. Opt. Soc. Am. A 25(7), 1772–1782 (2008). [CrossRef]  

18. 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(01), 1 (2013). [CrossRef]  

19. W. Krauze, P. Makowski, M. Kujawińska, and A. Kuś, “Generalized total variation iterative constraint strategy in limited angle optical diffraction tomography,” Opt. Express 24(5), 4924–4936 (2016). [CrossRef]  

20. J. Lim, A. B. Ayoub, E. E. Antoine, and D. Psaltis, “High-fidelity optical diffraction tomography of multiple scattering samples,” Light: Sci. Appl. 8(1), 1–12 (2019). [CrossRef]  

21. W. Krauze, “Optical diffraction tomography with finite object support for the minimization of missing cone artifacts,” Biomed. Opt. Express 11(4), 1919–1926 (2020). [CrossRef]  

22. L. Tian, J. Wang, and L. Waller, “3d differential phase-contrast microscopy with computational illumination using an led array,” Opt. Lett. 39(5), 1326–1329 (2014). [CrossRef]  

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

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

25. M. Chen, D. Ren, H.-Y. Liu, S. Chowdhury, and L. Waller, “Multi-layer born multiple-scattering model for 3d phase microscopy,” Optica 7(5), 394–403 (2020). [CrossRef]  

26. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, and P. van Mulbregt, and SciPy 1.0 Contributor, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nat. Methods 17(3), 261–272 (2020). [CrossRef]  

27. C. D. Kuglin, “The phase correlation image alignment methed,” in Proc. Int. Conference Cybernetics Society, (1975), pp. 163–165.

28. S. Preibisch, S. Saalfeld, and P. Tomancak, “Globally optimal stitching of tiled 3d microscopic image acquisitions,” Bioinformatics 25(11), 1463–1465 (2009). [CrossRef]  

29. P. Stępień, M. Ziemczonok, M. Kujawińska, M. Baczewska, L. Valenti, A. Cherubini, E. Casirati, and W. Krauze, “Numerical refractive index correction for the stitching procedure in tomographic quantitative phase imaging – supplementary animation,” figshare, 2022, https://doi.org/10.6084/m9.figshare.19650579 .

30. P. Stępień, M. Ziemczonok, M. Kujawińska, M. Baczewska, L. Valenti, A. Cherubini, E. Casirati, and W. Krauze, “Numerical refractive index correction for the stitching procedure in tomographic quantitative phase imaging – cross-section animations,” figshare2022, https://doi.org/10.6084/m9.figshare.20440440.

31. P. Stępień, “QPI-stitching-2D-3D,” Github, 2022, https://github.com/biopto/QPI-stitching-2D-3D.

32. P. Stępień, M. Ziemczonok, M. Kujawińska, M. Baczewska, L. Valenti, A. Cherubini, E. Casirati, and W. Krauze, “Numerical refractive index correction for the stitching procedure in tomographic quantitative phase imaging – dataset,” Zenodo, 2022, https://doi.org/10.5281/zenodo.6977026.

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 (4)

Fig. 1.
Fig. 1. Schematic view of the spatial positioning of the tomographic measurements. The measurements have been captured such that the volumes are located on a regular grid with an overlap, which is assumed for the presented baseline correction method. The averaged RI maps visualized below are utilized by the BC procedure. They are calculated by taking the mean along Z axis of each volume.
Fig. 2.
Fig. 2. a) Scheme showing the number of neighboring frames relative to the location of a given frame from the boundary, b) Visualization of the relative placement of the images with notations used in Eq. (3).
Fig. 3.
Fig. 3. Plot visualizing the minimization of the error expressed by the cost function (Eq. (4)) throughout the iterations of the IABC method in the simulated examples (regular case). The error is normalized – divided by the number of simulated tiles and the overlap percentage. As the plot is an average of 50 experiments, notice the narrow 95% confidence interval.
Fig. 4.
Fig. 4. Visualization of the processing steps. Note that the overlapping regions have been cropped out for the visual clarity. a) Raw 2D RI maps averaged along the Z axis, tiled in their nominal positions (before stitching), b) offset terms estimated within 200 iterations of IABC, with range of values marked with red lines on the colorbar (offsets colorbar range fixed to match the dynamic range of the RI colorbar), c) the 2D RI maps after baseline correction. d)-e) Stitched volume cross-section (central) after stitching with BigStitcher obtained with volumes d) without any preprocessing, e) after proposed baseline correction. For the animated version of the comparison visit [29] and [30].

Tables (1)

Tables Icon

Table 1. Comparison of the performance of opt-BC and IABC methods in terms of A and B performance metrics and execution time. The methods were applied to 50 randomized simulated datasets, hence the results are presented as an average with a standard deviation.

Equations (7)

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

Φ ( x , y ) = L Δ n ( x , y , z ) d l , t h e r e f o r e , Φ ( x , y ) + C = L ( Δ n ( x , y , z ) + C ) d l ,
Δ P i , j n = 1 2 E B i , j n + E T i , j n + E L i , j n + E R i , j n P i , j n + 1 = P i , j n + Δ P i , j n
E B i , j n = T i + 1 , j n B i , j n E T i , j n = B i 1 , j n T i , j n E L i , j n = R i , j 1 n L i , j n E R i , j n = L i , j + 1 n R i , j n
f n ( P ) = i = 1 k j = 1 l | Δ P i , j n |
E a b = ( I a + P a ) ( I b + P b ) = I a I b + P a P b
A = R M S E ( Δ R I avg H { Δ R I ~ avg } ) max ( Δ R I avg )
B = | Δ R I avg H { Δ R I ~ avg } Δ R I avg |
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.