1842
OPTICS LETTERS / Vol. 40, No. 8 / April 15, 2015
Digital resampling diversity sparsity constrained-wavefield reconstruction using single-magnitude image Yair Rivenson,†,* Maya Aviv (Shalev),† Aryeh Weiss, Hana Panet, and Zeev Zalevsky Faculty of Engineering, Bar-Ilan University, Ramat Gan 52900, Israel *Corresponding author: [email protected] Received January 13, 2015; accepted March 12, 2015; posted March 18, 2015 (Doc. ID 232338); published April 13, 2015 Phase retrieval is a well-established method for the recovery of an object wave from its magnitude-only measurements, with applications in fields such as material science, biology, and astronomy. In order to guarantee a stable and global solution for phase retrieval, measurement diversity is frequently used. However, this process requires taking several measurements, which hinders some of the advantages of phase retrieval compared with interferometric techniques. Here we propose a novel diversity framework that we refer to as digital resampling diversity. The proposed resampling diversity along with sparsity-constrained object reconstruction enables improved object reconstruction from a single squared magnitude image, without using object support constraint. We demonstrate this framework with simulation and experimental results. © 2015 Optical Society of America OCIS codes: (100.5070) Phase retrieval; (100.3190) Inverse problems; (100.2000) Digital image processing. http://dx.doi.org/10.1364/OL.40.001842
Phase retrieval is defined as the numerical recovery of a complex function from its magnitude measurement only: g jΦfj2 ;
(1)
where f ∈ CN is the signal, g ∈ RN is the measurement, and Φ ∈ CN×N is an operator relating the object’s wavefield to the measurement plane wavefiled. In the classical phase retrieval setting, Φ is the Fourier transform operator [1,2], but other types of measurement configurations may also be applied such as Fresnel, Fractional Fourier transforms [3], etc. Phase retrieval is mostly known as an alternative to interferometric-based methods [4] and is successfully applied in fields such as X-ray diffraction, astronomical imaging, and crystallography [5]. Typically, noninterferometric phase retrieval relies on some a priori information in order to facilitate phase recovery. Popular choices of such a priori knowledge are the signal support, non-negativity, or signal structure [2]. The reconstruction is then performed by employing classical algorithms that iterate back and forth between the measured intensity plane and the object plane [2], such as the Gerchberg–Saxton (GS), error reduction (ER) [2], hybrid input output (HIO) [2], relaxed averaged alternating reflections (RAAR) [6], etc. These algorithms are reported to be sensitive to the initial conditions and might lead to stagnation [7] in some cases. Often, these algorithms require trying several initializations in order to get a satisfying result. One approach to promote convergence to a global solution, to avoid stagnation, and to improve noise immunity is by using measurement diversity. The more popular approach is using axial diversity [8], which is carried out manually by taking a set of images at different distances from the object plane [8] or virtually using spatial light modulator while recording at a single camera plane [9]. Other forms of diversity mechanisms have also been demonstrated as well, such as transverse displacement 0146-9592/15/081842-04$15.00/0
[10] and polarization [11]. All these methods require taking several measurements, which ultimately complicate the hardware and make them prone to stability and alignment prediction errors, thus deviating from the main advantages of noninterferometric phase retrieval over interferometric methods. It also makes the measurement process less suitable for detection of dynamic events [12]. In this Letter, we introduce a new type of diversity mechanism that we refer to as digital resampling diversity which, when combined with sparsity promoting object constraint, can be used as a diversity mechanism for phase retrieval applications. Digital resampling was used recently in order to reduce noise in numerical hologram reconstruction [13]. In a closer relation to this undertaking, a similar idea was employed to reduce noise in fluorescence imaging, by resampling its Fourier pattern, and using the fact that reconstruction from each resampled pattern yields a correlated (desired) signal and uncorrelated noise that can then be reduced by averaging [14]. Thus, the main advantage of the proposed method is that it demonstrates superior convergence while requiring a single measurement and employing the object’s sparsity constraint only. We assume that the object can be sparsely represented using some mathematical operator, Ψ, i.e., if f Ψα, α ∈ CN is an S-sparse vector, meaning that it has only S ≪ N meaningful entries. Using this sparsity assumption, we wish to use the compressive sensing paradigm [15], which provides a mathematical and algorithmic framework for the recovery of S-sparse signals from M < Nonly measurements. Previous works that applied compressive sensing to the problem of phase retrieval have shown that accurate recovery is indeed possible from M < N phaseless measurements, by using the sparsity constraint as an integral part of the recovery scheme [16–18]. Following that line, when the entire diffraction pattern is captured (N pixels), a measurement redundancy is © 2015 Optical Society of America
April 15, 2015 / Vol. 40, No. 8 / OPTICS LETTERS
obtained. This stems from the fact that the S-sparse object can be recovered from only M < N pixels. The proposed framework is shown in Fig. 1, where an N pixel diffraction pattern is digitally resampled with r random binary masks, such that each one selects only M < N pixels. This practically creates measurement diversity. Each of the r subsampled measurements is used as constraint. The solution should satisfy the r sets of constraints. The major advantage of the proposed framework is its flexibility, i.e., the diversity mechanism can be easily controlled since it is all performed digitally, on a computer. Depending on the problem’s features, the digital resampling parameters can be tuned. That is, we can control the resampling depth (number of pixels, M, to be used in each resampling mask), the probability distribution (uniform, Gaussian, etc.) from which the pixels are being resampled, and the number of resampled patterns, r, to be generated. The digital resampling can be done for numerous types of measurements operators, Φ, including Fourier and Fresnel transforms. The sparsity-promoting constraint may include several types of sparsifying operators, Ψ, such as the wavelet transform. As a recovery algorithm, basic Fienup (HIO) algorithm or its variants can be employed by simply adding the sparsity constraints to the object space [19]. Different sparsifying operators can also be applied for phase and amplitude in case that the object is complex [20]. Another important feature of the proposed framework is that it enables signal recovery by using the sparsity constraint only, which is much more flexible than traditional constraints, such as support and non-negativity, although they could be easily assimilated to the framework as well. Here, we suggest the following hybrid input output— error reduction (HIO-ER) sparsity-constrained digitally resampled algorithm. At the beginning, an HIO scheme is applied, and after a fixed number of iterations, the ER algorithm is applied in order to refine the reconstruction.
Fig. 1. Digital resampling diversity concept. A single measurement is taken (at this illustration, far field), and it is being digitally randomly resampled r times.
1843
Algorithm: Input g—the recorded diffraction pattern, S—approximated sparse signal components. Rmask—resampling masks array of length r, THR—support change threshold, maxIterHIO, maxIterER— maximum number of iterations for the HIO/ER phases, respectively. Initialize g0 ← g × expfj2πpng where pn is drawn from a uniform distribution U0; 1, f 0 ← Φ g0 , l ← 0, Λ ← ∅, supportChange ← N Main HIO For l 1:maxIterHIO for k 1:r a. gk−1 ← Φf k−1 b. gk ← Rmaskfkg × g × gk−1 ∕jgk−1 j c. f k ← Φ gk d. αk ← Ψ f k —Project to sparsifying domain. e. αΛk ← Keep largest S elements of αk and update sparse signal support Λ. f. αk Λc ← αk−1 Λc − βαk Λc —reduce error outside sparse support region. g. f k ← Ψαk end for l←l1 end while Main ER While l < maxIterER and supportChange < THR for k 1:r h. gk−1 ← Φf k−1 i. gk ← Rmaskfkg × g × gk−1 ∕jgk−1 j j. f k ← Φ gk k. αk ← Ψ f k l. αsk ← Keep largest S elements of αk m. Update sparse signal support Λ n. f k ← Ψαk end for o. update supportChange ← Λl nΛl−1 l←l1 End
In steps (a)-(b) we update the phase of the measurement space and multiply it by the kth predefined random-resampling masks (× symbolizes Hadamard product). In step (c), we project the result to the object plane, and in step (d), we project the result to its sparsifying space. In step (e), we threshold the largest S coefficients and update their support. In step (f), we regularize the error outside the updated sparse support area with the β parameter. In step (g), we re-project the result to the object space. After a pre-defined number of loops, we perform the ER phase, where steps (h)–(n) are similar to the HIO approach, expect from the fact that we perform hard thresholding of the sparse coefficients, while in the HIO phase, the thresholding is regularized by the β parameter. From our experiments with the algorithm, we have witnessed that the signal support is obtained in the HIO step while the coefficients amplitudes are refined in the ER phase. At first we wish to demonstrate the benefits of the proposed algorithm with a set of numerical experiments, using Matlab. The target object was set to be a 90 × 90 pixels Shepp–Logan phantom object, which is shown in Fig. 2(a). This particular object has centrosymmetric support, which often leads to stagnation of the solution due to the twin-image problem [21]. The object was
1844
OPTICS LETTERS / Vol. 40, No. 8 / April 15, 2015
embedded into a 128 × 128 pixels array, which can be considered as undersampling of the object’s autocorrelation function [21] (in the vertical direction). Using the Haar wavelet as a sparsifying transform, Ψ, the number of nonzeros can be well estimated as S 1, 150 elements while the true number of nonzero components was 14,260 elements, where the difference can be regarded as “noise”. Thus, we set the threshold of the number of nonzeros in the algorithm to be 3S. We then recorded our measurement g, as the modulus of the Fourier transform of the zero-padded object. Reconstruction using the full diffraction pattern and HIO-ER sparsity-constrained algorithm is shown in Fig. 2(b). In Fig. 2(c) the reconstruction using the proposed resampling method is shown, with r 40 resampling masks drawn from a uniform probability density function, such that M∕N × 100 99% of the pixels are resampled from the diffraction pattern, for each one of the r constraints. Both simulations ran the same number of iterations (150 HIO and 200 ER iterations times the number of masks). Using a look-uptable, we made sure that no pixel is being resampled by the masks twice. The results shown in Fig. 2 are typical, and the obvious advantage of the proposed method is witnessed both visually and quantitatively. We have used β 0.06 for the reconstruction from the full and resampled diffraction patterns, which we numerically investigated to promote the best result. We ran the algorithm for 30 trials (using different initial conditions) and found the proposed digital resampling diversity sparsityconstrained HIO-ER converged 29 out of 30 trials, while no convergence was witnessed for the sparsityconstrained HIO-ER. This shows that the resampling mechanism is making the difference for the algorithm’s performance. We also compared the proposed algorithm to the standard HIO-ER algorithm with object nonnegativity constraints [21] and β 0.85 [Fig. 2(d)] and again in 30 trials no convergence was witnessed. In order to evaluate the performance of the algorithm for different sparsity ratios, S∕N obj , we have created a 64 × 64 object that is randomly sparsely populated with S (random amplitude) point sources. This time, we have zero padded the object to be 130 × 130, allowing the Nyquist sampling of the autocorrelation function. The number of resampling was set to r 10, with M∕N 0.99. We again repeated the same reconstruction procedures as elaborated in the previous paragraph, only
Fig. 2. (a) Original Shepp–Logan phantom. (b) Typical reconstruction result for the sparsity constrained HIO-ER algorithm (c) Typical reconstruction result for the resampling diversity HIO-ER algorithm PSNR 30.1 dB. Both (b) and (c) reconstructions were attempted with sparsity constraints only, i.e., no support or non-negativity constraints. (d) Typical result when attempting to solve the problem standard HIO-ER algorithm with no negativity and real constraints. The PSNR was calculated after image registration.
this time the canonical basis is the sparsifying operator, Ψ I. The average result over 10 trials is summarized in Fig. 3. It is evident from Fig. 3 that for a sparsity ratio till 0.1, the proposed resampling framework allows exact signal reconstruction MSE ∼ 1 × 10−6 , with merely the object sparsity constraint. For higher values of S, the proposed framework is still superior to the sparse HIO-ER method. The algorithm was then tested on experimental data. B16 melanoma cells were grown in RPMI medium supplemented with 10% FCS and antibiotics. 0.5 × 104 cells were seeded in 35-mm cell culture dish 72 h before imaging. These cells were imaged with a 40 × ∕NA 0.6 long working distance air objective, in a Nikon TE2000E fully automated inverted microscope. The condenser was first adjusted for Kohler illumination, and then the condenser diaphragm was closed to its minimum size to increase the spatial coherence of the illumination. The transmitted light path included an HQ535/50 535 nm 25 nm bandpass filter (Chroma Technology, Bellows Falls, Vermont, USA). The images were acquired through a 1× c-mount adapter, using a QImaging Retiga 2000R camera (QImaging, Surrey, BC, Canada), which has a pixel pitch of 7.4 μm. A Z-stack of images was acquired with a spacing of 2 μm between images. An image stack of an empty field was acquired for use in suppression of coherence artifacts. A region of interest of 512 × 512 pixels, containing a single cell, was cropped from the original 1600 × 1200 image. Fixed-pattern coherence artifacts were suppressed by dividing the data image stack by the background image stack. The reconstruction was carried out as follows. First, a reference phase was found by recoding 3 diffraction patterns separated 26 μm apart [e.g., Fig. 4(a)], such that the relation between the different diffraction patterns is given by the Fresnel transform, and then iteratively propagating through these diffraction patterns. The result is shown in Fig. 4(b). This result served as a reference to evaluate the
Fig. 3. Comparison of average MSE of uniformly at random randomly sparse populated object at different object sparsity ratios, S∕N obj . The reconstruction using the sparsityconstrained HIO-ER algorithm is shown with dashed stars, while reconstruction using the digital resampling-diversity HIO-ER is shown in the solid line with circles. Asbefore, both reconstructions were attempted with sparsity constraints only.
April 15, 2015 / Vol. 40, No. 8 / OPTICS LETTERS
1845
within popular approaches for phase retrieval. The proposed method outperformed sparsity-constrained phase retrieval and standard HIO-ER algorithm in both numerical simulations and optical experiments. The main limitation of this method is the requirement that the object’s wavefield should be sparsely represented in some mathematical basis, which is a fair assumption for many natural objects. There is still some research to be performed in order to optimize the different parameters of the algorithm (including M and r), which we leave for future work. † Both authors contributed equally to this work.
Fig. 4. Reconstruction results for the B16 melanoma cells. (a) One of the defocused images stack used for the phase reconstruction. (b) Reconstruction from three axially defocused images. (c) Reconstruction using sparsity-constrained HIO-ER algorithm. (d) Reconstruction using the digital resampling-diversity HIO-ER algorithm.
algorithm’s performance. We then took the diffraction pattern shown in Fig. 4(a) and applied the sparsityconstrained HIO-ER with and without the resampling diversity constraints, where we employed the Haar wavelet as a sparsifying basis. The recovery results are shown in Figs. 4(c) and 4(d), respectively. For the resampling diversity framework, we used r 30 random subsampling patterns where 99% of the pixels were randomly chosen for each diffraction pattern. For both algorithms, 100 HIO followed by 10 ER iterations were performed, such that the total number of constraints applied was the same as for the reconstruction from the 3 recorded diffraction patterns. For both algorithms, we evaluated the number of nonzero wavelet coefficients such that the S∕N obj 0.2 and β 0.75, while constraining the object to be pure phase. The MSE for the sparsity-constrained HIO-ER was 0.158, while for the digitally resampled sparsityconstrained HIO-ER, it was 0.0763. To conclude, we have proposed a new phase-retrieval framework that enables measurement diversity by digital resampling of the recorded diffraction pattern. The framework shows promising results in the recovery of objects without using any estimation of the support constraint. The framework is flexible and can be easily embedded
References 1. R. Gerchberg and W. Saxton, Optik 35, 237 (1972). 2. J. Fienup, Appl. Opt. 21, 2758 (1982). 3. Z. Zalevsky, D. Mendlovic, and R. G. Dorsch, Opt. Lett. 21, 842 (1996). 4. T. Kreis, Handbook of Holographic Interferometry, 1st ed. (Wiley-VCH, 2004). 5. J. Fienup, Appl. Opt. 52, 45 (2013). 6. D. R. Luke, Inverse Probl. 21, 37 (2005). 7. J. Fienup and C. Wackerman, J. Opt. Soc. Am. A 3, 1897 (1986). 8. R. A. Gonsalves, Opt. Eng. 21, 215829 (1982). 9. M. Agour, P. F. Almoro, and C. Falldorf, J. Eur. Opt. Soc. 7, 12046 (2012). 10. H. M. L. Faulkner and J. M. Rodenburg, Phys. Rev. Lett. 93, 023903 (2004). 11. B. R. Hunt, T. L. Overman, and P. Gough, Opt. Lett. 23, 1123 (1998). 12. Y. Shechtman, Y. C. Eldar, O. Cohen, and M. Segev, Opt. Express 21, 6327 (2013). 13. V. Bianco, M. Paturzo, P. Memmolo, A. Finizio, P. Ferraro, and B. Javidi, Opt. Lett. 38, 619 (2013). 14. M. Marim, E. Angelini, and J.-C. Olivo-Marin, Proc. SPIE Wavelets XIII 7446, 744605 (2009). 15. E. Candes and M. Wakin, IEEE Signal Proc. Mag. 25(2), 21 (2008). 16. M. L. Moravec, J. K. Romberg, and R. G. Baraniuk, Proc. SPIE 6701, 670120 (2007). 17. W. Chan, M. Moravec, R. Baraniuk, and D. Mittleman, Opt. Lett. 33, 974 (2008). 18. P. Schniter and S. Rangan, “Compressive phase retrieval via generalized approximate message passing,” in Proceeding of Allerton Conference on Communication, Control, and Computing, Monticello, Ill., October 2012. 19. S. Mukherjee and C. S. Seelamantula, “An iterative algorithm for phase retrieval with sparsity constraints: application to frequency domain optical coherence tomography,” IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Kyoto, Japan, 2012, pp. 553–556. 20. V. Katkovnik and J. Astola, J. Opt. Soc. Am. A 29, 44 (2012). 21. M. Giuzar-Sicairos and J. R. Fienup, J. Opt. Soc. Am. A 29, 2367 (2012).
OPTICS LETTERS / Vol. 40, No. 8 / April 15, 2015
Digital resampling diversity sparsity constrained-wavefield reconstruction using single-magnitude image Yair Rivenson,†,* Maya Aviv (Shalev),† Aryeh Weiss, Hana Panet, and Zeev Zalevsky Faculty of Engineering, Bar-Ilan University, Ramat Gan 52900, Israel *Corresponding author: [email protected] Received January 13, 2015; accepted March 12, 2015; posted March 18, 2015 (Doc. ID 232338); published April 13, 2015 Phase retrieval is a well-established method for the recovery of an object wave from its magnitude-only measurements, with applications in fields such as material science, biology, and astronomy. In order to guarantee a stable and global solution for phase retrieval, measurement diversity is frequently used. However, this process requires taking several measurements, which hinders some of the advantages of phase retrieval compared with interferometric techniques. Here we propose a novel diversity framework that we refer to as digital resampling diversity. The proposed resampling diversity along with sparsity-constrained object reconstruction enables improved object reconstruction from a single squared magnitude image, without using object support constraint. We demonstrate this framework with simulation and experimental results. © 2015 Optical Society of America OCIS codes: (100.5070) Phase retrieval; (100.3190) Inverse problems; (100.2000) Digital image processing. http://dx.doi.org/10.1364/OL.40.001842
Phase retrieval is defined as the numerical recovery of a complex function from its magnitude measurement only: g jΦfj2 ;
(1)
where f ∈ CN is the signal, g ∈ RN is the measurement, and Φ ∈ CN×N is an operator relating the object’s wavefield to the measurement plane wavefiled. In the classical phase retrieval setting, Φ is the Fourier transform operator [1,2], but other types of measurement configurations may also be applied such as Fresnel, Fractional Fourier transforms [3], etc. Phase retrieval is mostly known as an alternative to interferometric-based methods [4] and is successfully applied in fields such as X-ray diffraction, astronomical imaging, and crystallography [5]. Typically, noninterferometric phase retrieval relies on some a priori information in order to facilitate phase recovery. Popular choices of such a priori knowledge are the signal support, non-negativity, or signal structure [2]. The reconstruction is then performed by employing classical algorithms that iterate back and forth between the measured intensity plane and the object plane [2], such as the Gerchberg–Saxton (GS), error reduction (ER) [2], hybrid input output (HIO) [2], relaxed averaged alternating reflections (RAAR) [6], etc. These algorithms are reported to be sensitive to the initial conditions and might lead to stagnation [7] in some cases. Often, these algorithms require trying several initializations in order to get a satisfying result. One approach to promote convergence to a global solution, to avoid stagnation, and to improve noise immunity is by using measurement diversity. The more popular approach is using axial diversity [8], which is carried out manually by taking a set of images at different distances from the object plane [8] or virtually using spatial light modulator while recording at a single camera plane [9]. Other forms of diversity mechanisms have also been demonstrated as well, such as transverse displacement 0146-9592/15/081842-04$15.00/0
[10] and polarization [11]. All these methods require taking several measurements, which ultimately complicate the hardware and make them prone to stability and alignment prediction errors, thus deviating from the main advantages of noninterferometric phase retrieval over interferometric methods. It also makes the measurement process less suitable for detection of dynamic events [12]. In this Letter, we introduce a new type of diversity mechanism that we refer to as digital resampling diversity which, when combined with sparsity promoting object constraint, can be used as a diversity mechanism for phase retrieval applications. Digital resampling was used recently in order to reduce noise in numerical hologram reconstruction [13]. In a closer relation to this undertaking, a similar idea was employed to reduce noise in fluorescence imaging, by resampling its Fourier pattern, and using the fact that reconstruction from each resampled pattern yields a correlated (desired) signal and uncorrelated noise that can then be reduced by averaging [14]. Thus, the main advantage of the proposed method is that it demonstrates superior convergence while requiring a single measurement and employing the object’s sparsity constraint only. We assume that the object can be sparsely represented using some mathematical operator, Ψ, i.e., if f Ψα, α ∈ CN is an S-sparse vector, meaning that it has only S ≪ N meaningful entries. Using this sparsity assumption, we wish to use the compressive sensing paradigm [15], which provides a mathematical and algorithmic framework for the recovery of S-sparse signals from M < Nonly measurements. Previous works that applied compressive sensing to the problem of phase retrieval have shown that accurate recovery is indeed possible from M < N phaseless measurements, by using the sparsity constraint as an integral part of the recovery scheme [16–18]. Following that line, when the entire diffraction pattern is captured (N pixels), a measurement redundancy is © 2015 Optical Society of America
April 15, 2015 / Vol. 40, No. 8 / OPTICS LETTERS
obtained. This stems from the fact that the S-sparse object can be recovered from only M < N pixels. The proposed framework is shown in Fig. 1, where an N pixel diffraction pattern is digitally resampled with r random binary masks, such that each one selects only M < N pixels. This practically creates measurement diversity. Each of the r subsampled measurements is used as constraint. The solution should satisfy the r sets of constraints. The major advantage of the proposed framework is its flexibility, i.e., the diversity mechanism can be easily controlled since it is all performed digitally, on a computer. Depending on the problem’s features, the digital resampling parameters can be tuned. That is, we can control the resampling depth (number of pixels, M, to be used in each resampling mask), the probability distribution (uniform, Gaussian, etc.) from which the pixels are being resampled, and the number of resampled patterns, r, to be generated. The digital resampling can be done for numerous types of measurements operators, Φ, including Fourier and Fresnel transforms. The sparsity-promoting constraint may include several types of sparsifying operators, Ψ, such as the wavelet transform. As a recovery algorithm, basic Fienup (HIO) algorithm or its variants can be employed by simply adding the sparsity constraints to the object space [19]. Different sparsifying operators can also be applied for phase and amplitude in case that the object is complex [20]. Another important feature of the proposed framework is that it enables signal recovery by using the sparsity constraint only, which is much more flexible than traditional constraints, such as support and non-negativity, although they could be easily assimilated to the framework as well. Here, we suggest the following hybrid input output— error reduction (HIO-ER) sparsity-constrained digitally resampled algorithm. At the beginning, an HIO scheme is applied, and after a fixed number of iterations, the ER algorithm is applied in order to refine the reconstruction.
Fig. 1. Digital resampling diversity concept. A single measurement is taken (at this illustration, far field), and it is being digitally randomly resampled r times.
1843
Algorithm: Input g—the recorded diffraction pattern, S—approximated sparse signal components. Rmask—resampling masks array of length r, THR—support change threshold, maxIterHIO, maxIterER— maximum number of iterations for the HIO/ER phases, respectively. Initialize g0 ← g × expfj2πpng where pn is drawn from a uniform distribution U0; 1, f 0 ← Φ g0 , l ← 0, Λ ← ∅, supportChange ← N Main HIO For l 1:maxIterHIO for k 1:r a. gk−1 ← Φf k−1 b. gk ← Rmaskfkg × g × gk−1 ∕jgk−1 j c. f k ← Φ gk d. αk ← Ψ f k —Project to sparsifying domain. e. αΛk ← Keep largest S elements of αk and update sparse signal support Λ. f. αk Λc ← αk−1 Λc − βαk Λc —reduce error outside sparse support region. g. f k ← Ψαk end for l←l1 end while Main ER While l < maxIterER and supportChange < THR for k 1:r h. gk−1 ← Φf k−1 i. gk ← Rmaskfkg × g × gk−1 ∕jgk−1 j j. f k ← Φ gk k. αk ← Ψ f k l. αsk ← Keep largest S elements of αk m. Update sparse signal support Λ n. f k ← Ψαk end for o. update supportChange ← Λl nΛl−1 l←l1 End
In steps (a)-(b) we update the phase of the measurement space and multiply it by the kth predefined random-resampling masks (× symbolizes Hadamard product). In step (c), we project the result to the object plane, and in step (d), we project the result to its sparsifying space. In step (e), we threshold the largest S coefficients and update their support. In step (f), we regularize the error outside the updated sparse support area with the β parameter. In step (g), we re-project the result to the object space. After a pre-defined number of loops, we perform the ER phase, where steps (h)–(n) are similar to the HIO approach, expect from the fact that we perform hard thresholding of the sparse coefficients, while in the HIO phase, the thresholding is regularized by the β parameter. From our experiments with the algorithm, we have witnessed that the signal support is obtained in the HIO step while the coefficients amplitudes are refined in the ER phase. At first we wish to demonstrate the benefits of the proposed algorithm with a set of numerical experiments, using Matlab. The target object was set to be a 90 × 90 pixels Shepp–Logan phantom object, which is shown in Fig. 2(a). This particular object has centrosymmetric support, which often leads to stagnation of the solution due to the twin-image problem [21]. The object was
1844
OPTICS LETTERS / Vol. 40, No. 8 / April 15, 2015
embedded into a 128 × 128 pixels array, which can be considered as undersampling of the object’s autocorrelation function [21] (in the vertical direction). Using the Haar wavelet as a sparsifying transform, Ψ, the number of nonzeros can be well estimated as S 1, 150 elements while the true number of nonzero components was 14,260 elements, where the difference can be regarded as “noise”. Thus, we set the threshold of the number of nonzeros in the algorithm to be 3S. We then recorded our measurement g, as the modulus of the Fourier transform of the zero-padded object. Reconstruction using the full diffraction pattern and HIO-ER sparsity-constrained algorithm is shown in Fig. 2(b). In Fig. 2(c) the reconstruction using the proposed resampling method is shown, with r 40 resampling masks drawn from a uniform probability density function, such that M∕N × 100 99% of the pixels are resampled from the diffraction pattern, for each one of the r constraints. Both simulations ran the same number of iterations (150 HIO and 200 ER iterations times the number of masks). Using a look-uptable, we made sure that no pixel is being resampled by the masks twice. The results shown in Fig. 2 are typical, and the obvious advantage of the proposed method is witnessed both visually and quantitatively. We have used β 0.06 for the reconstruction from the full and resampled diffraction patterns, which we numerically investigated to promote the best result. We ran the algorithm for 30 trials (using different initial conditions) and found the proposed digital resampling diversity sparsityconstrained HIO-ER converged 29 out of 30 trials, while no convergence was witnessed for the sparsityconstrained HIO-ER. This shows that the resampling mechanism is making the difference for the algorithm’s performance. We also compared the proposed algorithm to the standard HIO-ER algorithm with object nonnegativity constraints [21] and β 0.85 [Fig. 2(d)] and again in 30 trials no convergence was witnessed. In order to evaluate the performance of the algorithm for different sparsity ratios, S∕N obj , we have created a 64 × 64 object that is randomly sparsely populated with S (random amplitude) point sources. This time, we have zero padded the object to be 130 × 130, allowing the Nyquist sampling of the autocorrelation function. The number of resampling was set to r 10, with M∕N 0.99. We again repeated the same reconstruction procedures as elaborated in the previous paragraph, only
Fig. 2. (a) Original Shepp–Logan phantom. (b) Typical reconstruction result for the sparsity constrained HIO-ER algorithm (c) Typical reconstruction result for the resampling diversity HIO-ER algorithm PSNR 30.1 dB. Both (b) and (c) reconstructions were attempted with sparsity constraints only, i.e., no support or non-negativity constraints. (d) Typical result when attempting to solve the problem standard HIO-ER algorithm with no negativity and real constraints. The PSNR was calculated after image registration.
this time the canonical basis is the sparsifying operator, Ψ I. The average result over 10 trials is summarized in Fig. 3. It is evident from Fig. 3 that for a sparsity ratio till 0.1, the proposed resampling framework allows exact signal reconstruction MSE ∼ 1 × 10−6 , with merely the object sparsity constraint. For higher values of S, the proposed framework is still superior to the sparse HIO-ER method. The algorithm was then tested on experimental data. B16 melanoma cells were grown in RPMI medium supplemented with 10% FCS and antibiotics. 0.5 × 104 cells were seeded in 35-mm cell culture dish 72 h before imaging. These cells were imaged with a 40 × ∕NA 0.6 long working distance air objective, in a Nikon TE2000E fully automated inverted microscope. The condenser was first adjusted for Kohler illumination, and then the condenser diaphragm was closed to its minimum size to increase the spatial coherence of the illumination. The transmitted light path included an HQ535/50 535 nm 25 nm bandpass filter (Chroma Technology, Bellows Falls, Vermont, USA). The images were acquired through a 1× c-mount adapter, using a QImaging Retiga 2000R camera (QImaging, Surrey, BC, Canada), which has a pixel pitch of 7.4 μm. A Z-stack of images was acquired with a spacing of 2 μm between images. An image stack of an empty field was acquired for use in suppression of coherence artifacts. A region of interest of 512 × 512 pixels, containing a single cell, was cropped from the original 1600 × 1200 image. Fixed-pattern coherence artifacts were suppressed by dividing the data image stack by the background image stack. The reconstruction was carried out as follows. First, a reference phase was found by recoding 3 diffraction patterns separated 26 μm apart [e.g., Fig. 4(a)], such that the relation between the different diffraction patterns is given by the Fresnel transform, and then iteratively propagating through these diffraction patterns. The result is shown in Fig. 4(b). This result served as a reference to evaluate the
Fig. 3. Comparison of average MSE of uniformly at random randomly sparse populated object at different object sparsity ratios, S∕N obj . The reconstruction using the sparsityconstrained HIO-ER algorithm is shown with dashed stars, while reconstruction using the digital resampling-diversity HIO-ER is shown in the solid line with circles. Asbefore, both reconstructions were attempted with sparsity constraints only.
April 15, 2015 / Vol. 40, No. 8 / OPTICS LETTERS
1845
within popular approaches for phase retrieval. The proposed method outperformed sparsity-constrained phase retrieval and standard HIO-ER algorithm in both numerical simulations and optical experiments. The main limitation of this method is the requirement that the object’s wavefield should be sparsely represented in some mathematical basis, which is a fair assumption for many natural objects. There is still some research to be performed in order to optimize the different parameters of the algorithm (including M and r), which we leave for future work. † Both authors contributed equally to this work.
Fig. 4. Reconstruction results for the B16 melanoma cells. (a) One of the defocused images stack used for the phase reconstruction. (b) Reconstruction from three axially defocused images. (c) Reconstruction using sparsity-constrained HIO-ER algorithm. (d) Reconstruction using the digital resampling-diversity HIO-ER algorithm.
algorithm’s performance. We then took the diffraction pattern shown in Fig. 4(a) and applied the sparsityconstrained HIO-ER with and without the resampling diversity constraints, where we employed the Haar wavelet as a sparsifying basis. The recovery results are shown in Figs. 4(c) and 4(d), respectively. For the resampling diversity framework, we used r 30 random subsampling patterns where 99% of the pixels were randomly chosen for each diffraction pattern. For both algorithms, 100 HIO followed by 10 ER iterations were performed, such that the total number of constraints applied was the same as for the reconstruction from the 3 recorded diffraction patterns. For both algorithms, we evaluated the number of nonzero wavelet coefficients such that the S∕N obj 0.2 and β 0.75, while constraining the object to be pure phase. The MSE for the sparsity-constrained HIO-ER was 0.158, while for the digitally resampled sparsityconstrained HIO-ER, it was 0.0763. To conclude, we have proposed a new phase-retrieval framework that enables measurement diversity by digital resampling of the recorded diffraction pattern. The framework shows promising results in the recovery of objects without using any estimation of the support constraint. The framework is flexible and can be easily embedded
References 1. R. Gerchberg and W. Saxton, Optik 35, 237 (1972). 2. J. Fienup, Appl. Opt. 21, 2758 (1982). 3. Z. Zalevsky, D. Mendlovic, and R. G. Dorsch, Opt. Lett. 21, 842 (1996). 4. T. Kreis, Handbook of Holographic Interferometry, 1st ed. (Wiley-VCH, 2004). 5. J. Fienup, Appl. Opt. 52, 45 (2013). 6. D. R. Luke, Inverse Probl. 21, 37 (2005). 7. J. Fienup and C. Wackerman, J. Opt. Soc. Am. A 3, 1897 (1986). 8. R. A. Gonsalves, Opt. Eng. 21, 215829 (1982). 9. M. Agour, P. F. Almoro, and C. Falldorf, J. Eur. Opt. Soc. 7, 12046 (2012). 10. H. M. L. Faulkner and J. M. Rodenburg, Phys. Rev. Lett. 93, 023903 (2004). 11. B. R. Hunt, T. L. Overman, and P. Gough, Opt. Lett. 23, 1123 (1998). 12. Y. Shechtman, Y. C. Eldar, O. Cohen, and M. Segev, Opt. Express 21, 6327 (2013). 13. V. Bianco, M. Paturzo, P. Memmolo, A. Finizio, P. Ferraro, and B. Javidi, Opt. Lett. 38, 619 (2013). 14. M. Marim, E. Angelini, and J.-C. Olivo-Marin, Proc. SPIE Wavelets XIII 7446, 744605 (2009). 15. E. Candes and M. Wakin, IEEE Signal Proc. Mag. 25(2), 21 (2008). 16. M. L. Moravec, J. K. Romberg, and R. G. Baraniuk, Proc. SPIE 6701, 670120 (2007). 17. W. Chan, M. Moravec, R. Baraniuk, and D. Mittleman, Opt. Lett. 33, 974 (2008). 18. P. Schniter and S. Rangan, “Compressive phase retrieval via generalized approximate message passing,” in Proceeding of Allerton Conference on Communication, Control, and Computing, Monticello, Ill., October 2012. 19. S. Mukherjee and C. S. Seelamantula, “An iterative algorithm for phase retrieval with sparsity constraints: application to frequency domain optical coherence tomography,” IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Kyoto, Japan, 2012, pp. 553–556. 20. V. Katkovnik and J. Astola, J. Opt. Soc. Am. A 29, 44 (2012). 21. M. Giuzar-Sicairos and J. R. Fienup, J. Opt. Soc. Am. A 29, 2367 (2012).
