Shock-wave tolerant phase reconstruction algorithm for Shack–Hartmann wavefront sensor data

Abstract. We develop a phase reconstruction algorithm for the Shack–Hartmann wavefront sensor (SHWFS) that is tolerant to phase discontinuities, such as the ones imposed by shock waves. In practice, this algorithm identifies SHWFS locations where the resultant tilt information is affected by the shock and improves the tilt information in these locations using the local SHWFS observation-plane irradiance patterns. The algorithm was shown to work well over the range of conditions tested with both simulated and experimental data. In turn, the reconstruction algorithm will enable robust wavefront sensing in transonic, supersonic, and hypersonic environments.


Introduction
The Shack-Hartmann wavefront sensor (SHWFS) is a common tool used to measure optical-path difference (OPD) for a variety of applications.The sensor consists of an array of lenslet subapertures that focus portions of the incident light onto a camera detector.Phase tilts in the subaperture's pupil plane manifest as irradiance patterns shifted from their on-axis locations in the observation plane.These shifts, along with the focal length of the lenslets, are then used to estimate the corresponding tilt across each sub-aperture pupil.2][3] For the purposes of this paper, all findings will be reported in terms of phase, ϕ, which can be related to OPD by ϕ ¼ ð−2π∕λÞ Ã OPD.
5][6][7][8] Shock waves, such as ones that form due to locally supersonic flow, can be modeled as near-instantaneous increases in density which, consequently, result in near-instantaneous changes in the optical phase for light propagating through the shocks.
Recent studies on the effects of aperture-intersecting shock waves on SHWFS data revealed that least-squares reconstruction yields accurate estimates of the change in phase across a shock, jΔϕj, when the magnitude of the phase discontinuity is between 0 and approximately 0.5π rad. 9,10When the magnitude of the phase discontinuity across the shock was greater than 0.5π rad, the shock-induced tilt across the measurement sub-apertures did not result in a proportional shift of the observation-plane irradiance pattern.Consequently, this led to an incorrect estimation of the OPD.This discrepancy was caused by the bifurcation of the irradiance pattern peaks in the observation plane resulting from the shock-induced phase discontinuity in the pupil plane.In other words, this bifurcation causes discrepancies between the pupil-plane tilt and the so-called centroid tilt (C-tilt) measured from the irradiance pattern shift in the observation plane.This behavior has also been observed experimentally. 11,12n the past, SHWFS data have been collected through shock waves; 13 however, the incorrect estimation of Δϕ resulting from the sub-aperture irradiance pattern bifurcation was not addressed.The purpose of this paper is to develop a first-of-its-kind, shock-tolerant phase reconstruction algorithm for the SHWFS utilizing findings presented in Ref. 10.This algorithm shows tremendous potential in extending the usefulness of the SHWFS into optical-turbulence environments where shock waves are present.In practice, the algorithm uses the following steps: 1. Calculate the tilts for each sub-aperture in the traditional sense using SHWFS data (calculate local x and y tilts from observation-plane irradiance pattern deflections).2. Identify sub-aperture locations affected by the shock using a beam-spread metric (identify corrupted tilt information to be replaced).3. Calculate the angle of the shock wave relative to the SHWFS x and y axes.4. Use the observation-plane irradiance pattern to obtain a better estimate of the local tilt across each identified sub-aperture.5. Replace the corrupted tilt information with the new tilt information found in step 4. The tilt data are then reconstructed into a continuous phase sheet using a least-squares reconstructor.
Each of these steps are discussed in Sec. 2. Simulated and experimental data are then used to demonstrate the algorithm, the results of which are presented in Sec. 3. Section 4 summarizes the paper.

Approach
In this section, each step of the algorithm is described in earnest.

Traditional Tilt Calculation
The first step of this algorithm, as outlined above, is to calculate the pupil-plane tilts from the SHWFS data in the traditional sense.Here, the observation-plane irradiance pattern shifts for each sub-aperture are calculated from their on-axis position and used along with the focal length of lenslets to estimate the local pupil-plane tilts.Without any extensions, these tilts are then used in a least-squares reconstructor to estimate continuous ϕ.Hereafter, this approach will be referred to as the "traditional" approach to which the "tolerant" approach, discussed in the next few sections, will be compared.

