## Abstract

We propose an optical correlation algorithm illustrating a new general method for reconstructing the phase skeleton of complex optical fields from the measured two-dimensional intensity distribution. The core of the algorithm consists in locating the saddle points of the intensity distribution and connecting such points into nets by the lines of intensity gradient that are closely associated with the equi-phase lines of the field. This algorithm provides a new partial solution to the inverse problem in optics commonly referred to as the phase problem.

© 2014 Optical Society of America

## 1.Introduction

Solving the phase problem in optics has attracted much attention, primarily in problems of diagnostics of an object structure within microscopy, pattern recognition, terrestrial telescopy, and biomedical optics [1, 2]. In general, the phase problem consists in deriving the spatial phase distribution for complex fields, including speckle fields from a measured intensity distribution. In general, the inverse problem has no universal solution. The-state-of-the-art in this area of investigations is reflected in a series of original publications and reviews [3–15]. (See specifically paper [16] containing a comprehensive list of references.) A novel promising approach in solving the phase problem arises from the concept of singular optics [17]. To be precise, the following propositions are assumed as a basis for this approach [18]: (*i*) amplitude zeroes (also named optical vortices, or wave front dislocations, or phase singularities) are the ‘reference’, structure-forming elements, whose set constitutes a singular skeleton of a field; (*ii*) spatially distributed amplitude zeroes obey the specific sign principle governing the characteristics (signs) of adjacent zeroes; (*iii*) the spatial distributions of intensity and phase in complex fields are interconnected. Therefore, knowing the locations and signs of amplitude zeroes, one can predict, at least in a qualitative manner, the behavior of a field (including spatial phase distribution, with accuracy not exceeding $\pi /2$) in all regions between the amplitude zeroes. Within this approach, solving the phase problem for scalar, *viz.* homogeneously polarized coherent optical fields, is reduced to (*i*) the development of reliable and practicable algorithms for location of amplitude of an intensity distribution with randomly located zeroes, and *(ii*) searching for the physically most attractive algorithm for reconstruction of the spatial phase distribution (spatial phase map) of a field. Routinely, the phase problem of this kind is solved efficiently by imposing a coherent reference wave onto the field to be analyzed [19–21]. But the use of a reference wave is impractical or even impossible in many cases, especially for distant diagnostics or in microscopy. Besides, the feasibility for obtaining an integrated phase map of a field based on the amplitude zeroes alone has not found substantial justification. At the same time, reconstruction of a phase skeleton of a complex optical field by finding the ‘reference’, structure-forming elements whose positions and characteristics provide reliable prediction of the behavior of the field parameters at all other areas is important for one more reason. As a matter of fact, reconstruction of the phase skeleton from a measured intensity distribution can provide vast possibilities for data compression within problems associated with optical telecommunications involving complex optical fields. That is why looking for new algorithms for obtaining at least partial solutions of the phase problem hinges not only to fundamental problem of modern optics, but can have important impact in applications of optical technologies.

The aim of this paper is to substantiate - proceeding from simple intuitive suppositions - the algorithm for reconstruction of the spatial phase distribution from a measured spatial intensity distribution of a complex (speckle) field. As a result, in contrast to the standard singular optical approach dealing with amplitude zeroes, we consider just the saddle points of a spatial intensity distribution as the structure-forming elements, as these points are interconnected by the lines of intensity gradient and really constitute (together with these lines) the phase skeleton of a field. In our opinion, the proposed algorithm might be considered as a partial case of a more general new method for solving the inverse problem in optics. We will represent the initial results of computer simulation within the scalar approximation, *viz.* assuming a homogeneously polarized field.

The proposed algorithm includes the following three actions: (*i*) bicubic spline interpolation [22]; (*ii*) location of the saddle points of intensity; (*iii*) connecting the saddle points of intensity by the gradient lines.

## 2. Location of the saddle points of intensity

A saddle point is the point from the function domain that is stationary not being a local extremum. Derivative of the function equals zero at this point in the transverse directions. That is why, primary analysis consists in choosing all points where the two derivatives simultaneously equal zero. As maximas and minimas of the function obey this condition, one must determine the criterion for selecting only the saddle points. The algorithm implemented by us is illustrated in Fig. 1. When passing a saddle point of intensity (blue point) one goes by turns through the points of alternate minima and maxima shown in Fig. 1 by green and red points, respectively. Therefore, the magnitudes of maxima (minima) are always larger (less) than the magnitude of a function at the saddle point. Thus, one initially identifies the stationary points of a field, where *dI/dx =* 0 and *dI/dy =* 0. Then, if passing a stationary point one meets alternating minima and maxima (with magnitudes larger and smaller than at the specified stationary point), this point is identified as a saddle point.

