Improved success rate and stability for phase retrieval by including randomized overrelaxation in the hybrid input output algorithm

doi:10.1364/OE.20.017093

« Show journal navigation

Improved success rate and stability for phase retrieval by including randomized overrelaxation in the hybrid input output algorithm

Martin Köhl, A. A. Minkevich, and Tilo Baumbach »View Author Affiliations

Optics Express, Vol. 20, Issue 15, pp. 17093-17106 (2012)
http://dx.doi.org/10.1364/OE.20.017093

View Full Text Article

Acrobat PDF (919 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. The phase retrieval problem
2. Reconstruction algorithms
3. Monitoring convergence of the reconstruction
4. Numerical examples
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (29)
Cited By
Figures (8)
Metrics

Abstract

In this paper, we study the success rate of the reconstruction of objects of finite extent given the magnitude of its Fourier transform and its geometrical shape. We demonstrate that the commonly used combination of the hybrid input output and error reduction algorithm is significantly outperformed by an extension of this algorithm based on randomized overrelaxation. In most cases, this extension tremendously enhances the success rate of reconstructions for a fixed number of iterations as compared to reconstructions solely based on the traditional algorithm. The good scaling properties in terms of computational time and memory requirements of the original algorithm are not influenced by this extension.

Inverse problems play an important role in many areas of physics. In an inverse problem we intend to reconstruct physical properties (“scattering potential”) of some object given the result of the interaction (“scattering”) with a well-defined input signal. Typically, propagating wave fields (e.g., acoustic waves, electromagnetic waves, . . .) are used to scan the object. However, for rapidly oscillating electromagnetic wave fields like x-rays, detection devices are only capable of recording the average intensity of the scattered wave field. Its phase information, which also contains important information about the scattering potential, is lost during the measurement process. Despite the lost phase information, x-rays are an important, non-destructive tool for the investigation of structural properties of matter because they exhibit high penetration depth and high spatial resolution [1

F. v. d. Veen and F. Pfeiffer, “Coherent x-ray scattering,” Phys J..: Condens. Matter 16, 5003–5030 (2004). [CrossRef]

–3

U. Pietsch, V. Holy, and T. Baumbach, High-Resolution X-ray Scattering From Thin Films to Lateral Nanostructures (Springer, New York, 2004).

]. Note, that 3^rd generation synchrotrons [4

G. N. Afanasiev, Vavilov-Cherenkov and Synchrotron Radiation: Foundations and Applications (Springer, Netherlands, 2004).

] provide brilliant x-ray beams with a very high degree of coherence and with very high flux, so that even dot- and wire-like nanostructures can be investigated [2

R. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Am. A 7, 394–411 (1990). [CrossRef]

I. Robinson and R. Harder, “Coherent x-ray diffraction imaging of strain at the nanoscale,” Nat Mater 8, 291–298 (2009). [CrossRef] [PubMed]

–11

A. A. Minkevich, E. Fohtung, T. Slobodskyy, M. Riotte, D. Grigoriev, M. Schmidbauer, A. C. Irvine, V. Novák, V. Holý, and T. Baumbach, “Selective coherent x-ray diffractive imaging of displacement fields in (Ga,Mn)As/GaAs periodic wires,” Phys. Rev. B 84, 054113 (2011). [CrossRef]

In this paper, we introduce a new algorithm – the (HIO+O_R)+ER-algorithm – which is capable of reconstructing the missing phase information in a variety of cases for which the traditionally and commonly used HIO+ER-algorithm [12

J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]

, 13

J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A 4, 118–123 (1986). [CrossRef]

] (reviewed in sec. 2) fails. Before we introduce this new algorithm in sec. 2, we shortly review the phase retrieval problem in sec. 1. In particular, we want to stress the fact that the (HIO+O_R)+ER-algorithm is based on randomized overrelaxation [14

A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984). [CrossRef]

,15

D. C. Youla and H. Webb, “Image restoration by the method of convex projections: part 1 – theory,” IEEE Trans. Med. Imaging 1, 81 –94 (1982). [CrossRef] [PubMed]

] and, thus, does not require additional mathematical constraints [16

A. A. Minkevich, T. Baumbach, M. Gailhanou, and O. Thomas, “Applicability of an iterative inversion algorithm to the diffraction patterns from inhomogeneously strained crystals,” Phys. Rev. B 78, 174110 (2008). [CrossRef]

] compared to the HIO+ER-algorithm. We reformulate our approach using projection polynomials which facilitate a straight forward investigation of randomization of the HIO-algorithm without correlations in the coefficients enforced by overrelaxation. Section 3 discusses numerical aspects how to measure the convergence of phase retrieval. Finally, sec. 4 illustrates the advantages of the (HIO+O_R)+ER-algorithm with respect to the traditional HIO+ER-algorithms in detail. For this purpose, we focus on three numerical examples each with different specific properties and compare the success rate of reconstructions based on our (HIO+O_R)+ER-algorithm and on the HIO+ER-algorithm. Moreover, we investigate the sensitivity of our algorithm to the choice of its parameters. In this paper, we focus on the algorithm itself and its improved convergence properties for specific, well-known problems related to phase retrieval based on the HIO+ER-algorithm. A separate publication will illustrate the improved features of the (HIO+O_R)+ER-algorithm for the reconstruction of the displacement field in inhomogeneously strained nanocrystals.

1. The phase retrieval problem

Consider some unknown object f_in(x⃗) in d-dimensional direct space (“position space”). We assume that (i) the shape 𝕊 of f_in(x⃗) (i.e., the smallest domain in direct space for which f_in(x⃗) ≠ 0) and (ii) the quantity A_k⃗ ≡ |FT{f_in(x⃗)}| (i.e., the amplitudes of its Fourier transform) are known. It has been proven that this information is sufficient to reconstruct the object f_in(x⃗) if the shape 𝕊 is finite and the dimension d of direct space is at least two [2

R. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Am. A 7, 394–411 (1990). [CrossRef]

, 17

R. P. Millane, “Multidimensional phase problems,” J. Opt. Soc. Am. A 13, 725–734 (1996). [CrossRef]

–19

R. H. T. Bates, “Fourier phase problems are uniquely solvable in more than one dimension,” Optik (Stuttgart) 61, 247–262 (1982).

]. Physically, the second condition is fulfilled in scattering experiments which can be described in lowest order Born-approximation (also called kinematical approximation). This is typically true for nanostructures illuminated by coherent x-ray beams.

Suppose we discretize a rectangular region 𝔻 ⊂ ℝ^d in direct space, which fully contains our object f_in(x⃗) (i.e., 𝕊 ⊂ 𝔻), with N_i, i ∈ {1, 2,...,d}, equidistant points. The discrete Fourier transform (DFT) maps this set of points to a mesh 𝕄 in reciprocal space which is related to the Fourier transform of our continuous direct space object f_in(x⃗) [20

L. Auslander and F. A. Grunbaum, “The Fourier transform and the discrete Fourier transform,” Inverse Prob. 5, 149 (1989). [CrossRef]

]. Hence, we have

f_{in} (\vec{x}) = \sum_{\vec{k} \in 𝕄} [\prod_{j = 1}^{d} exp (\frac{2 π i}{N_{j}} k_{j} x_{j})] A_{\vec{k}} exp (i Φ_{\vec{k}}) \forall \vec{x} \notin 𝕊,

(1)