Identify Affected Sub-Apertures
As shown in Ref. 10, when the jΔϕj caused by a shock wave intersecting the measurement aperture exceeds ∼0.5π rad, the shock-induced phase discontinuity in the pupil-plane of the SHWFS sub-aperture lenslets causes appreciable beam spreading in the resulting observation-plane irradiance patterns.Therefore, in order to identify which sub-aperture lenslets are affected by the shock, beam spread can be calculated for each of the sub-aperture observation-plane irradiance patterns.In this case, the second-moment beam width, otherwise referred to as D4σ, has proven promising.The equation for D4σ along the x axis, D4σ x , is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 7 ; 5 7 0 where n is the total number of pixels inside the sub-aperture irradiance image, I i is the observation-plane irradiance value associated with each pixel, i, inside the sub-aperture, and x is the centroid pixel location in the x direction.Here, the sub-aperture irradiance image has been converted to a raster-ordered vector to simplify the computation of D4σ x . 14A similar equation can be written for D4σ y .It is often convenient to report an overall D4σ, which can be accomplished by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 7 ; 4 7 2 Equation ( 2) is calculated for each sub-aperture observation-plane irradiance pattern and the results are normalized by the diffraction-limited spot size, D DL .For square apertures, D DL is given as D DL ¼ 2λf∕D, and for circular apertures, D DL is given as D DL ¼ 2.44λf∕D, where λ is the wavelength of light, f is the focal length of the lenslets, and D is the diameter of the lenslets.With appropriate thresholding, affected sub-apertures can be identified-i.e., high beam spread identifies a shock-affected sub-aperture, while low beam spread does not.
To illustrate this last point, a phase discontinuity of Δϕ ¼ −π rad was applied to a simulated, uniform, complex-optical field.The phase discontinuity was at a 30 deg angle (measured counterclockwise) from the vertical and offset from the middle of the measurement aperture by half a sub-aperture diameter.This complex-optical field was then applied to an SHWFS model consisting of 16 × 16 sub-apertures.The resultant SHWFS imagery is shown in the leftmost plot in Fig. 1.We see that along the line where the simulated shock occurs, the associated SHWFSirradiance patterns bifurcate.The degree of beam spreading calculated by Eq. ( 2) is shown in the middle plot of Fig. 1.Irradiance patterns that exhibit sufficient beam spread to qualify for being affected by the shock can be identified by thresholding the beam-spread values presented in the middle plot of Fig. 1.The flagged sub-apertures that met the thresholding criteria are presented in the rightmost plot of Fig. 1.In the sections that follow, we discuss the processing steps required to replace the corrupted tilt information in pupil locations flagged by the beam-spread calculation.

Calculate Shock Angle
After the sub-apertures affected by the shock have been identified, the angle that the shock intersects the SHWFS measurement aperture is calculated.This paper explores two methods to achieve this: 1.The first approach fits a line to the locations where the affected sub-apertures were identified.As such, this method assumes that only one shock intersects the SHWFS measurement aperture.Using a MATLAB function referred to as peaks2, 15 the two largest peaks from the observation-plane irradiance pattern were recorded.The accuracy of this method is limited to the number of sub-apertures present in the SHWFS arrangement.2. The second approach employs an image processing method known as a Hough transformation.A Hough transformation can be used to detect otherwise unnoticeable shapes and patterns (i.e., lines) in a noisy image. 16In this approach, a Hough transform is used to find the lines corresponding to the irradiance bifurcation in each sub-aperture.It is known that a line consists of co-linear points ðx i ; y i Þ and can be parameterized by the distance of a line perpendicular to the collection of points, ρ, and the angle, θ, at which ρ is defined, when considering the usual parameterization of a line, ρ ¼ x cosðθÞ þ y sinðθÞ.If θ is confined to an interval, say [0, 90 deg], then each parameterization is unique.However, if a line is parameterized by x i and y i instead of ρ and θ, the line can be identified as a set of co-linear points ðρ i ; θ i Þ, in the parameter space.In other words, a point in the image space corresponds to a line in the parameter space, and a point in the parameter space corresponds to a line in the image space.For the shock angle calculation, individual images of each affected sub-aperture irradiance pattern are first processed to isolate the irradiance peaks.Next, the points are transformed into the parameter space where they are represented as a group of lines.The point in the parameter space where the greatest number of lines intersect defines the ρ and θ of the hidden line in the image space.Since the θ determined through the Hough transform is defined as the angle of the line perpendicular to the line of interest and the irradiance peak bifurcation occurs orthogonal to the shock, the angle extracted from the transform is the angle of the shock.This method calculates the angle of the shock on a sub-aperture by sub-aperture basis and its accuracy is limited by the number of pixels corresponding to each sub-aperture.This approach allows for shock-tolerant phase reconstruction even if multiple shocks are present across the measurement pupil.