## 3. Connecting saddle points by intensity gradient lines

As a rule, *viz.* with probability 95%-98% as it will be argued later, see Fig. 4, the saddle points of intensity are located within the regions of rapid changing phase [18]. Therefore, the regions with small intensity gradients (smooth spatial changes of intensity) are the regions with rapid change of phase. In average, the modulo of the gradient of phase at the saddle points of intensity exceeds by $\sqrt{2}$ times its average magnitude [20, 21]. In turn, the gradient lines of intensity (whose reconstruction is described below) going from the saddle points correspond to equi-phase lines with the same probability. For this reason, *the saddle points of intensity* are chosen as ‘structure-forming’, steady-state points from which *the gradient lines of intensity* are reconstructed facilitating *phase mapping* of complex, spatially inhomogeneous optical fields.

From a saddle point of intensity, S (see Fig. 2), one draws the (dashed) line connecting the points with the maximal numerical gradient directed towards the specified saddle point. At the phase map, this line corresponds to the area of the smoothest spatial change of the phase. Passing the saddle point of intensity, one always meets two minima, which means that each saddle point is the origin of two gradient lines (lines 3-2-1-S and 3*-2*-1*-S). These lines approach the spatial area with increasing intensity. As a rule, such gradient lines for the intensity in complex fields intersect and eventually form a spatial net.

Figure 3 illustrates implementation of the described algorithm for a random speckle field. Fragments Figs. 3(a) and 3(b) characterize the intensity distribution and the phase distribution, respectively. The spatial phase distribution is represented as areas with magnitudes of phase within the limits: 0 to π/2, π/2 to π, π to 3π/2, and 3π/2 to 2π, associated with areas with different gray graduations, from white to dark-gray. The saddle points and the maxima of intensity are depicted by light-blue triangles and rhombuses, respectively; the phase singularities (amplitude zeroes) of opposite signs are depicted by red and dark-blue. The intensity gradient lines reconstructed following the described algorithm are shown in yellow. Zero lines for real and imaginary parts of the field’s complex amplitude are shown in red and dark-blue, respectively.

Figure 3 illustrates two peculiarities of the reconstructed gradient lines: (*i*) non-intersecting lines passing the saddle point connect phase singularities of opposite signs; (*ii*) the form of most of the gradient lines approximately reproduces the boundaries of areas of changing phase within the intervals: 0 to π/2, π/2 to π, π to 3π/2, and 3π/2 to 2π. (Note that the choice of the mentioned intervals of changing phase is, to a certain extent, conventional. This digitization of phase is the closest one to the Rayleigh’s criterion in classical optics. In some cases, depending on the required accuracy, this criterion can be replaced by a stronger one.)

These conclusions have been proven by the following analysis: We simulate a speckle field with specified phase distribution (within the far-field diffraction approximation). Reconstructing a pair of intensity gradient lines (as an analogue to the equi-phase lines for a field) originating from a saddle point, we compute the phase difference at the saddle point and at each point of the gradient line. Then, for the set of found magnitudes of phase difference, we establish histograms and estimate the confidence interval for deviation of a phase magnitude at each point of the gradient line from the phase at the initial saddle point. Analysis has been performed for four different speckle field distributions. The histograms for the corresponding phase distributions are shown in Fig. 4.

For versatile proving the proposed approach, we have performed simulation for several groups of objects including both random and quasi-deterministic sets of point coherent sources. Quasi-deterministic objects are modeled as sets of point sources randomly distributed within a small area but spatially repeated, so that the total distribution of sources takes up regularity, viz. periodicity. Figures 4(a)–4(c) correspond to speckle fields from a quasi-deterministic object whose linear size equals 1/20 of the shown speckle field. Therefore, the objects consist of 3 and 2 groups of sources for Figs. 4(a) and 4(c), respectively. Figures 4(b) and 4(d) correspond to speckle field from random set of point sources whose linear size equals 1/20, Fig. 4(b), and 1/10, Fig. 4(d), of the shown speckle field.

Analysis of the represented dependences shows that the width of the confidence interval does not exceed 4°, i.e. the reconstructed gradient lines of intensity are really the equi-phase lines of a field, within the accepted accuracy.