where x⃗ = (x₁, x₂,..., x_d)^T corresponds to the equidistantly sampled points in direct space defined by 𝔻 and Φ_k⃗ are the phases associated with the amplitudes A_k⃗ in reciprocal space.

In direct space, every point x⃗ inside the object domain 𝕊 corresponds to (up to) two real unknowns (magnitude and phase). However, outside the domain 𝕊, f_in(x⃗) is precisely zero for all positions x⃗. Hence, each point outside 𝕊 defines two real equations for the determination of the N real phases Φ_k⃗, if we require f_in(x⃗) = 0 ∀ x⃗ ∉ 𝕊. Therefore, we can hope to reconstruct the missing phases Φ_k⃗ on the sampled grid 𝕄 once

2 (N - N_{𝕊}) \geq N \leftrightarrow N \geq 2 N_{𝕊},

(2)

where N_𝕊 is the number of points inside the region 𝕊. The fraction σ = N/N_𝕊 is called oversampling ratio. From our elementary discussion, we can conclude that σ = 2 is a lower bound for a successful reconstruction, if no additional a priori knowledge is available. For more details on oversampling, we refer the reader to the investigations of Elser et Milliane [21

V. Elser and R. P. Millane, “Reconstruction of an object from its symmetry-averaged diffraction pattern,” Acta Crystallogr., Sect. A: Found. Crystallogr. 64, 273–279 (2008). [CrossRef]

] and of Miao et al. [22

J. Miao, D. Sayre, and H. N. Chapman, “Phase retrieval from the magnitude of the Fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A 15, 1662–1669 (1998). [CrossRef]

] (and the references therein).

With that knowledge, we restate the phase retrieval problem more precisely: Given the shape 𝕊 of some object f_in(x⃗) and the modulus A_k⃗ of the Fourier transform of f_in(x⃗) on a mesh 𝕄 that satisfies Eq. (2), we try to reconstruct the phases Φ_k⃗ and, thus, the object f_in(x⃗).

Note, that the reconstruction is only unique up to the unavoidable inherent symmetries of the Fourier transform [2

R. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Am. A 7, 394–411 (1990). [CrossRef]

, 17

R. P. Millane, “Multidimensional phase problems,” J. Opt. Soc. Am. A 13, 725–734 (1996). [CrossRef]

, 18

J. H. Seldin and J. R. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” J. Opt. Soc. Am. A 7, 412–427 (1990). [CrossRef]

], i.e., a global phase shift (typically irrelevant), a plane wave modulation e^ik⃗·x⃗ in reciprocal space (fixed by the position of the shape 𝕊 inside 𝔻) and the degeneracy between the object f_x⃗ and its complex conjugated, inversion symmetric object

f_{\vec{x}}^{*}

. This last ambiguity constitutes a severe problem if the shape 𝕊 is inversion symmetric, but the object f_in does not possess this symmetry [13

J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A 4, 118–123 (1986). [CrossRef]

,23

J. Fienup and C. Wackermann, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3 (1986). [CrossRef]

]. For the remainder of this paper, we restrict to domains 𝕊 which are not inversion symmetric. However, even if we exclude objects with inversion symmetric shape, it is highly non-trivial to formulate a phase retrieval algorithm that works automatically without manual supervision and tuning of parameters during the reconstruction. The (HIO+O_R)+ER-algorithm which we propose in the next section is capable of performing this task for a large class of objects f_in.

2. Reconstruction algorithms

In this section, we will shortly review the traditional combination of the hybrid input output and error reduction (HIO+ER) algorithm [12

J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]

, 13

J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A 4, 118–123 (1986). [CrossRef]

] and propose our extension to this algorithm. Throughout the entire reconstruction process, we assume that we know the exact geometrical shape 𝕊 of the object which we reconstruct.

2.1. Hybrid input output and error reduction

A fast and efficient iterative phase retrieval algorithm is based on alternating projections of a trial solution onto the constraints in direct space and in reciprocal space. Hence, a single iterative step (i) of this algorithm, which is called error reduction (ER) algorithm [12

J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]

], is defined by

f_{\vec{x}}^{(i + 1)} = P_{𝕊} P_{A} f_{\vec{x}}^{(i)} \equiv {\hat{H}}_{ER} f_{\vec{x}}^{(i)} .

(3)

Typically, initial phases

Φ_{\vec{k}}^{(0)}

are chosen randomly. Here, P_𝕊 and P_A are projection operators in direct space and reciprocal space respectively, i.e.,

P_{𝕊} f_{\vec{x}}^{(i)} = {\begin{array}{l} f_{\vec{x}}^{(i)} & if \vec{x} \in 𝕊, \\ 0 & if \vec{x} \notin 𝕊 \end{array}

(4a)

and

P_{A} g_{\vec{k}}^{(i)} = A_{\vec{k}} exp (iarg (g_{\vec{k}}^{(i)})) .

(4b)

For notational simplicity, we assume that any operand of P_A or P_𝕊 is transformed to the proper space (i.e., Fourier transform to reciprocal or direct space respectively) before the operator P_A or P_𝕊 is applied (e.g.,

P_{A} f_{\vec{x}}^{(i)}

means

P_{A} FT {f_{\vec{x}}^{(i)}}

The projection operator P_A is non-linear, non-convex and non-unique [24

H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Phase retrieval, error reduction algorithm, and fienup variants: a view from convex optimization,” J. Opt. Soc. Am. A 19, 1334–1345 (2002). [CrossRef]

]. The values

| g_{\vec{k}}^{(i)} |

are important for determining the difference to the (experimentally accessible) input data A_k⃗ (see Eq. (18)). The ER-algorithm is a local minimizer of a suitable chosen error metrics [12

J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]

,13

J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A 4, 118–123 (1986). [CrossRef]

]. However, in practice, phase retrieval problems involve many local minima. Hence, the ER-algorithm will in general not converge to the correct solution f_in(x⃗), but to some local minimum. In addition, the error metric may stagnate for many iterations before decreasing further. More information on stagnation problems and the ER-algorithm can be found in [12

J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]

,13

J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A 4, 118–123 (1986). [CrossRef]

,23

J. Fienup and C. Wackermann, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3 (1986). [CrossRef]

Due to these problems, further algorithms, which try to avoid stagnation and convergence to local minima, have been developed. A very important algorithm is the hybrid input output (HIO) algorithm proposed by Fienup [12

J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]

, 13

J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A 4, 118–123 (1986). [CrossRef]

, 23

J. Fienup and C. Wackermann, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3 (1986). [CrossRef]

]. This algorithm is also an iterative procedure which is defined by the map

