Optical Phase Unwrapping in the Presence of Branch Points

Strong turbulence causes phase discontinuities known as branch points in an optical field. These discontinuities complicate the phase unwrapping necessary to apply phase corrections onto a deformable mirror in an adaptive optics (AO) system. This paper proposes a non-optimal but effective and implementable phase unwrapping method for optical fields containing branch points. This method first applies a least-squares (LS) unwrapper to the field which isolates and unwraps the LS component of the field. Four modulo-2π-equivalent non-LS components are created by subtracting the LS component from the original field and then restricting the result to differing ranges. 2π phase jumps known as branch cuts are isolated to the non-LS components and the different non-LS realizations have different branch cut placements. The best placement of branch cuts is determined by finding the non-LS realization with the lowest normalized cut length and adding it to the LS component. The result is an unwrapped field which is modulo-2π-equivalent to the original field while minimizing the effect of phase cuts on system performance. This variable-range ‘φ LS + φnon−LS’ unwrapper, is found to outperform other unwrappers designed to work in the presence of branch points at a reasonable computational burden. The effect of improved unwrapping is demonstrated by comparing the performance of a system using a fixed-range ‘φLS + φnon−LS’ realization unwrapper against the variable-range ‘φLS + φnon−LS’ unwrapper in a closed-loop simulation. For the 0.5 log-amplitude variance turbulence tested, the system Strehl performance is improved by as much as 41.6 percent at points where fixedrange ‘φLS + φnon−LS’ unwrappers result in particularly poor branch cut placement. This significant improvement in previously poorly performing regions is particularly important for systems such as laser communications which require minimum Strehl ratios to operate successfully. © 2008 Optical Society of America OCIS codes: (010.1080) Active or adaptive optics; (010.7350) Wave-front sensing References and links 1. J. W. Hardy, Adaptive Optics for Astronomical Telescopes (Oxford University Press, 1998). 2. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping Theory, Algorithms, and Software (John Wiley & Sons, Inc., 1998). 3. D. L. Fried, “Branch point problem in adaptive optics,” J. Opt. Soc. Am. A 15 (1998). 4. D. L. Fried, “Adaptive optics wave function reconstruction and phase unwrapping when branch points are present,” Optics Communications 200 (2001). 5. J. Barchers, “The performance of wavefront sensors in strong turbulence,” Proc. SPIE 4839 (2003). #92309 $15.00 USD Received 30 Jan 2008; revised 23 Apr 2008; accepted 23 Apr 2008; published 1 May 2008 (C) 2008 OSA 12 May 2008 / Vol. 16, No. 10 / OPTICS EXPRESS 6985 6. D. L. Fried and J. L. Vaughn, “Branch cuts in the phase function,” Appl. Opt. 31 (1992). 7. D. L. Fried, “The Nature of the Branch Point Problem in Adaptive Optics,” Proc. SPIE 3381 (1998). 8. C. A. Primmerman, “Atmospheric-compensation experiments in strong-scintillation conditions,” Appl. Opt. 34 (1995). 9. M. C. Roggemann, “Branch-point reconstruction in laser beam projection through turbulence with finite-degreeof-freedom phase-only wavefront correction,” J. Opt. Soc. Am. A 17 (2000). 10. Wikipedia (2007). Website, URL http://en.wikipedia.org/wiki/Vector-potential . 11. J. Barchers, “Personal Communication,” (2007). 12. T. Rhoadarmer, “Development of a self-referencing interferometer wavefront sensor,” Proc. SPIE 5553 (2004). 13. L. C. Andrews and R. C. Philips, Laser Beam Propagation through Random Media (SPIE Optical Engineering Press, 1998). 14. T. J. Brennan and P. H. Roberts, AOTools The Adaptive Optics Toolbox-Users Guide version 1.2 (Optical Sciences Company, 2006).


Introduction/Background
Adaptive-optics (AO) systems have proven to be very effective at compensating for the effect of weak atmospheric turbulence. [1] Stronger turbulence in which intensity variation known as scintillation occurs is more challenging. [2] Significant scintillation can cause branch points where intensity is zero and phase is undefined. [3] Branch points cause lines of 2π phase discontinuity known as branch cuts in the field leading to the branch points.
These unavoidable branch cuts degrade AO system performance when correcting with a continuous surface deformable mirror due to the inability of the mirror to fit the required discontinuous phase. [4] Cut placement, however, affects the amount of degradation and cuts can be placed between branch points in many different ways. This paper proposes an effective method of unwrapping the phase of an optical field which places branch cuts where their negative impact on system performance is minimized.