To estimate the calculation time of implementation of the proposed algorithm here realized in Delphi software, note that this time mainly depends on the matrix size. There are estimations for processor Intel5, 2.6 GH, for intensity distribution for a set of 20 speckles: 200$\times $200 matrix – 0.4 sec, 400$\times $400 matrix – 2.1 sec ( + 425%), and 800$\times $800 matrix – 14.3 sec ( + 3725%). For comparison, for intensity distribution for the set of 60 speckles: 200$\times $200 matrix – 0.5 sec (plus 25% for 20 speckles), 400$\times $400 matrix – 2.4 sec ( + 380%) or plus 14% for 20 speckles, and 800$\times $800 matrix – 15.4 sec ( + 2980%) or plus 1% for 20 speckles.

Let us shortly discuss the sources of errors in this simulation. At first, both in recording of the speckle pattern using a CCD camera (as in this study) and in attempts to compensate for the effect of atmospheric turbulence in terrestrial telescopy by means of adaptive optics [6–8], one meets the serious (and universal) problem of digitization of ‘rough’ (intrinsically analogue) spatial distributions of intensity and phase. As a matter of fact, the edges and corners of cells in both cases inevitably produce additional phase singularities, whose number may be commensurable or even exceed the number of ‘true’ singularities intrinsic to the field *per se*. This circumstance can drastically distort the processed data and the result of processing. Thus, one must search for proper procedures for data smoothing, based either on bicubic spline interpolation [22] used by us, or based on Gaussian smoothing [23]. This circumstance serves for estimating the boundaries of applicability of the proposed algorithm. So, discreteness of the sensitive area of a detector and the necessity for data smoothing (spatial averaging) presume the structural elements of an analyzed field, viz. speckles be much larger than CCD camera pixels, at least by two orders of magnitude to provide correct data processing. In the described simulation, average linear speckle size exceeds the linear size of pixels by a factor of approximately 500. Secondly, available means for registration of complex intensity distributions (such as CCD cameras) are characterized by strong nonlinearities, especially within regions of vanishing intensity (amplitude zeroes). That is why, at least now, data proceedings based on CCD registration cannot be considered reliable enough. Correct compensation of recording nonlinearity presumes development of adequate data processing algorithms going beyond the framework of this paper.

## 4. Conclusions

We have introduced an optical correlation algorithm for reconstructing the phase skeleton of a complex optical field from a measured two-dimensional intensity distribution. The algorithm consists in location of the saddle points for the intensity distribution and subsequently connecting these points into nets by the lines of intensity gradient, which are closely related to the equi-phase lines of the field. It has been demonstrated that the set of saddle points of intensity and the intensity gradient lines facilitates the reconstruction of the phase skeleton of a complex scalar (homogeneously polarized) coherent optical field. The proposed algorithm provides a new partial solution to the phase problem in optics. Significance of this result follows from the possible use of complex two-dimensional optical fields within optical telecommunications. Actually, the proposed algorithm can be implemented for any polarization projection selected at the output of a detector even for an inhomogeneous polarized field. Moreover, finding out the skeleton of a speckle field consisting of a set of saddle points and connecting the lines of intensity gradient, being a rather time- and labor-consuming operation of processing of the measured field intensity distribution, nevertheless results in essential data compression undergoing transmission. So, even in the reported initial computer experiment such compression reaches a factor of 14.9, approximately. Namely, from 600x600 pixels of a frame, skeleton is reproduced by 13706 (3.8%), 22920 (6.4%), 30868 (8.6%), and 29290 (8.1%) pixels for frames of Figs. 4(a)–4(d), respectively. Note, the mentioned data compression is unprecedented for 2D signals of high complexity. Therefore, the behavior of a field at all other areas is reliably predicted (within the above mentioned accuracy), which provides saving coding with the option of high-quality recovering of the transmitted data. Consequently, the introduced algorithm might be considered as a partial case of a more general new method for solving the inverse problem in optics.

Generalization of the proposed approach with experimental verification will soon be reported elsewhere.

## Acknowledgments

This material is based upon work supported by the Ministry of Education and Science of Ukraine (grant No 0113U003242) and was partially supported by the Ministry of Education and Science of the Russian Federation (projects No 8703 and No 8877) and Joint Project SB RAS and NAN Ukraine (2013-2014). Steen G. Hanson acknowledges the financial support from the Danish Council for Technology and Innovation under the Innovation Consortium LICQOP, grant #2416669.

## References and links

**1. **R. H. T. Bates and M. J. McDonnell, *Image Restoration and Reconstruction* (Caledon Oxford, 1986).

**2. **T. Acharya and A. K. Ray, *Image Processing – Principles and Applications* (Wiley InterScience, 2006).

**3. **E. Abramochkin and V. Volostnikov, “Two-dimensional phase problem: differential approach,” Opt. Commun. **74**(3–4), 139–143 (1989). [CrossRef]