Improve Tilt Estimation Using Irradiance Patterns in Affected
Sub-Apertures This section discusses how the resulting irradiance patterns of the sub-apertures affected by the shock can be used to improve the local tilt information.In Ref. 10, it had been shown that the degree of irradiance pattern bifurcation is related to the pupil-plane jΔϕj.To investigate this, phase discontinuities of varying strengths ranging from −π to 0 rad in increments of 0.05π rad were simulated in the middle of a pupil and a Fresnel propagator was used to map the simulated pupil-plane phase to observation-plane irradiance patterns.For each resultant observation-plane irradiance pattern, the amplitude of the largest peak, I 1 , was divided by the amplitude of the second largest peak, I 2 .These results are shown in Fig. 2 as black circle markers where the x-axis is Δϕ∕π and the y-axis is I 1 ∕I 2 .
Fig. 2 Relationship between bifurcated irradiance pattern and phase discontinuity which caused it for square sub-apertures.
When there was no phase discontinuity in the pupil plane, I 1 ∕I 2 was the diffraction-limited peak ratio for a square aperture, 21.19. 17For this case, a one-dimensional slice of the twodimensional irradiance pattern is presented in the top-right plot of Fig. 2. Conversely, as jΔϕj increases, I 1 ∕I 2 decreases.The middle-right plot of Fig. 2 shows the one-dimensional irradiance pattern slice when a phase discontinuity of jΔϕj ¼ 0.67π rad was imposed across the center of the pupil.Here, we see that I 2 began to increase in amplitude relative to I 1 .When jΔϕj ¼ π rad, the phase upstream versus downstream of the discontinuity was perfectly out of phase, and a bifurcated irradiance pattern with equal amplitude peaks formed, i.e., I 1 ¼ I 2 .This last point is shown in the bottom-right plot of Fig. 2.
Conveniently, it has been shown that for square apertures, I 1 ∕I 2 can be related to the phase discontinuity in the pupil plane that caused the bifurcation through the following: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 7 ; 6 0 4 I 1 ∕I 2 ≈ 21.19 expðjΔϕjÞ: (3) In this equation, the constant is the diffraction-limited I 1 ∕I 2 ratio for square apertures. 17Using a fitting algorithm, an equivalent expression was found to work well for circular apertures, given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 7 ; 5 4 4 It is worth noting that in this expression, the constant is not the diffraction-limited I 1 ∕I 2 ratio for circular apertures.Although these equations were arrived at heuristically, they show good agreement with the numerical results obtained using the Fresnel propagator.The results using Eq. ( 3) are plotted as a green line in Fig. 2. For different phase discontinuity and SHWFS conditions, such as angled shocks and experimentally collected data, a different method for relating irradiance patterns to jΔϕj will need to be used.However, for the purposes of this paper and for introducing the theoretical construct behind a shock-wave tolerant phase reconstructor, this analytic development will suffice, and alternative approaches continue to be explored.
It should be noted that Eqs. ( 3) and ( 4) only discern the magnitude of the phase discontinuity in the pupil-plane between 0 and π, not whether Δϕ is positive or negative.As shown in Ref. 10, the same I 1 ∕I 2 ratio can be accomplished when the secondary peak is on either side of the primary peak.However, it was also shown that the location of the secondary peak in relation to the primary peak is directly related to the direction of the flow for the case of an intersecting shock, as well as whether jΔϕj is between 0 and π or π and 2π (c.f.figure 7 in Ref. 10).For the case of a shock, the density increases across the shock and in the direction of freestream flow causing the phase of the light downstream of the shock to lag relative to the phase of the light upstream of the shock.Therefore, by knowing the direction of the flow in relation to the SHWFS measurement aperture, the true Δϕ can be determined.For example, when the flow is moving from left to right across the measurement aperture and the secondary peak in the observation-plane irradiance pattern is to the right of the primary peak, the Δϕ calculated from either Eq.(3) or Eq. ( 4) is correct.However, if the secondary peak is to the left of the primary peak, then the true Δϕ is 2π − Δϕ.These behaviors are flipped if the flow direction is flipped.
After using Eq.(3) or Eq. ( 4) in conjunction with the direction of the flow to determine the phase discontinuity across the shock for the affected sub-apertures, this information can be used to replace the tilt estimates in each sub-aperture.This is accomplished in two steps.First, the calculated change in phase across the discontinuity is converted to a change in OPD using the relation, ΔOPD ¼ ð−λ∕2πÞ Ã Δϕ, which is then used to calculate the overall tilt, θ, caused by and orthogonal to the shock using E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 7 ; 1 8 4 Next, the x and y components of the corrected tilt are calculated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 7 ; 1 3 4 and DeFoor et al.: Shock-wave tolerant phase reconstruction algorithm. . .
where θx and θy are the global tilts and α is the shock wave angle in degrees.Here, the shock wave angle is calculated such that a horizontal shock is considered α ¼ 90 deg and a vertical shock is considered α ¼ 0 deg.θx and θy are calculated by averaging all x and y local tilts from the traditional approach, respectively.