Wavefront Sensors
There are many types of wavefront sensors (WFSs), but for strong turbulence conditions and the associated high levels of scintillation, interferometric sensors such as the self-referencing interferometer (SRI) are good choices. [5] Regardless of the source, this paper assumes amplitude and wrapped phase information is available for the field.

Irrotational and Rotational Fields
The phase of an optical field φ Tot (x, y) can be divided into irrotational φ Irr (x, y) and rotational φ Rot (x, y) components. [2] Under weak turbulence conditions, only the irrotational component of the phase is present because the total field is irrotational. [2] In the absence of detector noise, any rotational component is the result of branch points (Sec. 1.3) which are caused by strong turbulence.

Branch Points
Branch points are optical phenomena which occur at places of zero amplitude in an optical field. More specifically, they are defined as places about which a closed integral of phase gradients is non-zero. [2] In mathematical terms, ∇ × F(r) = 0 where F(r) = ∂ φ ∂ x i + ∂ φ ∂ y j is the gradient of the scalar function φ (x, y). Thus starting at a point A and integrating phase differential along a closed path around a branch point back to A yields a different phase than what was started with. Since the phase at point A has not changed, the integration along the closed path will necessarily be an integer multiple of 2π. A clockwise integration around a single positive branch point yields +2π while a negative branch point yields −2π. If multiple branch points are contained within an integration path, their effects sum with positive branch points cancelling out the effects of negative branch points. In principle, integrating phase gradients along a path around a point is the way the presence of branch points is determined. [3] Throughout this paper, 'x' and 'o' are used in figures to indicate the presence of positive and negative branch points, respectively.

Phase Cuts
A phase cut is a line across which there is an integer multiple of 2π phase discontinuities. In a sampled phase map, a phase cut is defined as wherever there is more than a π difference in phase between adjacent samples.
Segmented deformable mirrors (DMs) deal with phase cuts very easily since each pixel in a segmented DM is independent. Modulo-2π compensation is applied and cuts are essentially ignored. The phase is not required to be unwrapped prior to implementation on a segmented DM.
Typically, in the absence of branch cuts, continuous facesheet DMs outperform segmented DMs because they ramp smoothly from actuator to actuator, better matching the phase corrections needed between sample points. Phase cuts, however, degrade system performance because the DM cannot change shape abruptly and instead changes smoothly between actuators in attempting to match a phase cut. Regions between samples on either side of a cut are compensated poorly. As such, it is advantageous to eliminate phase cuts wherever possible and keep them short and through areas of low illumination when they cannot be eliminated.
Throughout this paper, phase cuts are depicted in figures by lines. The line colors are usually white but may vary from figure to figure in an effort to keep them distinct from the background.

Wrapping Cuts
Phase cuts take two forms, wrapping cuts and branch cuts. A wrapping cut is only due to the field being wrapped and proceeds from one edge of the optical field to another. It can be eliminated by adding or subtracting an integer multiple of 2π to the field on one side of the wrapping cut. Cuts which form a closed path within the field are also unwrapping cuts and can be eliminated similarly by either adding or subtracting an integer multiple of 2π to the interior or exterior of the cut path. As an example, Fig. 1 depicts a wrapped and unwrapped phase. Note that the wrapped phase is limited in range to [−π, π) while the unwrapped phase is not.   Figure 2 shows a phase with both wrapping and branch cuts. Unlike wrapping cuts, branch cuts do not extend across the entire field (or in a closed path) having at least one end terminating at a branch point. [6] They either connect branch points of opposite polarity or connect a branch point with the edge of the optical field (in effect placing a branch point of opposite sign just off the field at that point). By terminating at a branch point, they compensate for the non-zero curl of phase differential around the branch point. In a closed path around a single branch point, the phase differentials integrate to ±2π. As the line integral crosses the branch cut, however, ∓2π is added so that the closed line integral sums to zero as it would if there were not a branch point within the closed path. Branch cuts can be placed in a variety of ways, all of which will still compensate for the non-zero curl of branch points in the phase. Two examples of phase cut placement are shown in Fig. 3. The poor unwrap is created by simply unwrapping the field from left to right. The minimum cut distance unwrap was manually created to minimize the length of the branch cuts.