f_{\vec{x}}^{(i + 1)} = {\begin{array}{l} P_{A} f_{\vec{x}}^{(i)} & if \vec{x} \in 𝕊, \\ f_{\vec{x}}^{(i)} - β P_{A} f_{\vec{x}}^{(i)} & if \vec{x} \notin 𝕊, \end{array}

(5)

where β is a real parameter typically chosen in the range [0.5; 1.0] [23

J. Fienup and C. Wackermann, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3 (1986). [CrossRef]

]. Note, that this definition can be rewritten as

f_{\vec{x}}^{(i + 1)} = [1 - P_{𝕊} - β P_{A} + (1 + β) P_{𝕊} P_{A}] f_{\vec{x}}^{(i)} \equiv {\hat{H}}_{HIO} (β) f_{\vec{x}}^{(i)} .

(6)

The convergence properties of the HIO-algorithm have been investigated in more detail in [22

J. Miao, D. Sayre, and H. N. Chapman, “Phase retrieval from the magnitude of the Fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A 15, 1662–1669 (1998). [CrossRef]

,25

H. Takajo, T. Takahashi, R. Ueda, and M. Taninaka, “Study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A 15, 2849–2861 (1998). [CrossRef]

,26

H. Takajo, T. Takahashi, and T. Shizuma, “Further study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A 16, 2163–2168 (1999). [CrossRef]

This algorithm is much stronger in avoiding stagnation and local minima. However, the HIO-algorithm shows its full potential only if it is combined with the ER-algorithm. One step of this combined HIO+ER-algorithm consists of N_HIO iterations of the HIO-algorithm followed by N_ER iterations of the ER-algorithm. This combination is more successful than both algorithms on their own as it was observed in practice [13

J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A 4, 118–123 (1986). [CrossRef]

, 23

J. Fienup and C. Wackermann, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3 (1986). [CrossRef]

]. Although this algorithm is already very powerful, it is not yet satisfactory for many objects as we will show in sec. 4.

2.2. Hybrid input output with randomized overrelaxation ((HIO+O_R)+ER)

The HIO+ER-algorithm can be further improved by including randomized overrelaxation [14

A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984). [CrossRef]

, 15

D. C. Youla and H. Webb, “Image restoration by the method of convex projections: part 1 – theory,” IEEE Trans. Med. Imaging 1, 81 –94 (1982). [CrossRef] [PubMed]

,27

V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A 20, 40–55 (2003). [CrossRef]

]. Basically, overrelaxation corresponds to replacing a projection operator P_μ by

Q_{μ; λ_{μ}} \equiv 1 + λ_{μ} (P_{μ} - 1),

(7)

where λ_μ is a real constant called relaxation parameter [14

A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984). [CrossRef]

]. The limiting case Q_{μ;λ_μ} → P_μ is obtained for λ_μ → 1.

Overrelaxation (without any randomization) has been investigated for convex problems [15

D. C. Youla and H. Webb, “Image restoration by the method of convex projections: part 1 – theory,” IEEE Trans. Med. Imaging 1, 81 –94 (1982). [CrossRef] [PubMed]

] and in connection with the ER-algorithm for phase retrieval [14

A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984). [CrossRef]

]. Moreover, overrelaxation is also included in the difference map algorithm proposed by Elser in [27

V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A 20, 40–55 (2003). [CrossRef]

]. In the difference map algorithm with overrelaxation, the iterative step for our set of constraints is given by [28

S. Marchesini, “Invited article: A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007). [CrossRef] [PubMed]

]

f_{\vec{x}}^{(i + 1)} = [1 + β (P_{𝕊} Q_{A; λ_{A}} - P_{A} Q_{𝕊; λ_{𝕊}})] f_{\vec{x}}^{(i)} \equiv {\hat{H}}_{Diff} (β, λ_{A}, λ_{𝕊}) f_{\vec{x}}^{(i)},

(8)

where Elser proposes to choose the parameters as λ_A = λ_𝕊 = β⁻¹ [27

V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A 20, 40–55 (2003). [CrossRef]

]. A very nice comparison of several phase retrieval algorithms and a benchmark thereof can be found in [28

S. Marchesini, “Invited article: A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007). [CrossRef] [PubMed]

]. Marchesini demonstrated in [29

S. Marchesini, “Phase retrieval and saddle-point optimization,” J. Opt. Soc. Am. A 24, 3289–3296 (2007). [CrossRef]

] that for his particular numerical example an additional low-dimensional subspace saddle-point optimization was also able to overcome stagnation of the traditional HIO+ER-algorithm. In our extension, such an additional optimization is not necessary.

We propose to replace the (non-linear and non-convex) projection operator P_A in reciprocal space in the HIO-algorithm itself by its relaxed expression

Q_{A; λ_{A}} = 1 + λ_{A} (P_{A} - 1) .

(9)

The direct space assembly in Eq. (5) of the HIO-algorithm remains unchanged. This corresponds to changing Eq. (6) to

f_{\vec{x}}^{(i + 1)} = [1 - P_{𝕊} - β Q_{A; λ_{A}} + (1 + β) P_{𝕊} Q_{A; λ_{A}}] f_{\vec{x}}^{(i)} \equiv {\hat{H}}_{{HIO + O}_{R}} (β, λ_{A}) f_{\vec{x}}^{(i)} .

(10)

If we rewrite this expression in the projectors P_𝕊 and P_A, we get

{\hat{H}}_{HIO + O_{R}} (β, λ_{A}) \equiv [1 + β (λ_{A} - 1)] + [β - λ_{A} - β λ_{A}] P_{𝕊} - β λ_{A} P_{A} + [(1 + β) λ_{A}] P_{𝕊} P_{A} .

(11)

The deviation β (λ_A − 1) from the identity operator in the first term can neither be represented by the traditional HIO-algorithm for any value of β (see Eq. (6)) nor by the difference map algorithm for any (β, λ_A, λ_𝕊) (see Eq. (8)). In both cases, the previous iterative result

f_{\vec{x}}^{(i)}

is weighted with 1 or projected at least once in either direct (P_𝕊) or reciprocal space (P_A) before being included in the calculation for

f_{\vec{x}}^{(i + 1)}

Finally, we have to choose suitable values for the additional parameter λ_A in Eq. (9). If we restrict to a fixed set of values λ_A for all iterations, we risk simply exchanging one trap by another or one local minimum by another. Trying to overcome local minima and stagnation, we propose a new approach: For each iteration, we randomly select the parameter λ_A. More precisely, for the remainder of this paper, we investigate a uniform random distribution in the range [1 − ν, 1 +ν], ν ≥ 0, for λ_A which is reassigned each iteration. Unless stated otherwise, we choose ν = 0.5. Note, that a fixed value λ_A ≡ 1 for all iterations corresponds to the usual HIO-algorithm. In order to distinguish our new extension from the traditional HIO-algorithm, we call our extension HIO+O_R-algorithm. The power of this extension is strongly supported by the fact, that the HIO+O_R-algorithm is capable of reconstructing objects without including the ER-algorithm with significant success rates. If the HIO+O_R-algorithm is combined with ER, we call the algorithm (HIO+O_R)+ER-algorithm (see Fig. 1).

Fig. 1 Graphical illustration of the (HIO+O_R)+ER-algorithm and its building blocks.

If we look at Eq. (11) from a bird’s eye view, the most important difference of the operator Ĥ_{HIO+O_R} (β, λ_A) to the other algorithms which have been proposed [28

S. Marchesini, “Invited article: A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007). [CrossRef] [PubMed]

] is the fact, that the traps and tunnels [14

A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984). [CrossRef]

] of the operator in general depend on λ_A and, therefore, should change from iteration to iteration as a consequence of randomization. Thus, stagnation in a specific trap or tunnel is strongly reduced, as most traps and tunnels are no longer persistent throughout the iterative reconstruction process. However, the true solution f_in(x⃗) is a fixed point of the iterative procedure defined in Eq. (11) for all values of λ_A.

2.3. Projection polynomials