Replace and Reconstruct
After calculation of the improved x and y tilts for sub-apertures affected by an intersecting shock from Eqs. ( 6) and ( 7), this new tilt information is then used to replace the corrupted x and y tilt values found using the traditional approach.It is important to note that the tilt estimates for the unaffected sub-apertures remain unchanged.After this replacement is made, a least-squares reconstruction is applied.

Results
The following sections present the results of employing the algorithm described above on both simulated as well as experimentally collected data.

Simulated Data
In this section, three simple scenarios were simulated for both square and circular sub-aperture geometries to demonstrate the shock-tolerant approach outlined above.For the results that follow, a phase discontinuity was generated and varied from 0 to 2π rad in increments of 0.05π.For simplicity, it was assumed that the SHWFS sub-apertures were illuminated by unit-amplitude, collimated light.Furthermore, it was also assumed that the shock intersects the measurement aperture at α ¼ 0 deg (perpendicular to the x axis) and is located in the middle of a column of sub-apertures.All scenarios outlined in this section follow the procedure described in Ref. 18.The first scenario divided the measurement aperture into a 16 × 16 grid of sub-apertures each with diameter D ¼ 6.25 cm.The second scenario used 32 × 32 sub-apertures each with diameter D ¼ 3.125 cm.The third scenario used 64 × 64 sub-apertures each with diameter D ¼ 1.5625 cm.The entire grid size for each scenario was maintained at 2048 × 2048 grid points.For each case, a thin-lens transmittance function was then applied to each sub-aperture, and angularspectrum propagation was used to obtain an irradiance pattern at focus.These conditions were simulated for both square and circular SHWFS lenslet shapes.After generating simulated SHWFS frames for various Δϕ, multiple thresholding steps were taken to ensure minimal contamination from ambient light and diffraction inside the SHWFS sub-aperture locations.Next, both traditional and shock-tolerant approaches were used to generate x and y tilt matrices.The resulting tilts calculated using both approaches were then used in a least-squares reconstructor, in this case using the Southwell geometry, 3 to obtain a continuous OPD.The reconstructed OPD was then used to estimate the continuous ϕ from which Δϕ was calculated by subtracting the phase values on either side of the shock (Δϕ ¼ ϕ DOWNSTREAM − ϕ UPSTREAM ) for both the traditional and modified approaches.
Figure 3 presents the Δϕ results of the traditional and shock-tolerant approaches plotted as blue circle and red square markers, respectively.Here, the x axes are the input Δϕ∕π of the created phase discontinuities and the y axes are the output Δϕ∕π resulting from reconstructing the tilt values obtained using the traditional and modified approaches.Also plotted are black, dashed lines that represent the expected Δϕ∕π.The top row of Fig. 3 presents the results obtained using square sub-apertures and the bottom row presents the results obtained using circular subapertures.Each column presents results obtained using a different number of sub-apertures.Specifically, the 16 × 16, 32 × 32, and 64 × 64 sub-aperture results are presented in columns from left to right, respectively.
As mentioned earlier and described in Ref. 10, when the jΔϕj induced by the shock wave becomes greater than 0.5π, the traditional phase reconstruction begins to underestimate the true Δϕ across the shock.This is evident for all cases presented in Fig. 3.When the shock-tolerant reconstruction algorithm was employed, we see that across all cases, the reconstructed Δϕ is in DeFoor et al.: Shock-wave tolerant phase reconstruction algorithm. . .better agreement with the true Δϕ.Furthermore, the results presented in Fig. 3 also reveal that as the number of sub-apertures increase, both the shock-tolerant and traditional least-squares reconstruction algorithms deliver better estimations of Δϕ, in their own regard.With more sub-apertures, there is more local tilt information available and least-square fitting error is minimized.
The sinusoidal behavior of the traditional method, shown in blue circles in Fig. 3, is a direct result of the observation-plane irradiance pattern bifurcation.When this bifurcation increases, the irradiance pattern centroid deflection and the corresponding pupil-plane tilt are no longer in agreement for shock-affected sub-apertures.Rather, the energy associated with the pupil-plane phase discontinuity leads to beam spreading in the observation plane.For the case of a shockwave-induced phase discontinuity, the beam spreading is predictable, as shown in Sec.2.4.When the modified approach described above was employed, which leverages the predictability of the shock-induced observation-plane irradiance pattern beam spread, the reconstructed Δϕ∕π was in fair agreement with the input Δϕ∕π until Δϕ∕π ≈ 1.8, depending on the sub-aperture geometry.Above Δϕ∕π ≈ 1.8, the D4σ thresholding algorithm cannot discern which sub-aperture irradiance patterns are affected by the shock due to low beam spread.Therefore, the x and y tilt values calculated using the traditional approach are used instead.