Least-Squares Unwrappers
Least-squares (LS) unwrappers are very common methods of estimating the unwrapped phase of an optical field in AO systems designed for weak atmospheric turbulence. [2] There are two types, weighted and unweighted. For an N × N array of phases, an unweighted LS unwrapper is developed as where G is a 2N(N − 1) × N 2 transformation matrix that converts the N 2 vector of phases φ into a 2N(N − 1) vector of phase differentials in the x and y directions s and the inverse notation is taken to be the pseudo-inverse. In weak turbulence, s is most commonly the phase gradients provided by a Shack-Hartmann WFS. If actual phases φ Tot are available, the phase differentials s are developed as s = W (Gφ Tot ) where W () indicates the wrapping operation of limiting the differentials s to some 2π interval. An important point is that while creating an N 2 × N 2 pseudo inverse is computationally daunting, the problem is alleviated somewhat by G being sparse and fixed for a given AO system. Much of the work can be pre-computed a single time rather than having to be determined in real time during execution.

Weighted LS Unwrappers
Weighted LS unwrappers are sometimes used to minimize noise or emphasize certain parts of a field. In a weighted LS unwrapper, the slopes are weighted before applying the pseudo-inverse as It works essentially the same as an unweighted LS unwrapper, but the pseudo inverse cannot be pre-computed because the weighting matrix is not typically constant. This makes a weighted LS unwrapper difficult to implement in real-time systems.

LS Unwrappers and the Hidden Phase
In estimating the phases from the slopes, there is an implicit assumption that the sum of phase differentials is path independent, or that the field is irrotational. As a result, the phase estimate of an LS unwrapper is irrotational. The LS unwrapper does not reconstruct the rotational portion of the phase, which is why the rotational component of the phase is sometimes referred to as the "hidden phase". [7] This makes a simple LS unwrapper alone a non-optimal choice when compensating for strong turbulence. [8,9]

Non-LS Component of the Field
The non-LS component of the field is the difference between the original field and the output of a LS unwrapper. If the original field is irrotational, the output of the LS unwrapper will be modulo-2π-equivalent to the original field, and the non-LS component will be non-existent. If the original field has branch points and is rotational, the effects of those branch points will be isolated in the non-LS component. As such it is sometimes referred to as the rotational component. [2] Strictly speaking, the rotational component is not unique [10], so it is referred to here as the non-LS component, uniquely identifying it as the difference between the original field and the output of a LS unwrapper. For the purposes of this paper, the non-LS component will be wrapped to a particular 2π range.

Improved Unwrapper
The first step in unwrapping efficiently in the presence of branch points is generating the LS and non-LS components of the field through the use of an LS unwrapper, where LS() indicates applying an LS unwrapper operation to the vector of wrapped phases φ Tot and W () indicates wrapping the phase to some 2π range.  Thus total phase φ Tot adjusted to remove wrapping cuts while still retaining branch cuts can be determined as where U () indicates an unwrapping process which removes wrapping cuts (but not branch cuts). While removing any wrapping cuts, this unwrapped result is modulo-2π-equivalent to  φ Tot , maintaining both the irrotational and rotational components of the field. This has been covered in several texts [2] and is a common way of including rotational phase effects in the AO systems being developed to operate under strong turbulence conditions. [11] In general this approach is reasonably effective, as shown in Fig. 4. Here the field whose phases are in Fig. 2 has its wrapping cuts removed by the process depicted in Eq. (1). The resultant branch cuts are plotted over the intensity (instead of the phase as in previous figures) to show the effectiveness of the unwrap. In this case the branch cuts are reasonably short and seem to avoid the areas of high intensity, although they may not be optimal.
While generally effective, this unwrap method sometimes gives less appealing results as shown in Fig. 5. Here the branch cuts are much longer than they could be and go through areas of high intensity. Admittedly this is the worst realization encountered in the simulation, but poor results are all too often encountered.
Since after unwrapping the LS portion of the phase field φ LS is free of phase cuts, the non-LS portion φ Non−LS must be examined in order to reduce the impact of phase cuts. Being wrapped, φ Non−LS is restricted to some 2π range, say [0, 2π). If the range is changed to [−π/2, 3π/2) then all the points whose phase is in [3π/2, 2π) would have 2π subtracted from them. The resulting field would be modulo-2π equivalent to the original field, but would have branch cuts in different positions. The field depicted in Fig. 5 is re-depicted in Fig. 6 alongside unwraps for the same field with φ Non−LS having differing range restrictions. The unwrap with φ Non−LS restricted to [0, 2π) has terrible branch cut placement. The remaining three realizations depicted are much more reasonable, with the realization created by limiting φ Non−LS to [−π, π) having the lowest normalized cut length. It should be noted that the creation of four realizations is reasonable because the majority of the computational load is in executing the LS unwrapper which only has to be done once.