The extension of the HIO+ER-algorithm to the (HIO+O_R)+ER-algorithm is based on two ingredients: overrelaxation and randomization. In order to gain a deeper understanding of the (HIO+O_R)+ER-algorithm, it is useful to investigate both ingredients separately. Studying the HIO-algorithm extended by overrelaxation of P_A, but without randomization is straight forward: we simply keep a fixed, predefined overrelaxation parameter λ_A throughout the entire iterative reconstruction process. In this section, we present a framework for randomization of iterative projection algorithms performed in a way which is not based on overrelaxation of P_A. Although this framework can be applied for studying randomization of any iterative projection algorithm, we limit our numerical simulations (in section 4) to parameter values whose deterministic contribution coincides the traditional HIO-algorithm.

For this purpose, consider the projection polynomial operator

\begin{array}{l} {\hat{H}}_{Proj} (b, {\vec{c}}_{𝕊}, {\vec{c}}_{A}) = b 1 + \sum_{n = 1}^{n_{Max}^{even}} [c_{𝕊, 2 n} {(P_{𝕊} P_{A})}^{n} + c_{A, 2 n} {(P_{A} P_{𝕊})}^{n}] + \\ \sum_{n = 0}^{n_{Max}^{odd}} [c_{𝕊, 2 n + 1} P_{𝕊} {(P_{A} P_{𝕊})}^{n} + c_{A, 2 n + 1} P_{A} {(P_{𝕊} P_{A})}^{n}], \end{array}

(12a)

for obtaining the next iterative solution in the manner of Eq. (6) or Eq. (11).

n_{Max}^{even}

is given by the largest integer smaller or equal to

\frac{p}{2}

n_{Max}^{odd}

is determined by the largest integer smaller or equal to

\frac{p - 1}{2}

. p is the maximum number of successive projection operators which should be included in Ĥ_Proj.

This operator exploits the fact, that any combination of projection operators

P_{ξ_{1}}^{n_{1}} P_{ξ_{2}}^{n_{2}} \dots P_{ξ_{m}}^{n_{m}}

of two different kinds ξ_j ∈ {𝕊,A}, j ∈ {1,...,m}, reduces to one of the four building blocks P_𝕊 (P_AP_𝕊)ⁿ, P_A (P_𝕊P_A)ⁿ, (P_𝕊P_A)ⁿ and (P_AP_𝕊)ⁿ for some integer n ≥ 0. This is a consequence of the defining property

P_{ξ}^{2} = P_{ξ}

(idempotence) of projection operators (which implies

P_{ξ}^{n} = P_{ξ} \forall n \geq 1

). Note, that the product of two non-commutative idempotent operators (like P_A and P_𝕊) is no longer idempotent. The special case of the identity operator 1 and a single projection operator is included for the case n = 0 in those building blocks. However, the identity operator is included twice, namely for n = 0 for (P_𝕊P_A)ⁿ and (P_AP_𝕊)ⁿ. Therefore, one spurious coefficient in the projection polynomial is introduced if both building blocks are included for n = 0. This is avoided by restricting those two building blocks to n ≥ 1 and including the identity operator separately. Therefore, Eq. (12a) represents the most general polynomial expression which contains a maximum of q successive projections of two different kinds, if we incorporate the idempotence of projection operators.

In order to guarantee, that the true solution f_in(x⃗) is a fixed point of Ĥ_Proj for any choice of its parameters (b, c⃗_𝕊, c⃗_A), it is necessary to enforce one constraint on the coefficients: If the input for the projection operators P_𝕊 and P_A coincides with the true solution f_in(x⃗), the projection operators P_𝕊 and P_A reduce to the identity operator 1. Therefore, given the true solution f_in(x⃗) as input, the operator Ĥ_Proj simplifies to the identity operator 1 for any choice of its parameters (b, c⃗_𝕊, c⃗_A) if and only if the additional constraint

b = 1 - \sum_{n = 1}^{p} [c_{n, 𝕊} + c_{n, A}] .

(12b)

is fulfilled. Consequently, for a maximum of p successive projections in the projection polynomial Ĥ_Proj, only 2p free parameters appear in Ĥ_Proj by incorporating the idempotence of projection operators.

Given this fix point property for the true solution f_in(x⃗), we can perform proper randomization of this operator. For this purpose, we split our coefficients c_ξ,n, ξ ∈ {𝕊, A}, n ≥ 1, in a deterministic part

c_{ξ, n}^{(D)}

and a random part

r_{ξ, n} c_{ξ, n}^{(R)}

and assume r_ξ_,n to be uniformly distributed in the range [−1, 1], i.e.,

c_{ξ, n} = c_{ξ, n}^{(D)} + r_{ξ, n} c_{ξ, n}^{(R)} .

(13)

If we set some coefficients

c_{ξ, n}^{(R)} \neq 0

and draw statistically independent random values for r_ξ_,n, we can investigate randomization of the HIO-algorithm without the correlations in the coefficients c_ξ,n introduced by overrelaxation. b is fixed by Eq. (12b) and, thus, is fully correlated to the values of c_ξ,n.

In order to obtain the traditional HIO-algorithm in the limit

c_{ξ, n}^{(R)} \to 0

for all ξ and n, we have to choose (at least) q = 2 (see Eq. (6)). Therefore, Ĥ_Proj reduces to

{\hat{H}}_{{HIO}_{R}} (b, c_{𝕊, 1}, c_{𝕊, 1}, c_{A, 1}, c_{A, 2}) \equiv b 1 + c_{𝕊, 1} P_{𝕊} + c_{A, 1} P_{A} + c_{𝕊, 2} P_{𝕊} P_{A} + c_{A, 2} P_{A} P_{𝕊}

(14)

with the additional constraint

b = 1 - \sum_{n = 1}^{2} [c_{n, 𝕊} + c_{n, A}]

. The deterministic contribution

c_{ξ, n}^{(D)}

reproduces the traditional HIO-algorithm (see Eq. (6)) if we set

c_{𝕊, 1}^{(D)} = - 1, c_{A, 1}^{(D)} = - β, c_{𝕊, 2}^{(D)} = 1 + β, c_{A, 2}^{(D)} = 0 .

(15)

After replacing P_A by Q_{A;λ_A} (see Eq. (9)), the coefficients c_ξ_,n are parametrized by two parameters β and λ_A as

c_{𝕊, 1} = β - λ_{A} - β λ_{A}, c_{A, 1} = - β λ_{A}, c_{𝕊, 2} = (1 + β) λ_{A}, c_{A, 2} = 0,

(16a)

which corresponds to

c_{𝕊, 1} = - 1 - γ_{A} (1 + β), c_{A, 1} = - β (1 + γ_{A}), c_{𝕊, 2} = (1 + β) (1 + γ_{A}), c_{A, 2} = 0,

(16b)

if we substitute λ_A = 1 + γ_A. γ_A is uniformly distributed in [−ν, ν]. The deterministic contribution in Eq. (16b) is identical to Eq. (15) and reproduces the traditional HIO-algorithm in the limit ν → 0. However, the random contribution of each coefficient c_ξ_,n is determined solely by the value of γ_A and, thus, not statistically independent, but fully correlated. This is in strong contrast to a statistically independent randomization of each coefficient c_ξ,n separately.

Consequently, we can apply the framework of projection polynomials to investigate whether the benefits of randomization rely on the precise correlations evident in Eq. (16b) or if randomization of the traditional HIO-algorithm implies significant benefits also in absence of these correlations.

3. Monitoring convergence of the reconstruction