Experimental Data
This section demonstrates the shock-tolerant phase reconstruction algorithm on experimentally collected SHWFS imagery.In order to acquire the experimental data, a 532 nm laser source was first linearly polarized using a λ∕4 wave-plate, expanded, and collimated.The beam was then reflected off a spatial-light modulator (SLM).The SLM was used to impart phase discontinuities ranging from 0 to 2π rad in increments of 0.05π rad at a 0 deg angle in the middle of the beam.In combination with an iris and collimating lens after the SLM, a focus and tilt term were also applied to the SLM in order to decouple the zeroth-order diffractive term from the first-order diffractive term which contains the applied phase information.After which, the beam was re-imaged onto an SHWFS using a 4f relay consisting of two 300 mm lenses.A schematic of this experimental setup is presented in Fig. 4.
The SHWFS used 18 × 18 square, sub-apertures with a diameter of D ¼ 0.148 mm and a focal length of f ¼ 6.7 mm.The detector screen size was 360 × 360 pixels, so each lenslet was comprised of 20 × 20 pixels.The beam-spread threshold used for this data was D4σ∕D DL ¼ 2.9. Figure 5 shows the resultant SHWFS imagery for phase discontinuities of 0, π∕2, and π as well as the sub-apertures that qualify, which are highlighted by red squares.
Due to realistic experimental factors such as non-uniform pupil illumination and camera measurement noise, Eq. ( 3) could not be used to estimate the shock-induced Δϕ from the ratio of bifurcated irradiance pattern peaks.Instead, another equation was developed using a fitting algorithm to empirically collected data.After which, the experimental SHWFS measurements were processed using the traditional and shock-tolerant reconstruction approaches.These results are presented in Fig. 6.
Here, it can be seen that the traditional reconstruction results exhibit the same behavior observed in Fig. 3 as well as in Ref. 10. Furthermore, we see that the shock-tolerant reconstruction results presented in Fig. 6 are in fair agreement with the results presented in the top-right plot of Fig. 3, where the number of sub-apertures was comparable to the number of sub-apertures used in the experiment.We expect that if the number of sub-apertures increased, the Δϕ measured after employing the shock-tolerant reconstruction algorithm would likely be closer to the expected Δϕ.It is also worth noting that for Δϕ∕π values close to 1.35, the shock-tolerant reconstructed Δϕ∕π begins to level off.This is a result of the Gaussian beam used for the experiment.At Δϕ∕π values where the beam spread is low, the algorithm cannot discern between affected and unaffected sub-apertures based on the observation-plane irradiance pattern.This is especially true for those lenslets that lie close to the edge of the beam.The non-uniform illumination caused by the Gaussian beam profile leads to fewer detections of sub-apertures affected by the shock and consequently, more fitting error in the least-squares reconstruction.Nonetheless, the Fig. 4 Experimental setup for the SHWFS data collection (not drawn to scale).