- Email: [email protected]
Optical projection tomography via phase retrieval algorithms
Methods xxx (2017) xxx–xxx
Contents lists available at ScienceDirect
Methods journal homepage: www.elsevier.com/locate/ymeth
Optical projection tomography via phase retrieval algorithms Daniele Ancora a,b, Diego Di Battista a,c, Georgia Giasafaki a, Stylianos E. Psycharakis a,d, Evangelos Liapis a, Jorge Ripoll e,f, Giannis Zacharakis a,⇑ a
Institute of Electronic Structure and Laser, Foundation for Research and Technology – Hellas, 70013 Heraklion, Greece Department of Materials Science and Technology, University of Crete, 71003 Heraklion, Greece c Assing S.p.A, Monterotondo, 00015 Rome, Italy d School of Medicine, University of Crete, 71003 Heraklion, Greece e Department of Bioengineering and Aerospace Engineering, Universidad Carlos III de Madrid, 28911 Madrid, Spain f Instituto de Investigación Sanitaria del Hospital Gregorio Marañón, 28007 Madrid, Spain b
a r t i c l e
i n f o
Article history: Received 7 August 2017 Received in revised form 15 October 2017 Accepted 17 October 2017 Available online xxxx Keywords: Optical projection tomography Phase retrieval Registration-free imaging Autocorrelation imaging Hidden object tomography
a b s t r a c t We describe a computational method for accurate, quantitative tomographic reconstructions in Optical Projection Tomography, based on phase retrieval algorithms. Our method overcomes limitations imposed by light scattering in opaque tissue samples under the memory effect regime, as well as reduces artifacts due to mechanical movements, misalignments or vibrations. We make use of Gerchberg–Saxton algorithms, calculating first the autocorrelation of the object and then retrieving the associated phase under four numerically simulated measurement conditions. By approaching the task in such a way, we avoid the projection alignment procedure, exploiting the fact that the autocorrelation sinogram is always aligned and centered. We thus propose two new, projection-based, tomographic imaging flowcharts that allow registration-free imaging of opaque biological specimens and unlock three-dimensional tomographic imaging of hidden objects. Two main reconstruction approaches are discussed in the text, focusing on their efficiency in the tomographic retrieval and discussing their applicability under four different numerical experiments. Ó 2017 Elsevier Inc. All rights reserved.
1. Introduction Optical biomedical imaging modalities have become essential tools for disease detection and monitoring, occupying an increasing portion of medical imaging needs [1], characterized by the use of non-ionizing radiation, molecular specific contrast agents and minimum or non-invasiveness. Among these methods, Optical Projection Tomography (OPT) [2], considered the optical analogue of X-ray Computed Tomography (XCT), shares with the latter a broad set of reconstruction approaches, while the measurement scheme might change slightly depending on the experimental architecture. OPT relies on mechanical movement of the specimen rotating it around a fixed axis to acquire the full set of projections, rather than moving the measuring device around the object. However, mechanical translation and rotations are then needed to perform a three-dimensional reconstruction of the sample of interest, thus requiring a correct calibration of the experiment to obtain accurate reconstructions. Such calibration has to be often performed both at hardware and software level and might sometimes limit the acqui⇑ Corresponding author. E-mail address: [email protected] (G. Zacharakis).
sition procedure or introduce unwanted artifacts in the reconstructions. Many efforts have been made in finding new algorithmic ways to align, center or reduce misalignment effects at the tomographic level [3–6], but none has the advantage of being fully alignment-free. In this article, we study a more generic approach to confront with misaligned data reconstructions, proposing two algorithmic step-wise protocols that are able to retrieve quantitative reconstructions even in the presence of constant drifts, random vibrations or even overcoming the phase scrambling due to the presence of a scattering layer enclosing the object. The study focuses on two different approaches, both of them relying on the use of phase retrieval algorithms, to generate aligned autocorrelation sinograms that can be inverted in order to obtain tomographic reconstructions. We emphasize the fact that any positional information of the specimen is not important, since the autocorrelation is inherently always centered in its space of definition. Phase-Retrieved Tomography (PRT) was originally introduced in our previous work [7] to image experimentally the fluorescence distribution of necrotic cells in a human-breast tumor spheroid [8]. In this case the method was compared against classical reconstruction schemes from Single Plane Illumination Microscopy (SPIM) [9,10] and OPT, showing promising reconstruction abilities. Here
https://doi.org/10.1016/j.ymeth.2017.10.009 1046-2023/Ó 2017 Elsevier Inc. All rights reserved.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
2
D. Ancora et al. / Methods xxx (2017) xxx–xxx
we focus on a broader numerical study, better illustrating the potential of phase retrieval techniques applied for tomographic investigations. Alongside PRT, we propose another iterative technique to recover aligned sinograms to which we refer to as iterative Projection-Retrieved Tomography (iPRT). The two methods will be tested for the first time in reconstructing corrupted datasets generated by numerical simulations, focusing in particular on their ability to leave the user free of any alignment procedure. The possibility of performing hidden reconstructions will be further considered showing that, under memory effect conditions, blind tomography becomes naturally possible as a result of the incorporation of phase retrieval techniques into classical reconstruction schemes. 2. Material and methods In this article we generate the datasets numerically, analyzing the possibility offered by the methods under different measurement scenarios in order to have a known reference to compare with. The object under investigation is a previously retrieved reconstruction of the fluorescence distribution in a human-breast tumor spheroid [7], which we will use as a ground truth measurement. This volume has a size of 300 pixels (from now on abbreviated as px) in each spatial dimension, being characterized in the original experiment by a pixel size of 0:8 lm. Starting with this volume, we generate several datasets representing four different experimental measurement conditions. The first dataset is created by rotating the object around an axis parallel to z and positioned at the center of the camera plane, thus acquiring its intensity projections (Fig. 1, row A) as a function of the observation angle. This mathematically corresponds to the Radon transform of the sample volume and, experimentally, to an ideal OPT experiment. The second dataset (row B) is obtained by introducing a random twodimensional slight vibrations of 1 px during the rotation, leading
to the generation of a fuzzy sinogram. This could represent a situation in which we try to correct the axis of rotation to fit it onto the position of the barycenter of the object, or the case in which the stages slightly vibrate during the rotation. The third dataset (row C) is equal to the above except for the introduction of a constant drift in two dimensions while rotating the object. In fact, in some experimental acquisitions [7], we observed a similar perturbation due to the effect of gravity (sample slowly falling inside the medium) and due to mechanical misalignment of the rotating stage. The fourth, and last, dataset (row D) reproduces a measurement performed when the object is hidden behind a phase scrambling environment [11,12], such as a scattering envelope (e.g. an egg or a cocoon). This dataset is numerically calculated assuming the measurements to be in the memory effect regime [13], which states that a speckle pattern (generated by a single point source) translates in relation to the source movement behind the phasescrambling layer [14]. All the above described numerical experiments produced different sinograms (Fig. 1, column 2), that can be inverted (with an inverse Radon transformation) to obtain the volume reconstruction in each specific case. Column 3 of Fig. 1 shows, depending on the perturbation introduced in the measurement, how the reconstructed volume differs from the ideal case ultimately producing several kinds of unwanted reconstruction artifacts well known in the literature [4,15]. In the following, we will try to avoid the generation of alignment perturbations, exploiting the mathematical properties of the autocorrelation in combination with the phaseretrieval algorithm. 2.1. Autocorrelation sinogram From each of the previously obtained sinograms we can calculate the corresponding autocorrelation sinogram Ah (A-sinogram), which will be used for both the phase-retrieval methods. In this
Fig. 1. Four different numerical experiments, showing the resulting sinograms (column 2) of the three different rotations. A) The results of an aligned rotation, leading to optimal reconstructions of the object. B) A one-pixel random vibration was introduced while rotating the object, returning slight artifacts in the reconstruction. C) The measurement was affected by a two-dimensional drift (gravity falling and axis displacement) particularly evident in the corresponding reconstruction. D) Hidden object measurement, impossible to backproject with classical techniques.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
D. Ancora et al. / Methods xxx (2017) xxx–xxx
notation, h is the angle at which we acquire the projections that spans through the entire 360° rotation in steps of 1°. The calculation of the A-sinogram is simple for the direct projections and can be estimated by using the Wiener-Kinchin theorem:
Ah ¼ F 1 fjF fPh gj2 g:
ð1Þ
Instead, to obtain the A-sinogram in the case of hidden measurements we need to average the autocorrelation over different speckle-sinograms obtained by using different random phases, in analogy to the ensemble averages requested by the corresponding two-dimensional technique [11]. In our specific case, we averaged over 10 different random phases, but in general a higher number of phases corresponds to a better sampling of the object autocorrelation. In Fig. 2 column 1 (row A), we display the result of the calculation of the A-sinogram in the case of direct measurements (described in Fig. 1, rows A-C), which results to be identical in the three cases examined and then insensitive to any drift or vibration of the sample during the rotation. Fig. 2 column 1 (row C) shows the A-sinogram calculated from the hidden measurement (Fig. 1, row D), which exhibits the same features of the one calculated with direct measurements. Moreover, the A-sinogram resulting from the hidden experiment gets more accurate when we increase the ensemble average, and the noise (dark background of the A-sinogram in column 1, row C) can be further suppressed with a more accurate sampling. 2.2. Iterative projection retrieval tomography (iPRT) A straightforward approach, once we have calculated the Asinogram, is to proceed stepwise aiming at the retrieval of every single projection, thus rebuilding an aligned sinogram. To do so, we firstly retrieve the phase connected with the 2D autocorrelation of the initial projection, e.g. at angle h ¼ 0, with a classical phase
3
retrieval (PR) algorithm [16] starting with a random-phase guess. The position of the object retrieved in this way will be located randomly in space, due to the fact that the autocorrelation operation is invariant for any object translation. Running independent PR per each projection would return an unpredictable object shift in the camera plane, rendering a full reconstruction impossible. However, if we assume that the object is rotating in small steps, we could use the obtained phase at the starting angle /ðhÞ to guide the PR reconstruction of the next projection at angle h þ 1. In this way, the PR problem converges faster, since we already start with a guess that is very close to the phase solution /ðh þ 1Þ. Iterating this process until the last projection N and stacking every PR reconstruction with the same order of the A-sinogram, we end up having a seemingly aligned sinogram. Since the sinogram was build starting from a reconstruction located in a random spatial position, the rotation axis of such dataset is in general not centered in the camera plane, but it can be automatically found by comparing the reconstruction at 0° and 180°. Once centered, the obtained sinogram can be backprojected (e.g. with a filtered inverse Radon transformation) to retrieve the final reconstruction of the specimen. The whole process is described in the flowchart A of Fig. 3 and we refer to it as iterative Projection-Retrieved Tomography (iPRT). It is important to notice that we do not strictly have to insert the first projection into the iPRT in case we have direct access to the object image. In fact, by using the phase of the initial projection as the starting point of iPRT we will be able to align the sinogram with respect to the initial measurement, having no need to retrieve once again the phase of the initial step. Although this implies a further degree of simplicity (we do not have to run the first PR iterations), for the purpose of this study we will still force the retrieval of the phase of the first projection. We consider this situation to stress the fact that the overall phase is not a necessary aspect for realigning projections and applying iPRT directly to hidden measurements, where the initial object projection is not directly accessible.
Fig. 2. Autocorrelation sinograms (column 1) obtained for the four numerical experiments and their consequent iPRT sinogram reconstruction. A) PRT imaging for any of the measurement scheme (perfect rotation, vibration and drift), in column 2) it is shown the final three-dimensional object autocorrelation and in 3) the final reconstruction. B) iPRT volume reconstruction of the same datasets. C) Hidden PRT reconstruction, showing a nearly similar autocorrelation compared to that of A). D) Hidden iPRT reconstruction.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
4
D. Ancora et al. / Methods xxx (2017) xxx–xxx
Fig. 3. Block diagram of the phase retrieval methods for OPT reconstructions. A) The iPRT flowchart while B) The simpler PRT flowchart. Although iPRT has a more complicated flowchart, the method is less computational demanding, due to the fact that it relies on 2D Fourier transformations rather than 3D.
2.3. Phase-retrieved tomography (PRT) A different approach for the alignment of the dataset can be implemented by backprojecting directly the A-sinogram that was previously calculated. Due to the fact that the A-sinogram is perfectly centered, the volume reconstructed in this way will not suffer any misalignment, while it can be proven that is mathematically equivalent to the three-dimensional autocorrelation of the whole specimen [7]. The phase retrieval problem, then, can be extended to three dimensions by starting with a random 3Dphase, thus its solution would directly give the reconstruction of the object of interest with no further need of alignment procedures. Hence, we refer to this method as Phase-Retrieved Tomography (PRT) and its corresponding flowchart is shown in panel B of Fig. 3. We point out that although the flowchart is simpler than that of iPRT, the PRT method is computationally more demanding, since the phase retrieval algorithm employed for the solution is based on three-dimensional Fourier transforms, rather than twodimensional. 2.4. Phase retrieval algorithms Both of the methods proposed, base their actions on the usage of Gerchberg–Saxton phase retrieval (PR) algorithms [17] in order to retrieve the phase connected with the autocorrelation, or more appropriately, with the estimation of the modulus of the Fourier transform jUðkÞj. The retrieval of the correct phase, then, enables the reconstruction of either the single projections (iPRT) or the whole object (PRT). For both of the methods, we used a common combination of the hybrid input-output (HIO) and the error reduction (ER) [18] algorithms, starting with 2000 iteration of the HIO followed by 100 iterations of the ER. Both algorithms fall into the class of alternating projection methods [19,16], in which we pass
from the real (or object) space to the Fourier space applying appropriate constraining operations. The two approaches follow a simple sequence of operations composed by four steps, in which only the last one changes according to the algorithm used. Without losing in generality, in the following we will treat the one-dimensional case that can be trivially extended to two- or three-dimensions, respectively implemented into iPRT and PRT reconstruction protocols. As a first step, the PR algorithm starts with an initial random guess of the object g 0 ðxÞ then, at the j-th iteration, the algorithm acts with this sequence of operation: 1. Fourier transform the object estimation:
g j ðxÞ ! Gj ðkÞ ¼ F fg j ðxÞg
ð2Þ
2. Substitute its modulus with the originally measured one:
Gj ðkÞ ¼ jGj ðkÞjei/ðkÞ ! G0j ðkÞ ¼ jUðkÞjei/ðkÞ
ð3Þ
3. Calculate the inverse Fourier transform of this estimation: 0ðkÞ
Gj
! g 0j ðxÞ ¼ F 1 fG0j ðkÞg
ð4Þ
4. For the HIO, calculate the object update by using a feedback coming from the previous iterations in C, the region where the assumption of the image being real and positive are violated: (
for x R C
g 0j ðxÞ
g jþ1 ðxÞ ¼
g j ðxÞ bg 0j ðxÞ for x 2 C
ð5Þ
Or, for the final ER iterations, setting to zero the values in the region C:
(
g jþ1 ðxÞ ¼
g 0j ðxÞ for x R C 0
for x 2 C
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
ð6Þ
5
D. Ancora et al. / Methods xxx (2017) xxx–xxx
Ideally, such process is iterated until the modulus of the Fourier transform of the calculated object jGj ðkÞj converges to the measured one jUðkÞj, or when there are no points left within the region C. In practice, the convergence of the algorithm can be monitored calculating the Recovery Error Function (REF), defined as the Euclidean distance between the modulus of the Fourier transform and the measured one:
REF ¼
X jjGj ðkÞj jUðkÞjj2 :
ð7Þ
k
This quantity converges to zero, REF ! 0, when the solution approaches the exact one. Only for the case of iPRT at the second step of the sinogram reconstruction (h ¼ 1) we use an estimation of the object given by the previous reconstruction (instead of the random guess), and we consequently reduce the number of the steps of the HIO to 100. From now on, we will refer to PR algorithms as the two or three-dimensional processes that retrieve any phase from the corresponding autocorrelations; we will refer to PRT and iPRT as the overall imaging process, involving PR reconstructions at different steps described by the flowcharts of Fig. 3.
Considering only the translational variable g, we notice that it is possible to change the variable of integration to:
Z
Oðx þ g; y þ n; z þ vÞdg ¼
Z
Oðx þ g; y þ n; z þ vÞdx
¼ PO ðy þ n; z þ vÞ:
ð15Þ
In fact, the projection on g is independent on any x shift parallel to its direction. This ultimately implies that we can write:
PO ðy; zÞPO ðy þ n; z þ vÞ ¼
Z
Oðx; y; zÞPO ðy þ n; z þ vÞdx
ð16Þ
where it is now obvious that both the terms are equal, thus proving the initial assumption in Eq. (4). Since this operation can be repeated at each angle of observation, the backprojection of the A-sinogram leads to the reconstruction of the 3D autocorrelation of the specimen. We conclude that the association of the proper phase to this quantity, thus solving the three-dimensional phase retrieval problem, would enable us to retrieve the final reconstruction of the object.
3. Theory and calculation
4. Results and discussion
3.1. PRT theoretical background
In this section, we compare the results obtained for iPRT and PRT reconstruction schemes. We point out that in all retrieval processes studied in the present work, we always obtained consistent reconstructions starting from random guesses, differing only from each other by a random spatial translation. This is due to the fact that the autocorrelation is degenerate for any spatial translation of the sample, which consequently leads to an overall loss of information regarding the whole object positioning in space. Each reconstructed volume, then, was automatically translated back to the position of the ground truth volume by finding the absolute maximum of their cross-correlation. This procedure was necessary in order to compare the results between the different measuring schemes, but it is not required for the reconstruction itself. Fig. 2 shows the entire reconstruction process in a progressive diagram. Column 1 displays the A-sinograms (rows A and C) and the corresponding realigned sinograms calculated from those (row B and D). Column 2 visualizes the rendering of the three-dimensional autocorrelations calculated from the A-sinograms (only for rows A and C). Column 3 shows the final rendering of all the reconstructed volumes for iPRT and PRT. In more detail, Fig. 2 row A describes the PRT method applied to reconstruct misaligned datasets. The A-sinogram of column 1 was calculated from any of the direct sinogram observations of Fig. 1 column 2 (raw A to C), resulting perfectly aligned and identical for all the cases. Such A-sinogram is backprojected (hence, inverse Radon transformed) leading to the calculation of the autocorrelation volume of the specimen of interest, rendered in column 2 for visualization purposes. Subsequently, a three-dimensional phase retrieval algorithm was used to find the phase associated to this volume, enabling the imaging of the object. This finalizes the PRT reconstruction protocol of which the final result is displayed in column 3. In this case, the retrieved fluorophore distribution appears to be accurate and insensitive to any of the perturbations described in Fig. 1 (random vibrations and drifts), seemingly equivalent to the original ground truth (Fig. 1 column 1, row A). On the other hand, illustrated in Fig. 2 row B, the iPRT managed to reconstruct the sinogram from the autocorrelation sequence previously described, without relying on any spatial information of the object rotation. Although the retrieved sinogram (displayed in column 1, row B) seems to rotate smoothly, it contains a certain degree of misalignment and results to the production of artifacts when we backproject it for the final reconstruction. The result of the iPRT protocol is rendered in col-
Regardless of how the A-sinogram was calculated, the PRT method relies on the fact that its backprojection leads to the reconstruction of the three-dimensional autocorrelation of the sample. This can be mathematically verified since the problem can be reduced to the fact that the projection of the 3D autocorrelation must be equal to the 2D autocorrelation of the object projection per every angle of the measurement. Let us define the object of interest as Oðx; y; zÞ and its 3D autocorrelation Aðg; n; vÞ as:
Aðg; n; vÞ ¼ Oðx; y; zÞHOðx; y; zÞ Z ¼ Oðx; y; zÞOðx þ g; y þ n; z þ vÞdx dy dz:
ð8Þ
We define the observation angle h ¼ 0 to be parallel to the x axis (and consequently also to g), thus the projection of the object at this angle results equal to:
Z
PO ðy; zÞ ¼ Px ¼
ð9Þ
Oðx; y; zÞdx
and the projection of the autocorrelation:
PA ðn; vÞ ¼
Z
Aðg; n; vÞdg:
ð10Þ
Due to the above consideration, the problem is reduced to verify that the following relation is satisfied:
PO ðy; zÞHPO ðy; zÞ ¼ PA ðn; vÞ: Z Z
ð11Þ
We can write the terms in Eq. (9) as:
PO ðy; zÞPO ðy þ n; z þ vÞdy dz ¼ PO ðy; zÞPO ðy þ n; z þ vÞdydz ¼
Z
Z
Aðg; n; vÞdg
ð12Þ
Oðx; y; zÞOðx þ g; y þ n; z
þ vÞdx dy dz dg:
ð13Þ
Since in both the integrals we are integrating along y and z, we can then compare the two arguments, leading to the simplified formulation:
PO ðy;zÞPO ðy þ n;z þ vÞ ¼
Z
Oðx; y; zÞOðx þ g; y þ n; z þ vÞdxdg:
ð14Þ
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
6
D. Ancora et al. / Methods xxx (2017) xxx–xxx
umn 3. It is clear, in this case, that PRT performed better than iPRT retrieving a sharper reconstruction, although also iPRT converged into a meaningful (but noisier) result. Fig. 2 rows C and D show respectively the PRT and iPRT for the reconstructions of the hidden object, following exactly the same scheme described for the realignment of direct measurements. Both reconstructions were applied to the dataset of Fig. 1 row D, calculated assuming that the object was secreted by a scattering layer within the memory effect regime. In normal imaging modalities, the acquisition of a seemingly-random speckle pattern did not allow any classical reconstruction methods (volume in Fig. 1 column 3, row D), thus rendering the use of the phase retrieval protocols essential for any tomographic analysis. The A-sinogram (Fig. 2 column 1, raw C) results perfectly aligned and identical to the one obtained in the case of direct measurements, being characterized by a noisier background. The autocorrelations were calculated averaging 10 different random speckle patters, but we verified that a more accurate sampling would push further down the noise in the A-sinogram. In any case, its inversion returned a very good reconstruction of the three-dimensional autocorrelation of the sample (column 2, raw C), visually identical to the one reconstructed from direct measurements. To finalize the imaging process, the phase retrieval finds the phase to connect with such volume, thus leading to a very accurate reconstruction of the sample concealed by a scattering layer (column 3, raw C). From the Asinogram of the hidden measurements, we could generate back the direct sinogram of the object rotation by running the iPRT protocol (column 1, raw D). In this case, the resulting sinogram is very similar to the one obtained by direct measurements and leads to a consistent reconstruction of the specimen (column 3, raw D). In close analogy with the previous case, PRT returned better reconstructions also in the case of blind tomography if compared with iPRT, leading to sharp imaging of the fluorescence distribution. The resulting volume has a noisier background if compared to the PRT for the normal dataset, but still represents an excellent reconstruction for a hidden measurement. Also the iPRT converged into a correct solution, but this time leaving more artifacts if compared to the previous case. To have a more comprehensive understanding about the quality of the reconstruction, it is interesting to compare the projections along the three-spatial dimensions. Thus, calling any final recon-
structed volume as Vðx; y; zÞ, peak-normalized to unity, we calculate:
Z Px ¼
Z Vðx; y; zÞdx; Py ¼
Z Vðx; y; zÞdy; Pz ¼
Vðx; y; zÞdz;
ð17Þ
which lead to the production of three projected images per each reconstruction approach (Fig. 4). We can appreciate that PRT definitely solved the problem of misalignment (panel E), returning a nearly perfect reconstruction that matches that of a perfectly aligned measurement (shown in panel A). Instead, iPRT (panel F) used as a tool to realign the projection sinogram, returned a reconstruction having artifacts similar to the reconstruction perturbed by a constant drift (panel C). On the other hand, also, hidden PRT retrieved a perfectly aligned fluorophore distribution, but the reconstruction lies in a noisier environment if compared to that of the normal PRT (panel G). Lastly, the volume produced by the hidden iPRT reconstruction (panel H) seems to be affected by the same vibrational artifacts typical of a reconstruction in a vibrationalperturbed case (panel B). Fig. 5 confirms this fact, presenting the analysis of the residual misalignment of the reconstruction compared to the original, thus aligned, dataset. Here we use the aligned rotation as reference to measure, as a function of the angle, the position of the perturbed measurements and the reconstructed sinograms obtained with iPRT. Panel A presents the behavior of the displacement with respect to the horizontal axis (X in the sinogram): the residual of the iPRT reconstruction (violet curve) drifts away with respect to the aligned dataset (blue) more similarly to the drifted dataset (yellow) than to the vibrationally perturbed one. The reconstruction of hidden object with iPRT (green), instead, suffers a slight drift but stronger vibrational misalignment, making it more similar to the vibrationally perturbed sinogram (red). The residual alignment of the vertical axis Z (parallel to the rotation axis) shown in panel B presents another interesting situation, in which iPRT reconstruction does not suffer strong drift-effect but still a vibrational error affects the hidden iPRT reconstruction. In synthesis, iPRT reconstruction drifts gently away with respect to the original dataset in the horizontal direction, while there is no appreciable effect on the vertical axis. On the other hand, hidden iPRT clearly suffers vibrational perturbation on both axis thus leading to the visual artifacts discussed in Fig. 4.
Fig. 4. Reconstruction comparison, projection view. A) reconstruction from perfectly aligned measurement, B) vibrationally perturbed rotation, D) two-dimensional constant drift and E) hidden measurements. E-F) Corresponding reconstruction obtained respectively with PRT and iPRT by using the same dataset used for A-B). G-H) reconstructions obtained from hidden dataset.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
D. Ancora et al. / Methods xxx (2017) xxx–xxx
Fig. 5. Residual shifts of the perturbed and iPRT reconstructed sinograms with respect to the aligned datasets. A) Displacement with respect to the horizontal direction X, perpendicular to the rotation axis. B) Displacement with respect to the vertical direction Z, parallel to the axis of rotation.
Qualitatively, we observe that every reconstruction was able to retrieve a correct fluorophore distribution, but in general PRT works far better than iPRT in terms of sharpness and accuracy of the results. In any case, both the techniques overtook the limitations of performing hidden reconstruction (Fig. 4, panel D), opening novel promising imaging scenarios. To quantitatively compare the results, we focus on plot profiles calculated from projections onto two different axes, according to the following calculation:
Z Pxy ¼
Z Vðx; y; zÞdxdy; P yz ¼
Z Vðx; y; zÞdydz; Pzx ¼
Vðx; y; zÞdxdz; ð18Þ
From the analysis of the profiles in Fig. 6, we conclude that phase retrieval methods quantitatively retrieve the intensity distribution and that the results of PRT are as accurate as the ones from an aligned reconstruction. We can also observe the drift-effect introduced by iPRT reconstruction (yellow line, Fig. 6 panel D), similar to the effect of the reconstruction in the presence of drift (violet line, panel A, B). Of particular interest are also the results obtained for hidden PRT (violet line, panel D-F), in which a noisier background lifts up the projection profiles by a constant amount, still leaving the fluorophore distribution highly visible compared to the noise level. Interestingly then, we can conclude the analysis claiming that the PRT method performs better than iPRT for any measuring conditions and overcomes the limitations associated to the alignment of the sample movement. The latter, instead, suffer from two different alignment issues, depending on the kind of measurement initiating our analysis. Iterative reconstruction of subsequent single projections leads to a reconstruction of a sinogram affected
7
by a slight drift, while on the other hand, in the case of speckle measurement, it suffers artifacts typical of the vibrational noise. Finally, it is important to notice that all the methodologies proposed rely on PR algorithms to perform the reconstructions, which does not come free of drawbacks and limitations, potentially affecting their tomographic abilities. PR algorithms solve an inverse problem (i.e. recovery of the phase [16]) relying on the estimation of the autocorrelation (or, equivalently the Fourier modulus) of the object. The fact that autocorrelations are invariant to translation and to the transpose operator, imply that also the solutions are unique in a broader sense of the term [16,20]. Every solution obtained from a random initialization of the PR might differ for spatial positioning and diagonal flip, due to the mathematical properties of the autocorrelation. This is not a major problem, since the loss of spatial information does not interfere with the reconstruction performance. On the other hand, we have considered a fully transmitting sample that is identical under clockwise or counter-clockwise rotations, an assumption that might not be valid in experimental measurements. In this case, the sense of rotation matters and iPRT would be requested to start from a nontransposed version of the object in order to appropriately converge into meaningful reconstructions. This is easy to implement in case of direct observations but might be more difficult in the case of hidden measurements, thus requesting further studies specifically tuned to target this problematic. Another, but not less important, weakness of PR is the convergence rate in case of noisy measurements that might push the iteration towards a not optimal reconstruction. For this purpose, novel PR implementation (such as the Oversampling Smoothness [21] or Sparsity Phase Retrieval [22]) are promising choices to push the retrieval process towards robust reconstructions. In any case, we are planning to proceed ahead with the study of PR algorithms for tomographic implementations, aiming at an advanced understanding of other possible problematics and proposing new ways to tackle them.
5. Conclusions During this study, we described how optical projection tomography reconstructions could be aided by the usage of phase retrieval-based techniques [16]. With the specific purpose of tackling the limitations of OPT related to alignment and performing tomographic imaging of hidden features, the PR protocols presented in this work always allowed three-dimensional tomographic reconstructions. Interestingly, with both PRT and iPRT, we did not care for the center of rotation of the specimen [3] and neither for eventual sample movements [4] during the measurements, since the autocorrelations are insensitive to the absolute spatial positioning of the object. The methods proposed leave the user free of tedious post-alignment procedures [6], proposing an innovative solution to prevent artifacts generated by classical reconstruction schemes [15]. Moreover, the PR-based tomographic approaches proposed are naturally suitable for hidden threedimensional reconstructions, extending already existing 2D implementations [11,12] up to the 3D tomographic level. In this article, we have discussed how the PRT returns more accurate and alignment-free reconstructions, while iPRT cannot overcome the limitations on the alignment process of the classic methods. The key point on the success of PRT is associated to the fact that the autocorrelation of the object can be always estimated better than the object itself, due to its inherent properties of being aligned in the center of the space. Instead, iPRT exhibits instabilities in the reconstruction of the sinograms that can be compared with the effect of drifts and vibrational perturbations present in classic measurement schemes. A question then might arise at this point, concerning whether it is worth to use the iPRT method
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
8
D. Ancora et al. / Methods xxx (2017) xxx–xxx
Fig. 6. Integrated plot profiles comparison. In column 1) the classical reconstruction approaches, column 2) phase retrieval approaches. Each row shows the profiles obtained by integrating the other two missing variables (first row shows position in x coordinate having integrated in dydz and so on).
rather than its more accurate three-dimensional counterpart. Among the fact that no alignment is needed in any case, iPRT shows two favorable characteristics compared to PRT. First, it can be extremely faster, relying on two-dimensional Fourier transformations rather than three-dimensional, as required by the PRT. Second, the higher computational effort has to be made in the retrieval of the first projection of the sinogram sequence, while the other could be obtained by using a smaller number of iterations of the PR algorithm. This is due to a certain degree of continuity in the phase difference between two adjacent projections and the usage of the previously retrieved guess is an excellent prior information. Moreover PRT is computational-demanding when initiating with a random guess for the phase, thus iPRT reconstruction could be useful to feed the PRT with an already good prior information [16]. In these terms, iPRT might alleviate the effort requested for the finalization of PRT reconstruction, possibly needing less iterations to converge to meaningful results. On the other hand, due to the above reasons, iPRT is more likely to be used on its own as a tool for hidden two-dimensional imaging of a dynamically changing object, but this exceeds the purpose of this study. In fact, the hidden sinogram reconstructed in Fig. 2 row D can be turned into a movie that represents a hidden object rotating approximately around its center of mass, thus offering the possibility to image real-time the dynamics of objects hidden behind scattering curtains. Although we have already used PRT for alignment-free reconstructions to image deep into biological tissues [7], we did not use it yet for experimental tomographic reconstruction of hidden objects, even if we are currently working toward this goal. Moreover, due to their nature, both iPRT and PRT techniques might be also extended as reconstruction schemes for any other projection based reconstruction, such as bright field OPT, X-ray computed tomography (XCT) or magnetic resonance imaging (MRI). We already foresee that, with minor flowchart modifications, the algorithms proposed could aid other tomographic reconstruction techniques that rely on sample vertical translation such as helical OPT (hOPT) [23] or time-lapse OPT [24]. In our future studies we plan to approach these similar methodologies, with the explicit aim to further enhance the quality of tomographic reconstructions and, at the same time, leaving the user free of any calibration tasks. We believe that on this purpose, three-dimensional phase-retrieval reconstruction methods will play an important role further
enhanced by the possibility of massive parallelization of the algorithms, easily implementable with modern GPU computing. The limitations of the algorithms, at the moment, are merely technical and substantially related to memory handling of the volume reconstruction, but we are already working towards a more efficient parallel implementation of the phase retrieval iterations. Lastly, both methods could be potentially improved and tested under different experimental conditions, such as in the presence of noisier measurements. We did not examine yet this experimental regime where it is known that the phase retrieval encounters difficulties to find the appropriate convergence [16] while, among others, some new approaches have been proposed to tackle these problems. An efficient solution seems to be offered by the recent oversampling smoothness (OSS) protocol [21], that we are currently including into our protocol, together with ways to invert the autocorrelation sinogram. In fact, other techniques, such as simultaneous iterative reconstruction techniques (SIRT) and simultaneous algebraic reconstruction techniques (SART) on the GPU [25], are interesting alternatives to classical filtered backprojection reconstruction and are worth to be tested for the reconstruction of the three-dimensional autocorrelation. Nevertheless, phase retrieval algorithms could highly impact the field of computed tomography, offering increased flexibility under various experimental conditions and justifying our efforts in the development of novel PRbased reconstruction methodologies. Acknowledgements This work was supported by the EU Marie Curie ITN ‘‘OILTEBIA” PITN-GA-2012-317526 and the H2020 Laserlab Europe (EC-GA 654148). J. R. acknowledges support from the EC-FP7-CIG grant HIGH-THROUGHPUT TOMO and MINECO grant FIS2016-77892-R. The Authors also acknowledge support from the Institute of Molecular Biology and Biotechnology of the Foundation for Research and Technology – Hellas for the use of the cell culture facilities. References [1] V. Ntziachristos, Going deeper than microscopy, Nat. Methods 7 (2010) 603– 614. [2] J. Sharpe, U. Ahlgren, P. Perry, B. Hill, A. Ross, J. Hecksher-Sorensen, R. Baldock, D. Davidson, Optical projection tomography as a tool for 3D microscopy and gene expression studies, Science 296 (5567) (2002) 541–545.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
D. Ancora et al. / Methods xxx (2017) xxx–xxx [3] D. Dong, S. Zhu, C. Qin, V. Kumar, J.V. Stein, S. Oehler, C. Savakis, J. Tian, J. Ripoll, Automated recovery of the center of rotation in optical projection tomography in the presence of scattering, IEEE J. Biomed. Health. Inform. 17 (1) (2013) 198– 204. [4] U.J. Birk, M. Rieckher, N. Konstantinides, A. Darrell, A. Sarasa-Renedo, H. Meyer, N. Tavernarakis, J. Ripoll, Correction for specimen movement and rotation errors for in-vivo optical projection tomography, Biomed. Opt. Express 1 (1) (2010) 87–96. [5] M. Rieckher, U.J. Birk, H. Meyer, J. Ripoll, N. Tavernarakis, Microscopic optical projection tomography in vivo, PLOS one 6 (4) (2011) e18963. [6] A. Cheddad, C. Svensson, J. Sharpe, F. Georgsson, U. Ahlgren, Image processing assisted algorithms for optical projection tomography, IEEE Trans. Med. Imaging 31 (1) (2012) 1–15. [7] D. Ancora, D. Di Battista, G. Giasafaki, S.E. Psycharakis, E. Liapis, G. Zacharakis, Phase-retrieved tomography enables imaging of a tumor spheroid in mesoscopy regime, Sci. Rep. 7 (2017). [8] P.J. Verveer, J. Swoger, F. Pampaloni, K. Greger, M. Marcello, E.H.K. Stelzer, High-resolution three-dimensional imaging of large specimens with light sheet–based microscopy, Nat. Methods 4 (4) (2007) 311–313. [9] J. Huisken, J. Swoger, F. Del Bene, J. Wittbrodt, E.H.K. Stelzer, Optical sectioning deep inside live embryos by selective plane illumination microscopy, Science 305 (5686) (2004) 1007–1009. [10] J. Huisken, D.Y.R. Stainier, Selective plane illumination microscopy techniques in developmental biology, Development 136 (12) (2009) 1963–1975. [11] J. Bertolotti, E.G. van Putten, C. Blum, A. Lagendijk, W.L. Vos, A.P. Mosk, Noninvasive imaging through opaque scattering layers, Nature 491 (7423) (2012) 232–234. [12] O. Katz, P. Heidmann, M. Fink, S. Gigan, Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations, Nat. Photon. 8 (2014) 784–790. [13] S. Feng, C. Kane, P.A. Lee, D.A. Stone, Correlations and fluctuations of coherent wave transmission through disordered media, Phys. Rev. Lett. 61 (7) (1988) 834–837.
9
[14] D. Ancora, D. Di Battista, G. Giasafaki, S. Psycharakis, E. Liapis, A. Zacharopoulos, G. Zacharakis, Optical projection tomography via phase retrieval algorithms for hidden three dimensional imaging, Proc. SPIE 10074 (2017) 100741E. [15] J.R. Walls, J.G. Sled, J. Sharpe, R.M. Henkelman, Correction of artefacts in optical projection tomography, Phys. Med. Biol. 50 (19) (2005) 4645. [16] Y. Shechtman, Y.C. Eldar, O. Cohen, H.N. Chapman, J. Miao, M. Segev, Phase retrieval with application to optical imaging, IEEE Signal Process. Mag. 32 (3) (2015) 87–109. [17] R.W. Gerchberg, W.O. Saxton, A practical algorithm for the determination of the phase from image and diffraction plane pictures, Optik 35 (1972) 237. [18] J.R. Fienup, Phase retrieval algorithms: a comparison, App. Opt. 21 (15) (1982) 2758–2769. [19] J.R. Fienup, Phase retrieval algorithms: a personal tour [Invited], Appl. Opt. 52 (1) (2013) 45–56. [20] T. Bendory, R. Beinert and Y. C. Eldar, Fourier Phase Retrieval: Uniqueness and Algorithms, arXiv, p. arXiv:1705.09590, 2017. [21] J.A. Rodriguez, R. Xu, C.-C. Chen, Y. Zou, J. Miao, Oversampling smoothness: an effective algorithm for phase retrieval of noisy diffraction intensities, J. Appl. Crystallogr. 46 (2) (2013) 312–318. [22] K. Jaganathan, S. Oymak, B. Hassibi, Sparse Phase Retrieval: Uniqueness Guarantees and Recovery Algorithms, IEEE Trans. Signal Process. 65 (9) (2017) 2402–2410. [23] A. Arranz, D. Dong, S. Zhu, M. Rudin, C. Tsatsanis, J. Tian, J. Ripoll, Helical optical projection tomography, Optics Express 21 (22) (2013) 25912–25925. [24] A. Arranz, D. Dong, S. Zhu, C. Savakis, J. Tian, J. Ripoll, In-vivo optical tomography of small scattering specimens: time-lapse 3D imaging of the head eversion process in Drosophila melanogaster, Sci. Rep. 4 (2014) 7325. [25] D. Castaño-Díez, D. Moser, A. Schoenegger, S. Pruggnaller, A.S. Frangakis, Performance evaluation of image processing algorithms on the GPU, J. Struct. Biol. 164 (2008) 153–160.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
Contents lists available at ScienceDirect
Methods journal homepage: www.elsevier.com/locate/ymeth
Optical projection tomography via phase retrieval algorithms Daniele Ancora a,b, Diego Di Battista a,c, Georgia Giasafaki a, Stylianos E. Psycharakis a,d, Evangelos Liapis a, Jorge Ripoll e,f, Giannis Zacharakis a,⇑ a
Institute of Electronic Structure and Laser, Foundation for Research and Technology – Hellas, 70013 Heraklion, Greece Department of Materials Science and Technology, University of Crete, 71003 Heraklion, Greece c Assing S.p.A, Monterotondo, 00015 Rome, Italy d School of Medicine, University of Crete, 71003 Heraklion, Greece e Department of Bioengineering and Aerospace Engineering, Universidad Carlos III de Madrid, 28911 Madrid, Spain f Instituto de Investigación Sanitaria del Hospital Gregorio Marañón, 28007 Madrid, Spain b
a r t i c l e
i n f o
Article history: Received 7 August 2017 Received in revised form 15 October 2017 Accepted 17 October 2017 Available online xxxx Keywords: Optical projection tomography Phase retrieval Registration-free imaging Autocorrelation imaging Hidden object tomography
a b s t r a c t We describe a computational method for accurate, quantitative tomographic reconstructions in Optical Projection Tomography, based on phase retrieval algorithms. Our method overcomes limitations imposed by light scattering in opaque tissue samples under the memory effect regime, as well as reduces artifacts due to mechanical movements, misalignments or vibrations. We make use of Gerchberg–Saxton algorithms, calculating first the autocorrelation of the object and then retrieving the associated phase under four numerically simulated measurement conditions. By approaching the task in such a way, we avoid the projection alignment procedure, exploiting the fact that the autocorrelation sinogram is always aligned and centered. We thus propose two new, projection-based, tomographic imaging flowcharts that allow registration-free imaging of opaque biological specimens and unlock three-dimensional tomographic imaging of hidden objects. Two main reconstruction approaches are discussed in the text, focusing on their efficiency in the tomographic retrieval and discussing their applicability under four different numerical experiments. Ó 2017 Elsevier Inc. All rights reserved.
1. Introduction Optical biomedical imaging modalities have become essential tools for disease detection and monitoring, occupying an increasing portion of medical imaging needs [1], characterized by the use of non-ionizing radiation, molecular specific contrast agents and minimum or non-invasiveness. Among these methods, Optical Projection Tomography (OPT) [2], considered the optical analogue of X-ray Computed Tomography (XCT), shares with the latter a broad set of reconstruction approaches, while the measurement scheme might change slightly depending on the experimental architecture. OPT relies on mechanical movement of the specimen rotating it around a fixed axis to acquire the full set of projections, rather than moving the measuring device around the object. However, mechanical translation and rotations are then needed to perform a three-dimensional reconstruction of the sample of interest, thus requiring a correct calibration of the experiment to obtain accurate reconstructions. Such calibration has to be often performed both at hardware and software level and might sometimes limit the acqui⇑ Corresponding author. E-mail address: [email protected] (G. Zacharakis).
sition procedure or introduce unwanted artifacts in the reconstructions. Many efforts have been made in finding new algorithmic ways to align, center or reduce misalignment effects at the tomographic level [3–6], but none has the advantage of being fully alignment-free. In this article, we study a more generic approach to confront with misaligned data reconstructions, proposing two algorithmic step-wise protocols that are able to retrieve quantitative reconstructions even in the presence of constant drifts, random vibrations or even overcoming the phase scrambling due to the presence of a scattering layer enclosing the object. The study focuses on two different approaches, both of them relying on the use of phase retrieval algorithms, to generate aligned autocorrelation sinograms that can be inverted in order to obtain tomographic reconstructions. We emphasize the fact that any positional information of the specimen is not important, since the autocorrelation is inherently always centered in its space of definition. Phase-Retrieved Tomography (PRT) was originally introduced in our previous work [7] to image experimentally the fluorescence distribution of necrotic cells in a human-breast tumor spheroid [8]. In this case the method was compared against classical reconstruction schemes from Single Plane Illumination Microscopy (SPIM) [9,10] and OPT, showing promising reconstruction abilities. Here
https://doi.org/10.1016/j.ymeth.2017.10.009 1046-2023/Ó 2017 Elsevier Inc. All rights reserved.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
2
D. Ancora et al. / Methods xxx (2017) xxx–xxx
we focus on a broader numerical study, better illustrating the potential of phase retrieval techniques applied for tomographic investigations. Alongside PRT, we propose another iterative technique to recover aligned sinograms to which we refer to as iterative Projection-Retrieved Tomography (iPRT). The two methods will be tested for the first time in reconstructing corrupted datasets generated by numerical simulations, focusing in particular on their ability to leave the user free of any alignment procedure. The possibility of performing hidden reconstructions will be further considered showing that, under memory effect conditions, blind tomography becomes naturally possible as a result of the incorporation of phase retrieval techniques into classical reconstruction schemes. 2. Material and methods In this article we generate the datasets numerically, analyzing the possibility offered by the methods under different measurement scenarios in order to have a known reference to compare with. The object under investigation is a previously retrieved reconstruction of the fluorescence distribution in a human-breast tumor spheroid [7], which we will use as a ground truth measurement. This volume has a size of 300 pixels (from now on abbreviated as px) in each spatial dimension, being characterized in the original experiment by a pixel size of 0:8 lm. Starting with this volume, we generate several datasets representing four different experimental measurement conditions. The first dataset is created by rotating the object around an axis parallel to z and positioned at the center of the camera plane, thus acquiring its intensity projections (Fig. 1, row A) as a function of the observation angle. This mathematically corresponds to the Radon transform of the sample volume and, experimentally, to an ideal OPT experiment. The second dataset (row B) is obtained by introducing a random twodimensional slight vibrations of 1 px during the rotation, leading
to the generation of a fuzzy sinogram. This could represent a situation in which we try to correct the axis of rotation to fit it onto the position of the barycenter of the object, or the case in which the stages slightly vibrate during the rotation. The third dataset (row C) is equal to the above except for the introduction of a constant drift in two dimensions while rotating the object. In fact, in some experimental acquisitions [7], we observed a similar perturbation due to the effect of gravity (sample slowly falling inside the medium) and due to mechanical misalignment of the rotating stage. The fourth, and last, dataset (row D) reproduces a measurement performed when the object is hidden behind a phase scrambling environment [11,12], such as a scattering envelope (e.g. an egg or a cocoon). This dataset is numerically calculated assuming the measurements to be in the memory effect regime [13], which states that a speckle pattern (generated by a single point source) translates in relation to the source movement behind the phasescrambling layer [14]. All the above described numerical experiments produced different sinograms (Fig. 1, column 2), that can be inverted (with an inverse Radon transformation) to obtain the volume reconstruction in each specific case. Column 3 of Fig. 1 shows, depending on the perturbation introduced in the measurement, how the reconstructed volume differs from the ideal case ultimately producing several kinds of unwanted reconstruction artifacts well known in the literature [4,15]. In the following, we will try to avoid the generation of alignment perturbations, exploiting the mathematical properties of the autocorrelation in combination with the phaseretrieval algorithm. 2.1. Autocorrelation sinogram From each of the previously obtained sinograms we can calculate the corresponding autocorrelation sinogram Ah (A-sinogram), which will be used for both the phase-retrieval methods. In this
Fig. 1. Four different numerical experiments, showing the resulting sinograms (column 2) of the three different rotations. A) The results of an aligned rotation, leading to optimal reconstructions of the object. B) A one-pixel random vibration was introduced while rotating the object, returning slight artifacts in the reconstruction. C) The measurement was affected by a two-dimensional drift (gravity falling and axis displacement) particularly evident in the corresponding reconstruction. D) Hidden object measurement, impossible to backproject with classical techniques.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
D. Ancora et al. / Methods xxx (2017) xxx–xxx
notation, h is the angle at which we acquire the projections that spans through the entire 360° rotation in steps of 1°. The calculation of the A-sinogram is simple for the direct projections and can be estimated by using the Wiener-Kinchin theorem:
Ah ¼ F 1 fjF fPh gj2 g:
ð1Þ
Instead, to obtain the A-sinogram in the case of hidden measurements we need to average the autocorrelation over different speckle-sinograms obtained by using different random phases, in analogy to the ensemble averages requested by the corresponding two-dimensional technique [11]. In our specific case, we averaged over 10 different random phases, but in general a higher number of phases corresponds to a better sampling of the object autocorrelation. In Fig. 2 column 1 (row A), we display the result of the calculation of the A-sinogram in the case of direct measurements (described in Fig. 1, rows A-C), which results to be identical in the three cases examined and then insensitive to any drift or vibration of the sample during the rotation. Fig. 2 column 1 (row C) shows the A-sinogram calculated from the hidden measurement (Fig. 1, row D), which exhibits the same features of the one calculated with direct measurements. Moreover, the A-sinogram resulting from the hidden experiment gets more accurate when we increase the ensemble average, and the noise (dark background of the A-sinogram in column 1, row C) can be further suppressed with a more accurate sampling. 2.2. Iterative projection retrieval tomography (iPRT) A straightforward approach, once we have calculated the Asinogram, is to proceed stepwise aiming at the retrieval of every single projection, thus rebuilding an aligned sinogram. To do so, we firstly retrieve the phase connected with the 2D autocorrelation of the initial projection, e.g. at angle h ¼ 0, with a classical phase
3
retrieval (PR) algorithm [16] starting with a random-phase guess. The position of the object retrieved in this way will be located randomly in space, due to the fact that the autocorrelation operation is invariant for any object translation. Running independent PR per each projection would return an unpredictable object shift in the camera plane, rendering a full reconstruction impossible. However, if we assume that the object is rotating in small steps, we could use the obtained phase at the starting angle /ðhÞ to guide the PR reconstruction of the next projection at angle h þ 1. In this way, the PR problem converges faster, since we already start with a guess that is very close to the phase solution /ðh þ 1Þ. Iterating this process until the last projection N and stacking every PR reconstruction with the same order of the A-sinogram, we end up having a seemingly aligned sinogram. Since the sinogram was build starting from a reconstruction located in a random spatial position, the rotation axis of such dataset is in general not centered in the camera plane, but it can be automatically found by comparing the reconstruction at 0° and 180°. Once centered, the obtained sinogram can be backprojected (e.g. with a filtered inverse Radon transformation) to retrieve the final reconstruction of the specimen. The whole process is described in the flowchart A of Fig. 3 and we refer to it as iterative Projection-Retrieved Tomography (iPRT). It is important to notice that we do not strictly have to insert the first projection into the iPRT in case we have direct access to the object image. In fact, by using the phase of the initial projection as the starting point of iPRT we will be able to align the sinogram with respect to the initial measurement, having no need to retrieve once again the phase of the initial step. Although this implies a further degree of simplicity (we do not have to run the first PR iterations), for the purpose of this study we will still force the retrieval of the phase of the first projection. We consider this situation to stress the fact that the overall phase is not a necessary aspect for realigning projections and applying iPRT directly to hidden measurements, where the initial object projection is not directly accessible.
Fig. 2. Autocorrelation sinograms (column 1) obtained for the four numerical experiments and their consequent iPRT sinogram reconstruction. A) PRT imaging for any of the measurement scheme (perfect rotation, vibration and drift), in column 2) it is shown the final three-dimensional object autocorrelation and in 3) the final reconstruction. B) iPRT volume reconstruction of the same datasets. C) Hidden PRT reconstruction, showing a nearly similar autocorrelation compared to that of A). D) Hidden iPRT reconstruction.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
4
D. Ancora et al. / Methods xxx (2017) xxx–xxx
Fig. 3. Block diagram of the phase retrieval methods for OPT reconstructions. A) The iPRT flowchart while B) The simpler PRT flowchart. Although iPRT has a more complicated flowchart, the method is less computational demanding, due to the fact that it relies on 2D Fourier transformations rather than 3D.
2.3. Phase-retrieved tomography (PRT) A different approach for the alignment of the dataset can be implemented by backprojecting directly the A-sinogram that was previously calculated. Due to the fact that the A-sinogram is perfectly centered, the volume reconstructed in this way will not suffer any misalignment, while it can be proven that is mathematically equivalent to the three-dimensional autocorrelation of the whole specimen [7]. The phase retrieval problem, then, can be extended to three dimensions by starting with a random 3Dphase, thus its solution would directly give the reconstruction of the object of interest with no further need of alignment procedures. Hence, we refer to this method as Phase-Retrieved Tomography (PRT) and its corresponding flowchart is shown in panel B of Fig. 3. We point out that although the flowchart is simpler than that of iPRT, the PRT method is computationally more demanding, since the phase retrieval algorithm employed for the solution is based on three-dimensional Fourier transforms, rather than twodimensional. 2.4. Phase retrieval algorithms Both of the methods proposed, base their actions on the usage of Gerchberg–Saxton phase retrieval (PR) algorithms [17] in order to retrieve the phase connected with the autocorrelation, or more appropriately, with the estimation of the modulus of the Fourier transform jUðkÞj. The retrieval of the correct phase, then, enables the reconstruction of either the single projections (iPRT) or the whole object (PRT). For both of the methods, we used a common combination of the hybrid input-output (HIO) and the error reduction (ER) [18] algorithms, starting with 2000 iteration of the HIO followed by 100 iterations of the ER. Both algorithms fall into the class of alternating projection methods [19,16], in which we pass
from the real (or object) space to the Fourier space applying appropriate constraining operations. The two approaches follow a simple sequence of operations composed by four steps, in which only the last one changes according to the algorithm used. Without losing in generality, in the following we will treat the one-dimensional case that can be trivially extended to two- or three-dimensions, respectively implemented into iPRT and PRT reconstruction protocols. As a first step, the PR algorithm starts with an initial random guess of the object g 0 ðxÞ then, at the j-th iteration, the algorithm acts with this sequence of operation: 1. Fourier transform the object estimation:
g j ðxÞ ! Gj ðkÞ ¼ F fg j ðxÞg
ð2Þ
2. Substitute its modulus with the originally measured one:
Gj ðkÞ ¼ jGj ðkÞjei/ðkÞ ! G0j ðkÞ ¼ jUðkÞjei/ðkÞ
ð3Þ
3. Calculate the inverse Fourier transform of this estimation: 0ðkÞ
Gj
! g 0j ðxÞ ¼ F 1 fG0j ðkÞg
ð4Þ
4. For the HIO, calculate the object update by using a feedback coming from the previous iterations in C, the region where the assumption of the image being real and positive are violated: (
for x R C
g 0j ðxÞ
g jþ1 ðxÞ ¼
g j ðxÞ bg 0j ðxÞ for x 2 C
ð5Þ
Or, for the final ER iterations, setting to zero the values in the region C:
(
g jþ1 ðxÞ ¼
g 0j ðxÞ for x R C 0
for x 2 C
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
ð6Þ
5
D. Ancora et al. / Methods xxx (2017) xxx–xxx
Ideally, such process is iterated until the modulus of the Fourier transform of the calculated object jGj ðkÞj converges to the measured one jUðkÞj, or when there are no points left within the region C. In practice, the convergence of the algorithm can be monitored calculating the Recovery Error Function (REF), defined as the Euclidean distance between the modulus of the Fourier transform and the measured one:
REF ¼
X jjGj ðkÞj jUðkÞjj2 :
ð7Þ
k
This quantity converges to zero, REF ! 0, when the solution approaches the exact one. Only for the case of iPRT at the second step of the sinogram reconstruction (h ¼ 1) we use an estimation of the object given by the previous reconstruction (instead of the random guess), and we consequently reduce the number of the steps of the HIO to 100. From now on, we will refer to PR algorithms as the two or three-dimensional processes that retrieve any phase from the corresponding autocorrelations; we will refer to PRT and iPRT as the overall imaging process, involving PR reconstructions at different steps described by the flowcharts of Fig. 3.
Considering only the translational variable g, we notice that it is possible to change the variable of integration to:
Z
Oðx þ g; y þ n; z þ vÞdg ¼
Z
Oðx þ g; y þ n; z þ vÞdx
¼ PO ðy þ n; z þ vÞ:
ð15Þ
In fact, the projection on g is independent on any x shift parallel to its direction. This ultimately implies that we can write:
PO ðy; zÞPO ðy þ n; z þ vÞ ¼
Z
Oðx; y; zÞPO ðy þ n; z þ vÞdx
ð16Þ
where it is now obvious that both the terms are equal, thus proving the initial assumption in Eq. (4). Since this operation can be repeated at each angle of observation, the backprojection of the A-sinogram leads to the reconstruction of the 3D autocorrelation of the specimen. We conclude that the association of the proper phase to this quantity, thus solving the three-dimensional phase retrieval problem, would enable us to retrieve the final reconstruction of the object.
3. Theory and calculation
4. Results and discussion
3.1. PRT theoretical background
In this section, we compare the results obtained for iPRT and PRT reconstruction schemes. We point out that in all retrieval processes studied in the present work, we always obtained consistent reconstructions starting from random guesses, differing only from each other by a random spatial translation. This is due to the fact that the autocorrelation is degenerate for any spatial translation of the sample, which consequently leads to an overall loss of information regarding the whole object positioning in space. Each reconstructed volume, then, was automatically translated back to the position of the ground truth volume by finding the absolute maximum of their cross-correlation. This procedure was necessary in order to compare the results between the different measuring schemes, but it is not required for the reconstruction itself. Fig. 2 shows the entire reconstruction process in a progressive diagram. Column 1 displays the A-sinograms (rows A and C) and the corresponding realigned sinograms calculated from those (row B and D). Column 2 visualizes the rendering of the three-dimensional autocorrelations calculated from the A-sinograms (only for rows A and C). Column 3 shows the final rendering of all the reconstructed volumes for iPRT and PRT. In more detail, Fig. 2 row A describes the PRT method applied to reconstruct misaligned datasets. The A-sinogram of column 1 was calculated from any of the direct sinogram observations of Fig. 1 column 2 (raw A to C), resulting perfectly aligned and identical for all the cases. Such A-sinogram is backprojected (hence, inverse Radon transformed) leading to the calculation of the autocorrelation volume of the specimen of interest, rendered in column 2 for visualization purposes. Subsequently, a three-dimensional phase retrieval algorithm was used to find the phase associated to this volume, enabling the imaging of the object. This finalizes the PRT reconstruction protocol of which the final result is displayed in column 3. In this case, the retrieved fluorophore distribution appears to be accurate and insensitive to any of the perturbations described in Fig. 1 (random vibrations and drifts), seemingly equivalent to the original ground truth (Fig. 1 column 1, row A). On the other hand, illustrated in Fig. 2 row B, the iPRT managed to reconstruct the sinogram from the autocorrelation sequence previously described, without relying on any spatial information of the object rotation. Although the retrieved sinogram (displayed in column 1, row B) seems to rotate smoothly, it contains a certain degree of misalignment and results to the production of artifacts when we backproject it for the final reconstruction. The result of the iPRT protocol is rendered in col-
Regardless of how the A-sinogram was calculated, the PRT method relies on the fact that its backprojection leads to the reconstruction of the three-dimensional autocorrelation of the sample. This can be mathematically verified since the problem can be reduced to the fact that the projection of the 3D autocorrelation must be equal to the 2D autocorrelation of the object projection per every angle of the measurement. Let us define the object of interest as Oðx; y; zÞ and its 3D autocorrelation Aðg; n; vÞ as:
Aðg; n; vÞ ¼ Oðx; y; zÞHOðx; y; zÞ Z ¼ Oðx; y; zÞOðx þ g; y þ n; z þ vÞdx dy dz:
ð8Þ
We define the observation angle h ¼ 0 to be parallel to the x axis (and consequently also to g), thus the projection of the object at this angle results equal to:
Z
PO ðy; zÞ ¼ Px ¼
ð9Þ
Oðx; y; zÞdx
and the projection of the autocorrelation:
PA ðn; vÞ ¼
Z
Aðg; n; vÞdg:
ð10Þ
Due to the above consideration, the problem is reduced to verify that the following relation is satisfied:
PO ðy; zÞHPO ðy; zÞ ¼ PA ðn; vÞ: Z Z
ð11Þ
We can write the terms in Eq. (9) as:
PO ðy; zÞPO ðy þ n; z þ vÞdy dz ¼ PO ðy; zÞPO ðy þ n; z þ vÞdydz ¼
Z
Z
Aðg; n; vÞdg
ð12Þ
Oðx; y; zÞOðx þ g; y þ n; z
þ vÞdx dy dz dg:
ð13Þ
Since in both the integrals we are integrating along y and z, we can then compare the two arguments, leading to the simplified formulation:
PO ðy;zÞPO ðy þ n;z þ vÞ ¼
Z
Oðx; y; zÞOðx þ g; y þ n; z þ vÞdxdg:
ð14Þ
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
6
D. Ancora et al. / Methods xxx (2017) xxx–xxx
umn 3. It is clear, in this case, that PRT performed better than iPRT retrieving a sharper reconstruction, although also iPRT converged into a meaningful (but noisier) result. Fig. 2 rows C and D show respectively the PRT and iPRT for the reconstructions of the hidden object, following exactly the same scheme described for the realignment of direct measurements. Both reconstructions were applied to the dataset of Fig. 1 row D, calculated assuming that the object was secreted by a scattering layer within the memory effect regime. In normal imaging modalities, the acquisition of a seemingly-random speckle pattern did not allow any classical reconstruction methods (volume in Fig. 1 column 3, row D), thus rendering the use of the phase retrieval protocols essential for any tomographic analysis. The A-sinogram (Fig. 2 column 1, raw C) results perfectly aligned and identical to the one obtained in the case of direct measurements, being characterized by a noisier background. The autocorrelations were calculated averaging 10 different random speckle patters, but we verified that a more accurate sampling would push further down the noise in the A-sinogram. In any case, its inversion returned a very good reconstruction of the three-dimensional autocorrelation of the sample (column 2, raw C), visually identical to the one reconstructed from direct measurements. To finalize the imaging process, the phase retrieval finds the phase to connect with such volume, thus leading to a very accurate reconstruction of the sample concealed by a scattering layer (column 3, raw C). From the Asinogram of the hidden measurements, we could generate back the direct sinogram of the object rotation by running the iPRT protocol (column 1, raw D). In this case, the resulting sinogram is very similar to the one obtained by direct measurements and leads to a consistent reconstruction of the specimen (column 3, raw D). In close analogy with the previous case, PRT returned better reconstructions also in the case of blind tomography if compared with iPRT, leading to sharp imaging of the fluorescence distribution. The resulting volume has a noisier background if compared to the PRT for the normal dataset, but still represents an excellent reconstruction for a hidden measurement. Also the iPRT converged into a correct solution, but this time leaving more artifacts if compared to the previous case. To have a more comprehensive understanding about the quality of the reconstruction, it is interesting to compare the projections along the three-spatial dimensions. Thus, calling any final recon-
structed volume as Vðx; y; zÞ, peak-normalized to unity, we calculate:
Z Px ¼
Z Vðx; y; zÞdx; Py ¼
Z Vðx; y; zÞdy; Pz ¼
Vðx; y; zÞdz;
ð17Þ
which lead to the production of three projected images per each reconstruction approach (Fig. 4). We can appreciate that PRT definitely solved the problem of misalignment (panel E), returning a nearly perfect reconstruction that matches that of a perfectly aligned measurement (shown in panel A). Instead, iPRT (panel F) used as a tool to realign the projection sinogram, returned a reconstruction having artifacts similar to the reconstruction perturbed by a constant drift (panel C). On the other hand, also, hidden PRT retrieved a perfectly aligned fluorophore distribution, but the reconstruction lies in a noisier environment if compared to that of the normal PRT (panel G). Lastly, the volume produced by the hidden iPRT reconstruction (panel H) seems to be affected by the same vibrational artifacts typical of a reconstruction in a vibrationalperturbed case (panel B). Fig. 5 confirms this fact, presenting the analysis of the residual misalignment of the reconstruction compared to the original, thus aligned, dataset. Here we use the aligned rotation as reference to measure, as a function of the angle, the position of the perturbed measurements and the reconstructed sinograms obtained with iPRT. Panel A presents the behavior of the displacement with respect to the horizontal axis (X in the sinogram): the residual of the iPRT reconstruction (violet curve) drifts away with respect to the aligned dataset (blue) more similarly to the drifted dataset (yellow) than to the vibrationally perturbed one. The reconstruction of hidden object with iPRT (green), instead, suffers a slight drift but stronger vibrational misalignment, making it more similar to the vibrationally perturbed sinogram (red). The residual alignment of the vertical axis Z (parallel to the rotation axis) shown in panel B presents another interesting situation, in which iPRT reconstruction does not suffer strong drift-effect but still a vibrational error affects the hidden iPRT reconstruction. In synthesis, iPRT reconstruction drifts gently away with respect to the original dataset in the horizontal direction, while there is no appreciable effect on the vertical axis. On the other hand, hidden iPRT clearly suffers vibrational perturbation on both axis thus leading to the visual artifacts discussed in Fig. 4.
Fig. 4. Reconstruction comparison, projection view. A) reconstruction from perfectly aligned measurement, B) vibrationally perturbed rotation, D) two-dimensional constant drift and E) hidden measurements. E-F) Corresponding reconstruction obtained respectively with PRT and iPRT by using the same dataset used for A-B). G-H) reconstructions obtained from hidden dataset.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
D. Ancora et al. / Methods xxx (2017) xxx–xxx
Fig. 5. Residual shifts of the perturbed and iPRT reconstructed sinograms with respect to the aligned datasets. A) Displacement with respect to the horizontal direction X, perpendicular to the rotation axis. B) Displacement with respect to the vertical direction Z, parallel to the axis of rotation.
Qualitatively, we observe that every reconstruction was able to retrieve a correct fluorophore distribution, but in general PRT works far better than iPRT in terms of sharpness and accuracy of the results. In any case, both the techniques overtook the limitations of performing hidden reconstruction (Fig. 4, panel D), opening novel promising imaging scenarios. To quantitatively compare the results, we focus on plot profiles calculated from projections onto two different axes, according to the following calculation:
Z Pxy ¼
Z Vðx; y; zÞdxdy; P yz ¼
Z Vðx; y; zÞdydz; Pzx ¼
Vðx; y; zÞdxdz; ð18Þ
From the analysis of the profiles in Fig. 6, we conclude that phase retrieval methods quantitatively retrieve the intensity distribution and that the results of PRT are as accurate as the ones from an aligned reconstruction. We can also observe the drift-effect introduced by iPRT reconstruction (yellow line, Fig. 6 panel D), similar to the effect of the reconstruction in the presence of drift (violet line, panel A, B). Of particular interest are also the results obtained for hidden PRT (violet line, panel D-F), in which a noisier background lifts up the projection profiles by a constant amount, still leaving the fluorophore distribution highly visible compared to the noise level. Interestingly then, we can conclude the analysis claiming that the PRT method performs better than iPRT for any measuring conditions and overcomes the limitations associated to the alignment of the sample movement. The latter, instead, suffer from two different alignment issues, depending on the kind of measurement initiating our analysis. Iterative reconstruction of subsequent single projections leads to a reconstruction of a sinogram affected
7
by a slight drift, while on the other hand, in the case of speckle measurement, it suffers artifacts typical of the vibrational noise. Finally, it is important to notice that all the methodologies proposed rely on PR algorithms to perform the reconstructions, which does not come free of drawbacks and limitations, potentially affecting their tomographic abilities. PR algorithms solve an inverse problem (i.e. recovery of the phase [16]) relying on the estimation of the autocorrelation (or, equivalently the Fourier modulus) of the object. The fact that autocorrelations are invariant to translation and to the transpose operator, imply that also the solutions are unique in a broader sense of the term [16,20]. Every solution obtained from a random initialization of the PR might differ for spatial positioning and diagonal flip, due to the mathematical properties of the autocorrelation. This is not a major problem, since the loss of spatial information does not interfere with the reconstruction performance. On the other hand, we have considered a fully transmitting sample that is identical under clockwise or counter-clockwise rotations, an assumption that might not be valid in experimental measurements. In this case, the sense of rotation matters and iPRT would be requested to start from a nontransposed version of the object in order to appropriately converge into meaningful reconstructions. This is easy to implement in case of direct observations but might be more difficult in the case of hidden measurements, thus requesting further studies specifically tuned to target this problematic. Another, but not less important, weakness of PR is the convergence rate in case of noisy measurements that might push the iteration towards a not optimal reconstruction. For this purpose, novel PR implementation (such as the Oversampling Smoothness [21] or Sparsity Phase Retrieval [22]) are promising choices to push the retrieval process towards robust reconstructions. In any case, we are planning to proceed ahead with the study of PR algorithms for tomographic implementations, aiming at an advanced understanding of other possible problematics and proposing new ways to tackle them.
5. Conclusions During this study, we described how optical projection tomography reconstructions could be aided by the usage of phase retrieval-based techniques [16]. With the specific purpose of tackling the limitations of OPT related to alignment and performing tomographic imaging of hidden features, the PR protocols presented in this work always allowed three-dimensional tomographic reconstructions. Interestingly, with both PRT and iPRT, we did not care for the center of rotation of the specimen [3] and neither for eventual sample movements [4] during the measurements, since the autocorrelations are insensitive to the absolute spatial positioning of the object. The methods proposed leave the user free of tedious post-alignment procedures [6], proposing an innovative solution to prevent artifacts generated by classical reconstruction schemes [15]. Moreover, the PR-based tomographic approaches proposed are naturally suitable for hidden threedimensional reconstructions, extending already existing 2D implementations [11,12] up to the 3D tomographic level. In this article, we have discussed how the PRT returns more accurate and alignment-free reconstructions, while iPRT cannot overcome the limitations on the alignment process of the classic methods. The key point on the success of PRT is associated to the fact that the autocorrelation of the object can be always estimated better than the object itself, due to its inherent properties of being aligned in the center of the space. Instead, iPRT exhibits instabilities in the reconstruction of the sinograms that can be compared with the effect of drifts and vibrational perturbations present in classic measurement schemes. A question then might arise at this point, concerning whether it is worth to use the iPRT method
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
8
D. Ancora et al. / Methods xxx (2017) xxx–xxx
Fig. 6. Integrated plot profiles comparison. In column 1) the classical reconstruction approaches, column 2) phase retrieval approaches. Each row shows the profiles obtained by integrating the other two missing variables (first row shows position in x coordinate having integrated in dydz and so on).
rather than its more accurate three-dimensional counterpart. Among the fact that no alignment is needed in any case, iPRT shows two favorable characteristics compared to PRT. First, it can be extremely faster, relying on two-dimensional Fourier transformations rather than three-dimensional, as required by the PRT. Second, the higher computational effort has to be made in the retrieval of the first projection of the sinogram sequence, while the other could be obtained by using a smaller number of iterations of the PR algorithm. This is due to a certain degree of continuity in the phase difference between two adjacent projections and the usage of the previously retrieved guess is an excellent prior information. Moreover PRT is computational-demanding when initiating with a random guess for the phase, thus iPRT reconstruction could be useful to feed the PRT with an already good prior information [16]. In these terms, iPRT might alleviate the effort requested for the finalization of PRT reconstruction, possibly needing less iterations to converge to meaningful results. On the other hand, due to the above reasons, iPRT is more likely to be used on its own as a tool for hidden two-dimensional imaging of a dynamically changing object, but this exceeds the purpose of this study. In fact, the hidden sinogram reconstructed in Fig. 2 row D can be turned into a movie that represents a hidden object rotating approximately around its center of mass, thus offering the possibility to image real-time the dynamics of objects hidden behind scattering curtains. Although we have already used PRT for alignment-free reconstructions to image deep into biological tissues [7], we did not use it yet for experimental tomographic reconstruction of hidden objects, even if we are currently working toward this goal. Moreover, due to their nature, both iPRT and PRT techniques might be also extended as reconstruction schemes for any other projection based reconstruction, such as bright field OPT, X-ray computed tomography (XCT) or magnetic resonance imaging (MRI). We already foresee that, with minor flowchart modifications, the algorithms proposed could aid other tomographic reconstruction techniques that rely on sample vertical translation such as helical OPT (hOPT) [23] or time-lapse OPT [24]. In our future studies we plan to approach these similar methodologies, with the explicit aim to further enhance the quality of tomographic reconstructions and, at the same time, leaving the user free of any calibration tasks. We believe that on this purpose, three-dimensional phase-retrieval reconstruction methods will play an important role further
enhanced by the possibility of massive parallelization of the algorithms, easily implementable with modern GPU computing. The limitations of the algorithms, at the moment, are merely technical and substantially related to memory handling of the volume reconstruction, but we are already working towards a more efficient parallel implementation of the phase retrieval iterations. Lastly, both methods could be potentially improved and tested under different experimental conditions, such as in the presence of noisier measurements. We did not examine yet this experimental regime where it is known that the phase retrieval encounters difficulties to find the appropriate convergence [16] while, among others, some new approaches have been proposed to tackle these problems. An efficient solution seems to be offered by the recent oversampling smoothness (OSS) protocol [21], that we are currently including into our protocol, together with ways to invert the autocorrelation sinogram. In fact, other techniques, such as simultaneous iterative reconstruction techniques (SIRT) and simultaneous algebraic reconstruction techniques (SART) on the GPU [25], are interesting alternatives to classical filtered backprojection reconstruction and are worth to be tested for the reconstruction of the three-dimensional autocorrelation. Nevertheless, phase retrieval algorithms could highly impact the field of computed tomography, offering increased flexibility under various experimental conditions and justifying our efforts in the development of novel PRbased reconstruction methodologies. Acknowledgements This work was supported by the EU Marie Curie ITN ‘‘OILTEBIA” PITN-GA-2012-317526 and the H2020 Laserlab Europe (EC-GA 654148). J. R. acknowledges support from the EC-FP7-CIG grant HIGH-THROUGHPUT TOMO and MINECO grant FIS2016-77892-R. The Authors also acknowledge support from the Institute of Molecular Biology and Biotechnology of the Foundation for Research and Technology – Hellas for the use of the cell culture facilities. References [1] V. Ntziachristos, Going deeper than microscopy, Nat. Methods 7 (2010) 603– 614. [2] J. Sharpe, U. Ahlgren, P. Perry, B. Hill, A. Ross, J. Hecksher-Sorensen, R. Baldock, D. Davidson, Optical projection tomography as a tool for 3D microscopy and gene expression studies, Science 296 (5567) (2002) 541–545.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009
D. Ancora et al. / Methods xxx (2017) xxx–xxx [3] D. Dong, S. Zhu, C. Qin, V. Kumar, J.V. Stein, S. Oehler, C. Savakis, J. Tian, J. Ripoll, Automated recovery of the center of rotation in optical projection tomography in the presence of scattering, IEEE J. Biomed. Health. Inform. 17 (1) (2013) 198– 204. [4] U.J. Birk, M. Rieckher, N. Konstantinides, A. Darrell, A. Sarasa-Renedo, H. Meyer, N. Tavernarakis, J. Ripoll, Correction for specimen movement and rotation errors for in-vivo optical projection tomography, Biomed. Opt. Express 1 (1) (2010) 87–96. [5] M. Rieckher, U.J. Birk, H. Meyer, J. Ripoll, N. Tavernarakis, Microscopic optical projection tomography in vivo, PLOS one 6 (4) (2011) e18963. [6] A. Cheddad, C. Svensson, J. Sharpe, F. Georgsson, U. Ahlgren, Image processing assisted algorithms for optical projection tomography, IEEE Trans. Med. Imaging 31 (1) (2012) 1–15. [7] D. Ancora, D. Di Battista, G. Giasafaki, S.E. Psycharakis, E. Liapis, G. Zacharakis, Phase-retrieved tomography enables imaging of a tumor spheroid in mesoscopy regime, Sci. Rep. 7 (2017). [8] P.J. Verveer, J. Swoger, F. Pampaloni, K. Greger, M. Marcello, E.H.K. Stelzer, High-resolution three-dimensional imaging of large specimens with light sheet–based microscopy, Nat. Methods 4 (4) (2007) 311–313. [9] J. Huisken, J. Swoger, F. Del Bene, J. Wittbrodt, E.H.K. Stelzer, Optical sectioning deep inside live embryos by selective plane illumination microscopy, Science 305 (5686) (2004) 1007–1009. [10] J. Huisken, D.Y.R. Stainier, Selective plane illumination microscopy techniques in developmental biology, Development 136 (12) (2009) 1963–1975. [11] J. Bertolotti, E.G. van Putten, C. Blum, A. Lagendijk, W.L. Vos, A.P. Mosk, Noninvasive imaging through opaque scattering layers, Nature 491 (7423) (2012) 232–234. [12] O. Katz, P. Heidmann, M. Fink, S. Gigan, Non-invasive single-shot imaging through scattering layers and around corners via speckle correlations, Nat. Photon. 8 (2014) 784–790. [13] S. Feng, C. Kane, P.A. Lee, D.A. Stone, Correlations and fluctuations of coherent wave transmission through disordered media, Phys. Rev. Lett. 61 (7) (1988) 834–837.
9
[14] D. Ancora, D. Di Battista, G. Giasafaki, S. Psycharakis, E. Liapis, A. Zacharopoulos, G. Zacharakis, Optical projection tomography via phase retrieval algorithms for hidden three dimensional imaging, Proc. SPIE 10074 (2017) 100741E. [15] J.R. Walls, J.G. Sled, J. Sharpe, R.M. Henkelman, Correction of artefacts in optical projection tomography, Phys. Med. Biol. 50 (19) (2005) 4645. [16] Y. Shechtman, Y.C. Eldar, O. Cohen, H.N. Chapman, J. Miao, M. Segev, Phase retrieval with application to optical imaging, IEEE Signal Process. Mag. 32 (3) (2015) 87–109. [17] R.W. Gerchberg, W.O. Saxton, A practical algorithm for the determination of the phase from image and diffraction plane pictures, Optik 35 (1972) 237. [18] J.R. Fienup, Phase retrieval algorithms: a comparison, App. Opt. 21 (15) (1982) 2758–2769. [19] J.R. Fienup, Phase retrieval algorithms: a personal tour [Invited], Appl. Opt. 52 (1) (2013) 45–56. [20] T. Bendory, R. Beinert and Y. C. Eldar, Fourier Phase Retrieval: Uniqueness and Algorithms, arXiv, p. arXiv:1705.09590, 2017. [21] J.A. Rodriguez, R. Xu, C.-C. Chen, Y. Zou, J. Miao, Oversampling smoothness: an effective algorithm for phase retrieval of noisy diffraction intensities, J. Appl. Crystallogr. 46 (2) (2013) 312–318. [22] K. Jaganathan, S. Oymak, B. Hassibi, Sparse Phase Retrieval: Uniqueness Guarantees and Recovery Algorithms, IEEE Trans. Signal Process. 65 (9) (2017) 2402–2410. [23] A. Arranz, D. Dong, S. Zhu, M. Rudin, C. Tsatsanis, J. Tian, J. Ripoll, Helical optical projection tomography, Optics Express 21 (22) (2013) 25912–25925. [24] A. Arranz, D. Dong, S. Zhu, C. Savakis, J. Tian, J. Ripoll, In-vivo optical tomography of small scattering specimens: time-lapse 3D imaging of the head eversion process in Drosophila melanogaster, Sci. Rep. 4 (2014) 7325. [25] D. Castaño-Díez, D. Moser, A. Schoenegger, S. Pruggnaller, A.S. Frangakis, Performance evaluation of image processing algorithms on the GPU, J. Struct. Biol. 164 (2008) 153–160.
Please cite this article in press as: D. Ancora et al., Methods (2017), https://doi.org/10.1016/j.ymeth.2017.10.009