In order to investigate the success and convergence properties of the reconstruction, we monitor three parameters. First of all, we monitor the change of the reconstructed object f⁽ⁱ⁾ from iteration (i − 1) to iteration (i). Second, we measure the mathematical distance to the (sampled) true solution f_in(x⃗). And last, but not least, we measure the deviation in reciprocal space between the a priori given reciprocal input data A_k⃗ and the magnitude of the Fourier transform of the reconstructed object. All three parameters yield different information:

The distance of the reconstructed object f⁽ⁱ⁾ after iteration (i) to the true direct space object f_in is, of course, the best measure for the quality of the reconstructed image. We define this mathematical distance by an angle

φ^{(i)} = arccos (| 〈 f^{(i)}; f_{in} 〉 | / \sqrt{〈 f^{(i)}; f^{(i)} 〉 〈 f_{in}; f_{in} 〉}) .

(17)

The absolute value in the nominator of the argument of the arccos eliminates the influence of the undetermined global phase in f⁽ⁱ⁾. Moreover, φ⁽ⁱ⁾ is identical in direct and reciprocal space (invariance of scalar product upon unitary basis transformations like the Fourier transform).

A dimensionless, normalized error measure ε⁽ⁱ⁾ which depends only on the a priori known (measured) amplitudes A_k⃗ and not on the true direct space solution f_in(x⃗) is

ε^{(i)} = \frac{〈 | {\tilde{g}}^{(i)} | - A; | {\tilde{g}}^{(i)} | - A 〉}{〈 A; A 〉} = \frac{1}{{‖ A ‖}_{2}^{2}} \sum_{\vec{k}} {(| {\tilde{g}}^{(i)} (\vec{k}) | - A (\vec{k}))}^{2},

(18)

where ||·||_p is the p-norm. Since the construction of the absolute value is a non-linear, non-invertible map, the result can in general be quite different from the angle defined above. In fact, two quite different direct space objects f_in (hence quite different complex Fourier transforms) may posses an extremely similar distribution of the magnitude A in reciprocal space [16

,18

J. H. Seldin and J. R. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” J. Opt. Soc. Am. A 7, 412–427 (1990). [CrossRef]

In addition, from a numerical point of view, the change of the reconstructed object f⁽ⁱ⁾ from iteration to iteration is very important, too. This must not be confused with the change of the non-invertible map ε⁽ⁱ⁾ from iteration to iteration: Moving along a connected line of constant error metric does not change the error metric itself, but the reconstructed object f⁽ⁱ⁾ may change arbitrarily. Nevertheless, the algorithm may converge (i.e. no change from iteration to iteration), but not to the true solution f_in of the problem [14

A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984). [CrossRef]

, 16

]. For simulated data, we can judge if the algorithm converged to the true solution f_in by Eq. (17). A suitable choice for observing the change from iteration to iteration is given by the angle

χ^{(i)} = arccos (| 〈 f^{(i - 1)}; f^{(i)} 〉 | / \sqrt{〈 f^{(i - 1)}; f^{(i - 1)} 〉 〈 f^{(i)}; f^{(i)} 〉}) .

(19)

In contrast, the p-norm

δ^{(i)} = {‖ f^{(i - 1)} - f^{(i)} ‖}_{p}

of the difference is not optimal for measuring the convergence of the reconstruction process, because it does not eliminate global phase shifts from iteration (i − 1) to iteration (i).

4. Numerical examples

In this section, we demonstrate the improvements resulting from randomized overrelaxation for three carefully chosen numerical examples: First, we consider two pure phase objects

f (\vec{x}) = exp (i ξ (\vec{x})) Ω_{𝕊} (\vec{x}),

(20)

where ξ(x⃗) is a real function and Ω_𝕊(x⃗) is the shape function