**4. **E. Abramochkin and V. Volostnikov, “Relationship between two-dimensional intensity and phase in a Fresnel diffraction zone,” Opt. Commun. **74**(3–4), 144–148 (1989). [CrossRef]

**5. **V. Volostnikov, “Phase problem in optics,” J. Sov. Laser Research **11**(6), 601–626 (1990). [CrossRef]

**6. **M. Loktev and V. Volostnikov, “Singular wavefields and phase retrieval problem,” Proc. SPIE **3487**, 141–147 (1998). [CrossRef]

**7. **K. G. Larkin, D. J. Bone, and M. A. Oldfield, “Natural demodulation of two-dimensional fringe patterns. I. General background of the spiral phase quadrature transform,” J. Opt. Soc. Am. A **18**(8), 1862–1870 (2001). [CrossRef] [PubMed]

**8. **K. G. Larkin, “Natural demodulation of two-dimensional fringe patterns. II. Stationary phase analysis of the spiral phase quadrature transforms,” J. Opt. Soc. Am. A **18**(8), 1871–1881 (2001). [CrossRef]

**9. **F. Yu. Kanev, V. P. Lukin, and L. N. Lavrinova, “Correction of turbulent distortions based on the phase conjugation in the presence of phase dislocations in a reference beam,” Atmos. Oceanic Opt. **14**, 1132–1169 (2001).

**10. **V. P. Lukin and B. V. Fortes, “Phase-correction of turbulent distortions of an optical wave propagating under conditions of strong intensity fluctuations,” Appl. Opt. **41**(27), 5616–5624 (2002). [CrossRef] [PubMed]

**11. **V. A. Tartakovsky, V. A. Sennikov, P. A. Konyaev, and V. P. Lukin, “Wave reversal under strong scintillation conditions and sequential phasing in adaptive optics,” Atmos. Oceanic Opt. **15**, 1104–1113 (2002).

**12. **J. R. Fienup, “Lensless coherent imaging by phase retrieval with an illumination pattern constraint,” Opt. Express **14**(2), 498–508 (2006). [CrossRef] [PubMed]

**13. **M. Wielgus and K. Patorski, “Evaluation of amplitude encoded fringe patterns using the bidimensional empirical mode decomposition and the 2D Hilbert transform generalizations,” Appl. Opt. **50**(28), 5513–5523 (2011). [CrossRef] [PubMed]

**14. **M. Trusiak, K. Patorski, and M. Wielgus, “Adaptive enhancement of optical fringe patterns by selective reconstruction using FABEMD algorithm and Hilbert spiral transform,” Opt. Express **20**(21), 23463–23479 (2012). [CrossRef] [PubMed]

**15. **D. Barchiesi, “Numerical retrieval of thin aluminium layer properties from SPR experimental data,” Opt. Express **20**(8), 9064–9078 (2012). [CrossRef] [PubMed]

**16. **J. R. Fienup, “Phase retrieval algorithms: a personal tour [Invited],” Appl. Opt. **52**(1), 45–56 (2013). [CrossRef] [PubMed]

**17. **M. S. Soskin and M. V. Vasnetsov, “Singular optics,” Prog. Opt. **42**, 219–276 (2001).

**18. **I. Mokhun, “Introduction to linear singular optics,” in *Optical Correlation Techniques and Applications*, Ed. O. V. Angelsky, (2007), Chap. 1, TA 1630.A6, 1–133.

**19. **N. B. Baranova, A. V. Mamaev, H. F. Pilipetsky, V. V. Shkunov, and B. Y. Zel’dovich, “Wavefront dislocations: topological limitations for adaptive systems with phase conjugation,” J. Opt. Soc. Am. **73**(5), 525–528 (1983). [CrossRef]

**20. **N. Freund, Shvartsman, and V. Freilikher, “Optical dislocation networks in highly random media,” Opt. Commun. **101**(3–4), 247–264 (1993). [CrossRef]

**21. **Y. Galushko and I. Mokhun, “Characteristics of scalar random field and its vortex networks. Recovery of the optical phase,” J. Opt. A: Pure Appl. Opt. **11**094017 (2009).

**22. **R. Keys, “Cubic convolution interpolation for digital image processing,” IEEE Trans. Signal Processing Acoust. Speech Signal Processing **29**(6), 1153–1160 (1981). [CrossRef]

**23. **V. I. Vasil’ev and M. S. Soskin, “Analysis of dynamics of topological peculiarities of varying random vector fields,” Ukr. J. Phys. **52**, 1123–1129 (2007).