Unwrapping Metric -Normalized Cut Length
Having developed multiple modulo-2π-equivalent phase realizations, it is necessary to compare different branch cut placements so that the best one can be chosen. Short cuts through regions of minimal illumination have the least impact on system performance. [6] As such, the metric used in this work is 'normalized cut length' which is the line integral of field intensity along any phase cuts divided by the average intensity of the field. It is an indication of what proportion of light in the system is along phase cuts. Since light along branch cuts is erroneously corrected For a discretely sampled field, normalized cut length is determined by first isolating the phase cuts within the field. This can be accomplished by taking the difference between adjacent pixels first up and down and then side to side. The intensities on either side of the cuts are then summed and divided by two to account for the average intensity along the cuts. Finally the result is normalized by dividing by the sum of the field's intensities.
The advantage of normalized cut length is that it can be computed during system execution and is highly correlated to system performance. In order to show the correlation of normalized cut length to system performance, a 256 × 256 complex 'Fine' field is developed from the 32 × 32 field by interpolation from the coarser field. Similarly a 32 × 32 array of phases from an unwrapper is converted into a 256 × 256 arrays of phases. This models an idealized DM whose surface varies smoothly between actuators. The DM model is translated into the complex domain byû + iv =Â exp(iφ ) whereû andv are the real and imaginary estimates of the field, A is the estimated amplitude of the field andφ is the estimated phase of the field. The field- where F is the 'Fine' field, E is the estimated DM field and * is the conjugation operator. [12] By developing both the 'Fine' field and estimated DM field by interpolating from the same coarser field, degradations in the field estimation Strehl are isolated solely to the effect of phase cuts in the field. This allows direct comparison between normalized cut length and the effect of phase cuts on field estimation Strehl.
Normalized cut length is plotted against field estimation Strehl ratio for the various fields and unwrapping methods examined during this work in Fig. 7 and has a correlation of −0.9982. Normalized cut length is shown to be a good measure of the impact of phase cuts on field estimation Strehl ratio, and thus on system performance.