Ω_{𝕊} (\vec{x}) = {\begin{matrix} 0 & if \vec{x} \notin 𝕊, \\ 1 & if \vec{x} \in 𝕊 . \end{matrix}

(21)

The reconstruction of the first phase object (see Fig. 2(a)) is plagued by (at least) one strong local minimum far away from the true solution f_in which results in severe problems for the traditional HIO+ER-algorithm. The reconstruction of the second phase object (see Fig. 2(b)) is plagued by many local minima near the perfect solution f_in, in particular, by stripes [23

J. Fienup and C. Wackermann, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3 (1986). [CrossRef]

]. Both phase objects have a smooth variation of their phase, i.e., phase variations of 2π extend over several sampling points. As a third test object, we investigate a purely real object f_in(x⃗) with strong variation of its magnitude over very short length scales and a different shape (see Fig. 3(a)). In fact, the success rate of the traditional HIO-algorithm without randomized overrelaxation was lowest for this third example (see Fig. 3(b) and Fig. 4). From the numerical results, we conclude that randomized overrelaxation is a powerful extension for a wide class of objects and not only for objects possessing some particular features.

Fig. 2 Pure phase objects used for the investigation of the convergence of the (HIO+O_R)+ER-algorithm. The magnitude of both objects is constant. The phase is plotted using a HSV color-bar, however, the region outside the shape 𝕊 has been set to black. The oversampling ratio is σ = 8.456.

Fig. 3 Investigation of the (HIO+O_R)+ER-algorithm for a purely real object f_in with strong variation of its magnitude over short length scales (oversampling ratio σ = 5.06). In Fig. (b), continuous lines represent the (HIO+O_R)+ER-algorithm, isolated dots the HIO+ER-algorithm. A pure HIO+O_R-calculation without ER is included as black, dash-dotted curve.

Fig. 4 Comparison of the success rate of reconstructions of pure phase objects (see Eq. (20), (22) and (23)) with the HIO+ER- and the (HIO+O_R)+ER-algorithm. The parameter β was fixed to 0.85. Continuous lines represent results of the (HIO+O_R)+ER-algorithm, isolated dots of the HIO+ER-algorithm. A pure HIO+O_R-calculation without ER is included as black, dash-dotted curve.

The first pure phase objects which we used in our numerical simulations is given by

ξ_{1} (x, y) = {(2 π)}^{2} [{(\frac{x}{b_{1}})}^{2} + {(\frac{y}{c_{1}})}^{2}],

(22)

where ξ₁ is defined according to Eq. (20) and the parameters were chosen as b₁ = 1.5515 and c₁ = −1.835. The shape 𝕊 has been restricted to the triangular region defined by (0, 0), (0, a) and (a, 0) with a = 0.97. The object has been sampled on an equidistant grid with N₁ × N₂ = 256 × 256 data points in the interval (x, y) = (−1, −1) to (x, y) = (1, 1). The resulting object is depicted in Fig. 2(a).

The second object is defined by

ξ_{2} (x, y) = (2 π) [{(x - b_{2})}^{3} + {(y - c_{2})}^{2} + \frac{x^{2} y^{3}}{c_{2}^{2} b_{2}^{3}}]

(23)

with the same shape 𝕊 and sampling parameters as in the first example. The parameters were set to b₂ = −0.3515 and c₂ = 0.535. The phase field in this case is depicted in Fig. 2(b).

A reconstruction is considered successful once the angle φ⁽ⁱ⁾ to the input object (see Eq. (17)) falls below some given limit φ_Converged (=0.05° unless stated otherwise). We repeat the reconstruction process for N_Real initial trials (random phases at each k⃗-point, amplitudes

A_{\vec{k}}^{(0)}

equal to given amplitudes A_k⃗). Finally, we calculate the success rate which tells us the percentage of reconstruction processes that have been successful up to iteration (i). A good reconstruction algorithm

should reach a success rate of almost 100%,
should not depend on the starting point (as long as no good starting point is available),
should perform the reconstructions with little computational effort and
should possess these properties for a wide range of its internal parameters (i.e., β, ν, N_HIO and N_ER in case of the (HIO+O_R)+ER-algorithm).

Each iteration in the HIO+ER-algorithm and its (HIO+O_R)+ER-extension can be implemented very efficiently exploiting the good scaling properties of the FFT-algorithm. Therefore, as long as the number of iterations is not too high, the third requirement is clearly met by both, the HIO+ER and (HIO+O_R)+ER-algorithm.

In order to investigate the other points, we plotted the respective success rates for the HIO+ER- and (HIO+O_R)+ER-algorithm for different internal parameters in Fig. 4. We normalized the computational effort, i.e., one iteration of the (HIO+O_R)+ER-algorithm with N_HIO = 130 and N_ER = 10 corresponds to two successive iterations with N_HIO = 50 and N_ER = 20. This prefactor of 2 can be found in the legend of the graph in front of the square brackets for each calculation.

The inclusion of overrelaxation improves the success rate for the first phase object defined by Eq. (22) for each combination of N_HIO and N_ER: in contrast to the HIO based calculations, it quickly converged to 100%. In fact, overrelaxation is so powerful that stagnation in the HIO-algorithm is eliminated entirely and intermediate ER-iterations are no longer necessary (see black dash-dotted curve in Fig. 4). The success rate of the traditional HIO+ER-algorithm stagnates on some level far below 100%. This stagnation is significant in the long term behavior, as we illustrate in Fig. 5 by plotting the success rate of the traditional HIO+ER-algorithm for 50 times as many iterations as in Fig. 4(a). Note, that all trials that did not converge after 2500 iterative steps using the traditional HIO+ER-algorithm with N_HIO = 130 and N_ER = 10 in Fig. 5 are separated from the true solution by an angle greater than φ = 55°.

Fig. 5 Long-term stagnation of the success rate of the traditional HIO-algorithm (without overrelaxation and without randomization) for the first phase object (see Eq. 22) for different choices of the internal parameters (N_HIO, N_ER), but fixed β = 0.85.

If we consider the second phase object defined by Eq. (23), the traditional HIO+ER-algorithm performs better here than for the first phase object. For the one specific set of parameters N_HIO = 130 and N_ER = 10, the success rate of the traditional HIO+ER-algorithm is even slightly superior to the (HIO+O_R)+ER-extension, but for this specific set of parameters, the success rate converges to 100% very quickly for both algorithms. However, for the set of parameters (N_HIO = 50, N_ER = 20) and (N_HIO = 40, N_ER = 30), the success rate quickly converges to 100% only for the (HIO+O_R)+ER-algorithm, but not for the HIO+ER-algorithm. Hence, the (HIO+O_R)+ER-algorithm is much more stable with respect to the particular choice of the internal parameters N_HIO and N_ER including, in particular, the case N_ER = 0.

The traditional HIO+ER-algorithm also has severe problems performing a successful reconstruction for the real object depicted in Fig. 3(a): The success rate for the reconstruction process is very poor if we apply our criterion φ_Converged ≤ 0.05° to determine successful reconstructions. The (HIO+O_R)+ER-algorithm is much more successful and reaches success rates of 100% with few iterations (see Fig. 3(b)). Furthermore, reaching the success rate of 100% does not sensitively depend on the internal parameters N_HIO and N_ER. Again, keeping N_ER small compared to N_HIO or eliminating ER-iterations completely yields the best results. Keep in mind, that neither the reality constraint nor the positivity constraint of the object have been exploited during reconstruction.

For the second and third test object, all trials that failed to converge to the true solution f_in basically evolved towards the correct direction. However, at some point, stagnation in form of stripes close to the true solution f_in appeared. These stripes are very persistent and prevent the usual HIO+ER-algorithm from further convergence to the true solution f_in. Typically, the mathematical distance φ of the objects containing stripes to the true solution f_in is in the order of 0.05 to 5.0 degrees. In case of the purely real test object, after 150 iterations all reconstructed objects reached a distance φ ≤ 0.5° (= 10 · φ_Converged) for HIO₁₃₀ − ER₁₀ and 31% for HIO₅₀ − ER₂₀. The remaining 71% of the objects for HIO₅₀ − ER₂₀ reached a distance φ ≤ 5° (= 100 · φ_Converged) using traditional HIO+ER-algorithm and 150 iterations. Note, that the change χ⁽ⁱ⁾ (see Eq. (19)) from iteration to iteration does not vanish, i.e. numerically, the algorithm did not converge within the given number of iterations.

Figure 6(a) shows the behavior of the success rate of the (HIO+O_R)+ER-algorithm for different values of the parameter ν. In addition, Fig. 6(b) depicts the success rate for different values of the parameter β. Within the range β ∈ [0.5, 1.0], which is typically used in the traditional HIO+ER-algorithm, the (HIO+O_R)+ER-extension does not sensitively depend on the particular choice of β. Therefore, the (HIO+O_R)+ER-algorithm does not depend sensitively on any of its internal parameters N_HIO, N_ER, β and ν. Hence, it fulfills all four requirements on a good reconstruction algorithm which we stated above.

Fig. 6 Investigation of the sensitivity of the (HIO+O_R)+ER-algorithm with N_HIO = 50 and N_ER = 20 on the choice of the parameter β and on the range of the uniform distribution determining the relaxation parameter λ_A for the first phase object (see Eq. 22).

Two concepts (overrelaxation and randomization) have been exploited to modify the HIO-algorithm: Fig. 7 depicts the success rate, if overrelaxation is performed with a static parameter λ_A. We see extremely sensitive behavior of the success rate on the precise value of λ_A if it is kept fixed. Whereas increasing λ_A up to 1.02 improved the success rate for the first object, it completely vanished for λ_A = 1.01 for the second phase object. For a deviation from 1.00 greater than ten percent, almost no successful reconstructions have been observed for any of the three test objects. In addition, we do not a priori know the (non-universal) optimum value for λ_A which depends on the object f_in. Consequently, static overrelaxation does not fulfill the requirements for a good reconstruction algorithm stated above.

Fig. 7 Success rate of the HIO+ER-algorithm for fixed overrelaxation λ_A (no randomization). Parameters are chosen as (N_HIO = 130, N_ER = 10) and β = 0.85.

Based on the discussion in sec. 2.3, we test randomization of the HIO-algorithm for β = 0.85 without employing overrelaxation by setting the amplitudes

c_{ξ, n}^{(R)}

for the random contribution to the coefficients c_ξ,n in the operator Ĥ_{HIO_R} defined in Eq. (14) equal to

c_{𝕊, 1}^{(R)} = c_{A, 1}^{(R)} = c_{𝕊, 2}^{(R)} = c_{A, 2}^{(R)} = R = 0.2

. The deterministic contribution

c_{ξ, n}^{(D)}

is set according to Eq. (15). The result for all three test objects is illustrated in Fig. 8 and compared to the result for the HIO+O_R-algorithm. Clearly, pure randomization without overrelaxation performs equally well as the HIO+O_R-algorithm for all three test objects. However, in our opinion, the HIO+O_R-algorithm should be preferred over randomization without overrelaxation, because (i) it requires only two Fourier transformations for each iteration (lower computational effort) and (ii) it has less degrees of freedom (one uniform random distributions instead of four). Note, however, that due to the large number of degrees of freedom of projection polynomials, further optimization of this approach is possible (including more advanced statistical correlations in the coefficients of the projection polynomial). In addition, more elaborate random distributions might improve both approaches further.

Fig. 8 Comparison of the success rate of both frameworks which provide a generalization of the traditional HIO-algorithm based on randomization. Continuous lines illustrate the behavior for randomized overrelaxation of P_A (see Eq. (10)), whereas dots represent the behavior of the success rate resulting from independent randomization of the coefficients in a projection polynomial (see Eq. (13) and Eq. (14)). In both cases, the deterministic contribution is equivalent to the traditional HIO-algorithm with parameters N_HIO = 140 and β = 0.85. No ER has been performed.

In summary, we succeeded in overcoming stagnation of the traditional HIO+ER-algorithm for various objects. For this purpose, the concept of overrelaxation was applied to the reciprocal space projection in the HIO-algorithm which produced additional degrees of freedom. Randomization of these additional degrees of freedom from iteration to iteration was exploited to prevent stagnation in local minima, especially in traps or tunnels. In particular, this was demonstrated for two pure phase objects with smooth phase variation (with different behavior with respect to the traditional HIO+ER-algorithm) and a purely real object with strong amplitude variations over short length scales. The improved algorithm is far less sensitive to the specific choice of the internal parameters N_HIO and N_ER than the traditional HIO+ER-algorithm. Moreover, the (HIO+O_R)+ER-algorithm remains almost insensitive to the particular value of the internal parameter β in the range β ∈ [0.5, 1.0]. Furthermore, the the good scaling properties in terms of computational time and memory consumption of the HIO+ER-algorithm are fully maintained in the (HIO+O_R)+ER-algorithm. In conclusion, a reconstruction based on the (HIO+O_R)+ER-algorithm exhibits significantly enhanced stability and success probability in comparison with the traditional and commonly used HIO+ER-algorithm.

References and links

1.	F. v. d. Veen and F. Pfeiffer, “Coherent x-ray scattering,” Phys J..: Condens. Matter 16, 5003–5030 (2004). [CrossRef]
2.	R. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Am. A 7, 394–411 (1990). [CrossRef]
3.	U. Pietsch, V. Holy, and T. Baumbach, High-Resolution X-ray Scattering From Thin Films to Lateral Nanostructures (Springer, New York, 2004).
4.	G. N. Afanasiev, Vavilov-Cherenkov and Synchrotron Radiation: Foundations and Applications (Springer, Netherlands, 2004).
5.	I. Robinson and R. Harder, “Coherent x-ray diffraction imaging of strain at the nanoscale,” Nat Mater 8, 291–298 (2009). [CrossRef] [PubMed]
6.	J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature 400, 342–344 (1999). [CrossRef]
7.	M. A. Pfeifer, G. J. Williams, I. A. Vartanyants, R. Harder, and I. K. Robinson, “Three-dimensional mapping of a deformation field inside a nanocrystal,” Nature 442, 63–66 (2006). [CrossRef] [PubMed]
8.	A. Biermanns, A. Davydok, H. Paetzelt, A. Diaz, V. Gottschalch, T. H. Metzger, and U. Pietsch, “Individual GaAs nanorods imaged by coherent x-ray diffraction,” J. Synchrotron Radiat. 16, 796–802 (2009). [CrossRef] [PubMed]
9.	A. A. Minkevich, M. Gailhanou, J.-S. Micha, B. Charlet, V. Chamard, and O. Thomas, “Inversion of the diffraction pattern from an inhomogeneously strained crystal using an iterative algorithm,” Phys. Rev. B 76, 104106 (2007). [CrossRef]
10.	A. A. Minkevich, E. Fohtung, T. Slobodskyy, M. Riotte, D. Grigoriev, T. Metzger, A. C. Irvine, V. Novák, V. Holý, and T. Baumbach, “Strain field in (Ga,Mn)As/GaAs periodic wires revealed by coherent x-ray diffraction,” Europhys. Lett. 94, 66001 (2011). [CrossRef]
11.	A. A. Minkevich, E. Fohtung, T. Slobodskyy, M. Riotte, D. Grigoriev, M. Schmidbauer, A. C. Irvine, V. Novák, V. Holý, and T. Baumbach, “Selective coherent x-ray diffractive imaging of displacement fields in (Ga,Mn)As/GaAs periodic wires,” Phys. Rev. B 84, 054113 (2011). [CrossRef]
12.	J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]
13.	J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A 4, 118–123 (1986). [CrossRef]
14.	A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984). [CrossRef]
15.	D. C. Youla and H. Webb, “Image restoration by the method of convex projections: part 1 – theory,” IEEE Trans. Med. Imaging 1, 81 –94 (1982). [CrossRef] [PubMed]
16.	A. A. Minkevich, T. Baumbach, M. Gailhanou, and O. Thomas, “Applicability of an iterative inversion algorithm to the diffraction patterns from inhomogeneously strained crystals,” Phys. Rev. B 78, 174110 (2008). [CrossRef]
17.	R. P. Millane, “Multidimensional phase problems,” J. Opt. Soc. Am. A 13, 725–734 (1996). [CrossRef]
18.	J. H. Seldin and J. R. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” J. Opt. Soc. Am. A 7, 412–427 (1990). [CrossRef]
19.	R. H. T. Bates, “Fourier phase problems are uniquely solvable in more than one dimension,” Optik (Stuttgart) 61, 247–262 (1982).
20.	L. Auslander and F. A. Grunbaum, “The Fourier transform and the discrete Fourier transform,” Inverse Prob. 5, 149 (1989). [CrossRef]
21.	V. Elser and R. P. Millane, “Reconstruction of an object from its symmetry-averaged diffraction pattern,” Acta Crystallogr., Sect. A: Found. Crystallogr. 64, 273–279 (2008). [CrossRef]
22.	J. Miao, D. Sayre, and H. N. Chapman, “Phase retrieval from the magnitude of the Fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A 15, 1662–1669 (1998). [CrossRef]
23.	J. Fienup and C. Wackermann, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3 (1986). [CrossRef]
24.	H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Phase retrieval, error reduction algorithm, and fienup variants: a view from convex optimization,” J. Opt. Soc. Am. A 19, 1334–1345 (2002). [CrossRef]
25.	H. Takajo, T. Takahashi, R. Ueda, and M. Taninaka, “Study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A 15, 2849–2861 (1998). [CrossRef]
26.	H. Takajo, T. Takahashi, and T. Shizuma, “Further study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A 16, 2163–2168 (1999). [CrossRef]
27.	V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A 20, 40–55 (2003). [CrossRef]
28.	S. Marchesini, “Invited article: A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007). [CrossRef] [PubMed]
29.	S. Marchesini, “Phase retrieval and saddle-point optimization,” J. Opt. Soc. Am. A 24, 3289–3296 (2007). [CrossRef]

OCIS Codes
(100.3190) Image processing : Inverse problems
(100.5070) Image processing : Phase retrieval
(290.3200) Scattering : Inverse scattering
(100.3200) Image processing : Inverse scattering
(110.3200) Imaging systems : Inverse scattering

ToC Category:
Image Processing

History
Original Manuscript: June 18, 2012
Manuscript Accepted: June 25, 2012
Published: July 11, 2012

Citation
Martin Köhl, A. A. Minkevich, and Tilo Baumbach, "Improved success rate and stability for phase retrieval by including randomized overrelaxation in the hybrid input output algorithm," Opt. Express 20, 17093-17106 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-15-17093

Sort: Author | Year | Journal | Reset

References

F. v. d. Veen and F. Pfeiffer, “Coherent x-ray scattering,” Phys J..: Condens. Matter16, 5003–5030 (2004). [CrossRef]
R. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Am. A7, 394–411 (1990). [CrossRef]
U. Pietsch, V. Holy, and T. Baumbach, High-Resolution X-ray Scattering From Thin Films to Lateral Nanostructures (Springer, New York, 2004).
G. N. Afanasiev, Vavilov-Cherenkov and Synchrotron Radiation: Foundations and Applications (Springer, Netherlands, 2004).
I. Robinson and R. Harder, “Coherent x-ray diffraction imaging of strain at the nanoscale,” Nat Mater8, 291–298 (2009). [CrossRef] [PubMed]
J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature400, 342–344 (1999). [CrossRef]
M. A. Pfeifer, G. J. Williams, I. A. Vartanyants, R. Harder, and I. K. Robinson, “Three-dimensional mapping of a deformation field inside a nanocrystal,” Nature442, 63–66 (2006). [CrossRef] [PubMed]
A. Biermanns, A. Davydok, H. Paetzelt, A. Diaz, V. Gottschalch, T. H. Metzger, and U. Pietsch, “Individual GaAs nanorods imaged by coherent x-ray diffraction,” J. Synchrotron Radiat.16, 796–802 (2009). [CrossRef] [PubMed]
A. A. Minkevich, M. Gailhanou, J.-S. Micha, B. Charlet, V. Chamard, and O. Thomas, “Inversion of the diffraction pattern from an inhomogeneously strained crystal using an iterative algorithm,” Phys. Rev. B76, 104106 (2007). [CrossRef]
A. A. Minkevich, E. Fohtung, T. Slobodskyy, M. Riotte, D. Grigoriev, T. Metzger, A. C. Irvine, V. Novák, V. Holý, and T. Baumbach, “Strain field in (Ga,Mn)As/GaAs periodic wires revealed by coherent x-ray diffraction,” Europhys. Lett.94, 66001 (2011). [CrossRef]
A. A. Minkevich, E. Fohtung, T. Slobodskyy, M. Riotte, D. Grigoriev, M. Schmidbauer, A. C. Irvine, V. Novák, V. Holý, and T. Baumbach, “Selective coherent x-ray diffractive imaging of displacement fields in (Ga,Mn)As/GaAs periodic wires,” Phys. Rev. B84, 054113 (2011). [CrossRef]
J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt.21, 2758–2769 (1982). [CrossRef] [PubMed]
J. Fienup, “Reconstruction of a complex-valued object from the modulus of its Fourier transform using a support constraint,” J. Opt. Soc. Am. A4, 118–123 (1986). [CrossRef]
A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A1, 932–943 (1984). [CrossRef]
D. C. Youla and H. Webb, “Image restoration by the method of convex projections: part 1 – theory,” IEEE Trans. Med. Imaging1, 81 –94 (1982). [CrossRef] [PubMed]
A. A. Minkevich, T. Baumbach, M. Gailhanou, and O. Thomas, “Applicability of an iterative inversion algorithm to the diffraction patterns from inhomogeneously strained crystals,” Phys. Rev. B78, 174110 (2008). [CrossRef]
R. P. Millane, “Multidimensional phase problems,” J. Opt. Soc. Am. A13, 725–734 (1996). [CrossRef]
J. H. Seldin and J. R. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” J. Opt. Soc. Am. A7, 412–427 (1990). [CrossRef]
R. H. T. Bates, “Fourier phase problems are uniquely solvable in more than one dimension,” Optik (Stuttgart)61, 247–262 (1982).
L. Auslander and F. A. Grunbaum, “The Fourier transform and the discrete Fourier transform,” Inverse Prob.5, 149 (1989). [CrossRef]
V. Elser and R. P. Millane, “Reconstruction of an object from its symmetry-averaged diffraction pattern,” Acta Crystallogr., Sect. A: Found. Crystallogr.64, 273–279 (2008). [CrossRef]
J. Miao, D. Sayre, and H. N. Chapman, “Phase retrieval from the magnitude of the Fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A15, 1662–1669 (1998). [CrossRef]
J. Fienup and C. Wackermann, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A3 (1986). [CrossRef]
H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Phase retrieval, error reduction algorithm, and fienup variants: a view from convex optimization,” J. Opt. Soc. Am. A19, 1334–1345 (2002). [CrossRef]
H. Takajo, T. Takahashi, R. Ueda, and M. Taninaka, “Study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A15, 2849–2861 (1998). [CrossRef]
H. Takajo, T. Takahashi, and T. Shizuma, “Further study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A16, 2163–2168 (1999). [CrossRef]
V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A20, 40–55 (2003). [CrossRef]
S. Marchesini, “Invited article: A unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum.78, 011301 (2007). [CrossRef] [PubMed]
S. Marchesini, “Phase retrieval and saddle-point optimization,” J. Opt. Soc. Am. A24, 3289–3296 (2007). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper

Figures


Fig. 1	Fig. 2	Fig. 3


Fig. 4	Fig. 5	Fig. 6


Fig. 7	Fig. 8

« Previous Article | Next Article »

OSA is a member of CrossRef.

OpticsInfoBase

Optics Express

Improved success rate and stability for phase retrieval by including randomized overrelaxation in the hybrid input output algorithm

Article Tools

Article Navigation

Abstract

1. The phase retrieval problem

2. Reconstruction algorithms

2.1. Hybrid input output and error reduction

2.2. Hybrid input output with randomized overrelaxation ((HIO+O_R)+ER)

2.3. Projection polynomials

3. Monitoring convergence of the reconstruction

4. Numerical examples

References and links

References

Acta Crystallogr., Sect. A: Found. Crystallogr.

Appl. Opt.

Europhys. Lett.

IEEE Trans. Med. Imaging

Inverse Prob.

J. Opt. Soc. Am. A

J. Synchrotron Radiat.

Nat Mater

Nature

Optik (Stuttgart)

Phys J..: Condens. Matter

Phys. Rev. B

Rev. Sci. Instrum.

Other

Cited By

Figures

Select a Conference
Session or Paper #
Year

OpticsInfoBase

Optics Express

Improved success rate and stability for phase retrieval by including randomized overrelaxation in the hybrid input output algorithm

Article Tools

Article Navigation

Abstract

1. The phase retrieval problem

2. Reconstruction algorithms

2.1. Hybrid input output and error reduction

2.2. Hybrid input output with randomized overrelaxation ((HIO+OR)+ER)

2.3. Projection polynomials

3. Monitoring convergence of the reconstruction

4. Numerical examples

References and links

References

Acta Crystallogr., Sect. A: Found. Crystallogr.

Appl. Opt.

Europhys. Lett.

IEEE Trans. Med. Imaging

Inverse Prob.

J. Opt. Soc. Am. A

J. Synchrotron Radiat.

Nat Mater

Nature

Optik (Stuttgart)

Phys J..: Condens. Matter

Phys. Rev. B

Rev. Sci. Instrum.

Other

Cited By

Figures

2.2. Hybrid input output with randomized overrelaxation ((HIO+O_R)+ER)