Simulation and Results
In order to test the unwrapper, a simulation was created that isolates and unwraps 32 × 32 sections from a 513 × 513 optical simulation generated test screen. The simulation was run on two test screens depicting fields with intensity log-amplitude variances of 0.4 and 0.8. Each had a scaling of 16 pixels per the atmospheric coherence diameter r 0 . The log-amplitude variance of intensity is a measure of the scintillation of the field and a reasonable indication of the turbulence strength. [ Fig. 5 have an unwrapped version where a cut passes through or close to a region of high intensity. This is certainly an unwrapping solution to avoid and justifies creating an unwrapper which can choose the best of four unwrap realizations.
For each of the 232,324 possible realizations, the integrated cut intensity metric is recorded for the four φ non−LS ranges. The average and maximum score for all realizations is then determined for each of the four ranges. These data show how the unwrapper would perform if the range was fixed to a particular range. The integrated cut intensity is also recorded for the φ non−LS range which gives the lowest score. The average and maximum is determined for this best of four φ non−LS ranges and compared against the average and maximum scores from the fixed ranges.
The results of the unwrapper using both unweighted and weighted (by field intensity) LS unwrappers to separate out the rotational component are given in Tables 1 and 2. Compared to limiting the non-LS component to a single range, the variable-range 'φ LS + φ non−LS ' mean normalized cut length is reduced in both cases. Perhaps more importantly, the worst realizations are avoided in a variable-range 'φ LS + φ non−LS ' unwrapper so that the maximum normalized cut length is dramatically reduced. The weighted variable-range 'φ LS + φ non−LS ' unwrapper has the effect of influencing the LS portion of the field towards the areas of higher intensity. The non-LS portion of the field is then influenced towards the areas of lower intensity and branch cuts are forced into darker portions of the field. While a weighted LS unwrapper has the best performance, the computational cost of a weighted unwrapper is significant (see section 4).

Comparison to Other Unwrappers
In order to evaluate the worth of the variable-range 'φ LS + φ non−LS ' method, it was compared to other unwrappers designed to work with branch points. The other unwrappers are the fixedrange 'φ LS + φ non−LS ' unwrap, Goldstein's branch cut placement unwrap method [2], Waveprop's xphase [14], and Fried's smoothphase. [4] The fixed-range 'φ LS + φ non−LS ' unwrapper is the same as the variable-range 'φ LS + φ non−LS ' but only develops a single unwrap realization instead of choosing the best of four realizations. Goldstein's branch cut placement method attempts to determine minimum length branch cuts that connect branch points. [2] Xphase is a Matlab unwrapping function from the AOTools Matlab toolbox. It is designed to work with fields containing branch points and attempt to place branch cuts in low intensity regions of the field. [14] It should be noted that xphase required the 32 × 32 field to be zero-padded to 64 × 64 in order to work properly. Fried's smoothphase unwrapper separates the field into rotational and irrotational components by first determining the rotational component (after balancing the number of branch points by adding additional branch points along the edge of the field as necessary). Once separated, the irrotational component can be unwrapped and then recombined with the rotational component of the field. [4] The comparison between unwrapping methods is given in Tables 3 and 4 for log-amplitude variances of 0.4 and 0.8, respectively. Execution time is the time needed to execute an unwrapper in Matlab on a Pentium 4 CPU (3.2GHz) with 2.0 GB of RAM over the 230,000+ frames tested. While execution times may depend on Matlab implementation, indications from this simulation are that the variable-range 'φ LS + φ non−LS ' unwrapper using an unweighted LS gives the best performance at a reasonable computation burden. The variable-range 'φ LS + φ non−LS ' unwrapper using a weighted LS improves performance still more, but at an unreasonable computational burden. AOTools xphase unwrapper gave slightly improved results compared to a variable-range 'φ LS + φ non−LS ' using an unweighted LS unwrapper, but at over six times the computational burden.

Impact on System Performance
The purpose of developing an improved unwrapper is to improve the performance of a closedloop AO system encountering strong turbulence. As such, a 1000 frame closed-loop AO simulation was performed under 0.5 log-amplitude variance strong turbulence in order to compare the effect of the unwrapping on system performance.
With the exception of the log-amplitude variance, simulation conditions were purposely benign in order to isolate the unwrapping as the dominate factor on system performance. The remaining simulation conditions were r 0 = 4D SA where D SA is the diameter of a sub-aperture, sample rate = 223 f G where f G is the Greenwood frequency of the atmosphere, average SNR 200. The simulation used a leak-free integrator controlled feedback with a error signal gain of 0.4. The control law was applied immediately before the unwrapper, whose output then went to the DM.    age improvement of the variable-range 'φ LS + φ non−LS ' unwrapper over a fixed-range 'φ LS + φ non−LS ' unwrapper is difficult to determine. In order to develop an average improvement, the simulations were extended to 10,000 frames to provide each fixed-range of the 'φ LS + φ non−LS ' unwrappers with areas of both good and bad performance. The results of the simulation were put into histograms and then summed to form cumulative distribution functions (CDFs) shown in Fig. 9. The CDFs show how the variable range φ non−LS unwrapper improves performance. The CDF of the variable range φ non−LS unwrapper is shifted to the right when compared to the CDF of the fixed range φ non−LS unwrapper. Not only does this indicate improved average performance, but indicates more significant improvement for systems such as laser communication where performance thresholds which inhibit operation below a certain Strehl ratio.

Conclusion
In the presence of branch points, unwrapping the phase is a difficult problem. Isolating the rotational component by using a LS unwrapper to separate the field into its LS and non-LS components seems an excellent approach. The wrapping phase cuts of the irrotational component are automatically eliminated by the LS unwrapper. Altering the range of the rotational component is a simple and effective way of varying the placement of the branch cuts associated with the rotational phase component, and computing the normalized cut length is an effective way of comparing the effectiveness of branch cut placements. Choosing the best of four branch cut realizations not only improves average cut placement but, perhaps more significantly, eliminates the worst cut placements which would significantly degrade AO system performance. The improved unwrapping eliminates regions of degraded performance where previous unwrappers yielded poor branch cut placements. The reduced areas of poor performance not only improves average performance, but may significantly improve systems such as laser communications where falling below a performance threshold will cause signal fading.

Disclaimer
The opinions and views expressed by the authors are not necessarily those of the Department of Defense or the United States Air Force.