Iterative parallel registration of strongly misaligned wavefront segments

: The paper presents an algorithm for the precise registration of multiple wavefront segments containing large misalignment and phase differences. The measurement of a wavefront with huge dynamics or a large aperture size can be carried out in multiple Shack-Hartmann sensor measurements of segments of the wavefront. The registration algorithm is flexible with respect to the shape of the wavefront and can reconstruct a plane as well as divergent wavefronts, making it suitable for freeform wavefronts. The algorithm enables parallel registration of the wavefront segments which is carried out in an iterative manner to compensate for large misalignment errors. A simulative analysis of the proposed algorithm compares its performance to a fast parallel registration (FPR) algorithm and the established iterative closest point (ICP) algorithm. For a sensor misalignment of up to 100 µ m and 3 mrad the algorithm registers a plane and a divergent wavefront with a precision that is a factor 4 and 12 better than the registration precision of the FPR and ICP algorithm.


Introduction
A Shack-Hartmann sensor (SHS) is well suited for the evaluation of optical systems, as it provides a vibration-insensitive and reference free measurement of the optical wavefront with a large dynamic range [1,2]. This makes the SHS a frequently used device in ophthalmology [3], adaptive optics [4,5], free-space optical communication [6], optical system alignment [7] and production of optical systems and components [2]. However, wavefronts exceeding the sensors dynamic range or aperture size can only be captured in one measurement by using additional supporting null optics [8]. The drawback of supporting optics is that they are additional sources of errors, limiting the measurement quality [9]. In an alternative concept the SHS is combined with a positioning system to enable a measurement of the wavefront beyond the dynamic range or aperture size of the sensor [10][11][12]. In particular, segments of the wavefront are measured by the SHS part by part. Adjacent wavefront segments are measured with a spatial overlap, enabling the reconstruction of the entire wavefront using registration algorithms [13,14] that are also capable of handling deviations of the sensor from the intended measurement positions. As registration errors limit the measurement quality of the wavefront, a precise registration algorithm is necessary for the assessment of high-end optical systems. In the last decades freeform optics grew in popularity because of their high optical performance [15]. The shape of wavefronts generated by freeform optics can be of any type, demanding algorithms for the registration of wavefronts beyond plane wavefronts. In [16] a fast and parallel registration (FPR) algorithm is reported, where wavefront segments are registered in parallel within less than a second. The algorithm is analysed with respect to sensor misalignment up to 5 µm and 200 µrad and shows high quality results. However, a poor calibration of the positioning system or application-specific requirements with respect to measurement time, travel range, etc. might entail even larger sensor misalignment. The registration performance of the FPR algorithm deteriorates in case of larger sensor misalignment, as the used local approximation of the global mismatch metric limits its applicability to measurements with only moderate sensor misalignment.
The contribution of this paper is the development and evaluation of an algorithm that enables high-quality registration of plane as well as divergent wavefronts also in the presence of large sensor misalignment of the order of 100 µm and 3 mrad. Section 2. introduces the algorithm and discusses its properties. Section 3. presents a simulative analysis of the algorithm and Section 4. concludes the paper.

Measurement concept
In the proposed measurement concept the SHS measures segments of the entire wavefront at specific sensor positions [11,12]. Typically, the sensor deviates from its nominal position and alignment, because of uncertainties and errors in the positioning system (see Fig. 1), resulting in misaligned wavefront segments in the global frame (FG). The size of a wavefront segment is limited to the size of the sensor aperture. If adjacent measurements overlap, the misaligned wavefront segments can be registered to reconstruct the entire wavefront. In particular, the wavefront segments are registered by minimizing their overlap mismatch using rigid body transformation and wavefront propagation as illustrated in Fig. 2. The latter one is used to minimize phase differences between the wavefront segments caused by the misalignment of the sensor or a scan trajectory deviating from the wavefront of a specific phase. The sensor aperture of an SHS consists of a lenslet array, where at each lenslet the local gradient of the incident wavefront is measured [17]. From the measured gradients the corresponding wavefront segment can be reconstructed in form of a point cloud, which is a set of threedimensional points contained in the wavefront segment. For wavefront segments with large dynamics the corresponding point clouds can be determined from the phase distribution on the lenslet array [16], as the phase gradients are directly determined from the SHS measurement [18]. There are several algorithms available for the reconstruction of a distribution from a discrete set of local gradients, which can be divided into zonal and modal reconstruction algorithms [19,20]. Zonal reconstruction is typically preferred, as it better preserves details of the wavefront [21], which are important for a successful registration of the segments. The normal vector of the wavefront segment at each point of the point cloud is directly determined from the gradient measurements and is necessary for the proposed registration algorithm.

Iterative parallel registration algorithm
After the wavefront reconstruction the point cloud of each wavefront segment i = 1..U is determined in the local coordinate system of the SHS. For the reconstruction of the entire wavefront the sensor position and alignment in FG is necessary for each measurement. However, because of uncertainties in the positioning system, the exact position of the sensor is subject to uncertainty. The initial guess for the sensor position of the measurement of wavefront segment i is defined by the translation vector T 0i ∈ R 3 and the rotation matrix R 0i ∈ R 3×3 . Conveniently, the nominal position of the sensor is used for T 0i and R 0i . FSi is defined as the local coordinate system of the SHS positioned in FG with T 0i and R 0i . The reconstructed point cloud of segment i is directly represented in FSi, denoted as P i 0i with elements x i 0ij ∈ R 3 and normal vectors n i 0ij ∈ R 3 . The upper index defines the coordinate system in which an object is represented and j is an index to specify an individual point. Transformation from FSi into FG is given by the rigid body transformation x 0ij = R 0i x i 0ij + T 0i ∈ P 0i and n 0ij = R 0i n i 0ij , as illustrated in Fig. 3. The upper index is omitted for objects represented in FG. As the actual position deviates from the assumed nominal position, P 0i is not correctly positioned and the point clouds have an overlap mismatch. Additionally, the actual sensor position might be at different phases. To remove the phase differences, the wavefront segments have to be propagated and as an initial guess S 0i ∈ R can be defined for the propagation distance of segment i. To register the wavefront segments in parallel, each point cloud is transformed by where θ i ∈ R 3 defines a rotation with R(θ i ) as the corresponding rotation matrix. k i ∈ R 3 defines a translation and s i ∈ R the distance along which the wavefront segment is propagated additionally to S 0i . The parameters are with respect to FSi and collected in the vector a T i = (k T i , θ T i , s i ) ∈ R 7 . A metric for the overlap mismatch between two transformed segments (i,k) is where the functions W i i (·, a i , S 0k ) = W i i (a i , S 0k ) and W i k (·, a k , S 0i ) = W i k (a k , S 0i ) denote the transformed segments in FSi. q kin ∈ R 2 is a sampling point that belongs to the overlapping region of the segments and lies in the x-y plane of FSi. The squared differences between the segments are added up. Instead ofM ik (a i , a k ) the metric is considered, where both wavefront segments are back-propagated by S 0i providing the advantage that only one segment has to be propagated by presumed propagation data. Despite backpropagation, Eq. (4) can be used as an alternative to Eq. (3) for registration of the segments [14]. With Eq. (2) the following point clouds are contained in the segment functions of Eq. (4): transformed from FSk to FSi given by With Eq. (4) a global mismatch metric for the entire overlap mismatch in the set of wavefront segments is given by which is the sum over all overlapping wavefront segments and N denoting the total number of sampling points. A T = (a T 1 , .., a T U ) ∈ R 7 U are the transformation parameters of all segments. For the global mismatch metric of Eq. (7), the FPR algorithm considers for each overlapping segment pair (i,k) the point clouds of Eq. (5) with initial positioning data, i.e.
The FPR algorithm is carried out in three steps. In Eq. (8) the two point clouds are represented in the local coordinate system of one point cloud, i.e. FSi. In the first step of the FPR algorithm, the point cloud P i 0i of the local coordinate system and the corresponding normal vectors n i 0ij are interpolated for subpixel registration, leading to the interpolants F i (x, y) ∈ R for the point cloud and N i (x, y) ∈ R 3 for the normal vectors. As P i 0i denotes the unpropagated point cloud with respect to FSi, it has the same shape for any presumed registration data, defined by R 0i , T 0i and S 0i . Hence, for any change of the presumed registration data the same interpolants can be used meaning that they have to be determined only once which is the advantage of using the global mismatch metric of Eq. (7). χ i kin ∈ P i k (0, ∆S 0ki ) define the points that belong to the overlapping region with P i 0i . The x-y components of χ i kin define the sampling points q kin ∈ R 2 of the metric for the overlap mismatch in Eq. (4). In the second step of the FPR algorithm the corresponding point in The segment functions at the sampling points in Eq. (4) can then be approximated by where z ikn ,z ikn ∈ R are the z components of χ i kin andχ i kin . C ikn ,C ikn ∈ R 7 are determined by χ i kin andχ i kin and the corresponding normal vectors. With Eq. (10) the global mismatch metric of Eq. (7) can be approximated by with B ikn =z ikn − z ikn . In the third step of the FPR algorithm the parameters A * that register the point clouds are determined by minimizing the approximated global mismatch metric of Eq. (11). The minimization of this expression can be carried out by solving a matrix equation where Q ∈ R V×7 U and B ∈ R V are determined by the coefficients of Eq. (11) with V denoting the total number of squared terms. As registration determines the entire wavefront up to a rigid body transformation and a wavefront propagation, the seven transformation parameters of one point cloud are set to 0 to make Q T Q invertible. Conveniently, this is the point cloud of the wavefront segment in the center of the set of segments to keep registration parameters small. If the parameters A * are large, owing to large sensor misalignment, the quality of the approximation of the wavefront segments at the sampling points in Eq. (10) decreases and the FPR algorithm shows registration errors. The iterative fast parallel registration (IFPR) algorithm is based on the FPR algorithm but not limited by the approximation of Eq. (10).
Similarly to the FPR algorithm, the IFPR algorithm is initialised with T 0i , R 0i and S 0i for the presumed sensor position and phase data of the measurements. The first three steps of the IFPR algorithm are the same as for the FPR algorithm. First, the interpolation of P i 0i and the corresponding normal vectors is carried out leading to F i (x, y) and N i (x, y). Second, for each overlap the corresponding points and normal vectors are determined based on Eq. (9). Third, the parameters that minimize Eq. (11), i.e. A * = (a * T 1 , .., a * T U ), are determined by solving Eq. (12).
Inserting Eq. (2) into Eq. (1) using A * leads to The can be considered as a better guess for the registration data. In the fourth step of the IFPR algorithm the relative change of the global mismatch metric with respect to the improved registration data is considered to evaluated whether the registration is sufficient. The relative change of the global mismatch metric is given by where M 0g and M 1g are M g (A = 0) (see Eq. (7)) evaluated for T 0i , R 0i , S 0i and T 1i , R 1i , S 1i , respectively. The approximation of the global mismatch metric in Eq. (11) is exact for A = 0, hence it can be used to compute Eq. (14). If there is a significant relative change of the global mismatch metric by using the improved registration data, there might be the chance for an additional improvement of the registration data by applying the FPR algorithm with T 1i , R 1i , S 1i as the presumed registration data. This is equivalent to repeating step two and three of the IFPR algorithm.
Step one is not repeated, owing to the fact that the same interpolants are used to approximate Eq. (11) for any presumed registration data. The iteration including step two to step four is then repeated until the relative change of the global mismatch metric is smaller than a threshold denoted by ε indicating the convergence to the correct registration data T * i , R * i and S * i and serving as a stopping condition for the iteration. The steps and the iteration of the IFPR algorithm are illustrated in Fig. 4.
The complexity of the first iteration in dependence of the number of segments (U) and the number of points per overlap (PPO) used for the registration can be divided into the following terms The first term describes the complexity of the interpolation (step 1), which is independent of the number of PPO. The second term describes the complexity of step 2 and 4, where α ∈ R with αU being the number of overlaps. The last two terms describe the complexity of step 3. In particular, they describe the complexity of the symmetric product Q T Q and the Cholesky decomposition used to solve Eq. (12) [22]. In conclusion, Eq. (15) shows that the computational effort of one iteration has a linear dependence on the number of PPO and a cubic dependence on the number of segments.

Algorithm analysis
The performance of a registration algorithm is influenced by several quantities, e.g. sensor misalignment, measurement noise, etc. The performance of the IFPR algorithm is evaluated with respect to these quantities in a simulation-based analysis. For a comparison of the IFPR algorithm to other algorithms, the FPR algorithm [16] and the established iterative closest point (ICP) algorithm [23] are considered.

Simulation setting
A plane and a divergent wavefront, depicted in Fig. 5(a) and 5(b) respectively, is considered for the evaluation of the registration performance of the algorithms. The plane wavefront is generated by collimating a spherical wavefront with a meniscus lens [14,24]. It has a diameter of 50 mm and a peak to valley (PV) of 11 µm containing mainly spherical aberration (PV= 11 µm) and secondary spherical aberration (PV= 0.5 µm). The second wavefront is a spherical wavefront (diameter of 30 mm) with a divergence of 140 • and contains mainly spherical aberration (PV= 6 µm), astigmatism (PV= 7 µm) and coma (PV= 6 µm) resulting in a total PV of 13 µm. The plane wavefront is measured at 25 sensor positions in the x-y plane arranged in a chessboard pattern and the sensor aperture is a square with a side length of 13 mm. There are 43 sensor positions for the measurement of the divergent wavefront, where a circular sensor aperture with a diameter of 7 mm is used. The size of the lenslets, contained by the sensor aperture, is set to 130 × 130 µm 2 meaning that the total number of lenslets per sensor aperture is 10.000 for the plane wavefront and 2.200 for the divergent wavefront. The wavefronts and the SHS measurements are simulated with a custom raytracing software implemented in MATLAB (The MathWorks Inc., Natick, MA, USA) and using OpticStudio (Zemax LLC, Kirkland, WA, USA). With the software sensor misalignment, measurement noise as well as systematic measurement errors are simulated. After the measurement, the wavefront segments are reconstructed by a spline-based zonal reconstruction algorithm [25]. Then the reconstructed wavefront segments are registered by the algorithms. The algorithms are running on a personal computer with 6 cores and a processor frequency of 2.6 GHz.
With the ICP algorithm the wavefront segments are registered sequentially, meaning that, starting from the wavefront segment in the center, the wavefront segments are consecutively added in a spiral way [14]. For the IFPR and FPR algorithm, the point clouds are interpolated with cubic interpolation and the normal vectors with linear interpolation [16]. The threshold for the stopping condition of the IFPR algorithm (see Fig. 4) is set to ε = 1 3 , which turns out to be a suitable value, as during convergence the metric is typically decreased by orders of magnitude with one iteration. The result of an algorithm is evaluated in three steps. First, the registered wavefront is fitted into the original wavefront and the difference between the wavefronts is computed. Second, measurement noise and systematic measurement errors are removed from the difference to determine only the registration errors. Third, the root mean square (RMS) and PV values of the registration errors are determined.

Reference configuration
Sensor misalignment defines a deviation of the actual sensor positioning data from the nominal positioning data and is reflected by k i ∈ R 3 and θ i ∈ R 3 of Eq. (2). In particular, k i reflects translational and θ i rotational misalignment of wavefront segment i. For the simulation of misalignment, the parameters are randomly distributed between predefined misalignment ranges, which is [−50, 50] µm for the components of k i and [−1500, 1500] µrad for the components of θ i . Measurement noise is simulated by overlaying the point clouds of the segments with a normal noise distribution with zero mean and a standard deviation of 10 nm. The overlap area of adjacent measurements is set to 20 % of the area of the sensor aperture. In the measurement of the divergent wavefront the overlap size is 20 % of the area of the sensor aperture or larger because of the complex shape of the wavefront. The average number of PPO used for registration is 928 for the plane wavefront and 315 for the divergent wavefront. In the following sections the algorithms are analysed with respect to misalignment ranges, noise standard deviation, overlap size and the number of PPO, where the values used for these quantities in this section define the reference configuration.
The registration errors of the IFPR, FPR and ICP algorithm for the reference configuration are depicted in Fig. 6 and Fig. 7 for the plane and the divergent wavefront, respectively. The IFPR algorithm attains high quality registration results of both wavefronts with an RMS registration error of 14 nm for the plane and an RMS registration error of 31 nm for the divergent wavefront. Registration of the divergent wavefront leads in general to larger errors, since a larger number of wavefront segments has to be registered and a smaller number of PPO is used as compared to the plane wavefront. The IFPR algorithm requires 3 to 4 iterations to register the wavefronts of the reference case as illustrated in Fig. 8 demonstrating its fast convergence. Depending on the initial guess for the registration data a certain amount of iterations is needed, as the approximation of transformed segments of the FPR algorithm (see Eq. (10)) is less qualitative the larger the necessary transformation parameters are. This explains the larger RMS registration errors of the FPR algorithm of 24 nm for the plane wavefront and 122 nm for the divergent wavefront depicted in Fig. 6(b), 7(b) and 7(c). The IFPR algorithm registers the plane wavefront at least about a factor 2 better than the other algorithms. For the divergent wavefront the improvement is a factor 4. The ICP algorithm has the largest registration errors, as the wavefront segments are sequentially registered, leading to an accumulation of the registration errors. Moreover, the ICP algorithm can not propagate the wavefront segments and does not compensate phase differences between them, leading to enlarged registration errors for the divergent wavefront. Analysis shows that the IFPR algorithm attains comparable results with respect to freeform wavefronts. Increasing the dynamics of the plane wavefront to a PV of 1 mm the RMS registration error is still around 10 − 20 nm.

Influence of misalignment
Sensor misalignment is caused by uncertainties and errors in the positioning system. The misalignment is divided into translational and rotational misalignment reflected by parameters k i and θ i for segment i. In the simulation the components of k i and θ i are randomly picked for a defined misalignment range.
The influence of the misalignment range on the registration quality of the algorithms with respect to the plane and the divergent wavefront is depicted in Fig. 9. Translational misalignment ranges are considered up to ±100 µm. For the rotational misalignment two ranges, i.e. ±1.5 and ±3 mrad, are exemplarily considered. The values for the misalignment ranges are on a realistic order with respect to a multi-axis positioning system [26]. The plots show the RMS registration error on a logarithmic scale in dependence of the misalignment range. The PV registration error is typically larger than the RMS registration error a factor 5 to 8. The IFPR algorithm has for all considered misalignment ranges an RMS registration error smaller 25 nm for the plane wavefront and smaller 50 nm for the divergent wavefront. For large misalignment of ±100 µm and ±3 mrad, the IFPR algorithm has registration errors a factor 4 and a factor 12 smaller than the other algorithms for the plane and the divergent wavefront respectively.
In general, the IFPR algorithm attains better and more robust registration performance than the other algorithms. Only for the divergent wavefront and small misalignment the FPR algorithm has slightly smaller registration errors than the IFPR algorithm. This is explained, as for wavefront segments with weak features, the wavefront segments might be shifted more against each other   than necessary with many iterations. Hence, there are two reasons for the stopping condition of the iteration. First, it avoids unnecessary iterations saving computation time. Second, in case of wavefront segments with weak features, it prevents the IFPR algorithm from too large in-plane shifts of the wavefront segments. For translational misalignment smaller ±30 µm and rotational misalignment smaller ±1.5 mrad, the FPR algorithm obtains results comparable to the IFPR algorithm, as the approximation of Eq. (10) is still sufficiently good.

Influence of noise and systematic error
Background light, readout and dark currents are sources of noise in the measurement with an SHS [27]. The total measurement noise is simulated by adding a normal noise distribution to the wavefront segments. Moreover, an SHS measurement might contain systematic errors. In this section a systematic measurement error is simulated by an additional error distribution added to the wavefront segments with a PV = 5 nm. The distribution of the systematic error is depicted in Fig. 10 and based on a typical systematic error distribution of a calibrated SHS [28]. The influence of measurement noise on the algorithms performance in presence of the systematic measurement error is shown in Fig. 11. For both wavefronts and the considered noise standard deviations the IFPR algorithm attains the best registration results of the algorithms with a minimum RMS registration error of 6 nm for the plane and 23 nm for the divergent wavefront. Compared to the results of the reference configuration where no systematic error is simulated, the IFPR and the FPR algorithm do not decrease in registration performance, while the ICP algorithm, especially for the divergent wavefront, has larger registration errors explained by the reduced robustness of sequential registration to measurement errors as compared to parallel registration. Especially for wavefront segments with weak features and large misalignment, the wavefront segments might be shifted to far by the FPR algorithm, as Eq. (10) qualitatively describes the transformed segments only for sufficiently small transformation parameters. Basically, noise prevents the segments from getting shifted too far, as the overlap mismatch increases. Hence, the registration errors might increase with a smaller noise standard deviation.

Influence of points per overlap
The points in the overlap regions contain the surface information that enables registration of the segments. Typically, more PPO improve the registration result, because more information of the segments shape is available. Nevertheless, a smaller number of PPO has the benefit of a decreased computation time. In Fig. 12 the influence of the average number of PPO used for registration on the computation time of the algorithms and on their registration results is depicted with a uniform distribution of the selected points in the overlap region.
As expected, the computation time of the algorithms decreases with a smaller number of PPO. Only the computation time of the ICP algorithm for the registration of the plane wavefront increases despite a decrease of the number of PPO from 150 to 100, which is explained by more iterations used by the algorithm. The IFPR algorithm has the smallest registration errors for all considered numbers of PPO, demonstrating its robustness despite a small amount of surface information. For 100 PPO the IFPR algorithm has a computation time between 200 ms and 300 ms with an RMS registration error of 21 nm for the plane and 64 nm for the divergent wavefront. Increasing the number of PPO to 300 reduces the RMS registration error for the divergent wavefront to 32 nm while the computation time increases to 470 ms. For the plane wavefront the RMS registration error of the IFPR algorithm increases when the number of PPO is increased from 200 to 300 showing that more PPO might reduce the registration quality in some cases. With 600 PPO, the IFPR algorithm registers the plane wavefront in 410 ms with an RMS registration error of 15 nm. Despite requiring three iterations to register the wavefronts, the computation time of the IFPR algorithm is not three times the computation time of the FPR algorithm. The reason for this is that some parts of the FPR algorithm have to be carried out only in the first iteration of the IFPR algorithm, e.g. the interpolation of the point clouds, the determination of the points belonging to an overlap, and need not to be repeated in following iterations.

Influence of overlap size
Besides increasing the number of PPO, more surface information is obtained by increasing the overlap size. Additionally, the algorithms get more sensitive to out of plane angles between the wavefront segments. The drawback of an increased overlap size is that more wavefront segments have to be measured leading to an increased measurement time. The RMS registration error in dependence of the overlap size, in percent of the area of the sensor aperture, is shown in Table 1. For all overlap sizes the average number of PPO used for registration is set to a constant value of 300. By increasing the overlap size from 20 % to 40 % the RMS registration error of the IFPR algorithm decreases by a factor of 3 to 4. In particular, the RMS registration error gets smaller than 10 nm for both wavefronts making the IFPR algorithm applicable for the evaluation of high-end optical systems. The RMS registration error of the ICP algorithm decreases by a factor of 15 when the overlap size is increased from 20 % to 40 %, demonstrating the high sensitivity of sequential registration to available surface information. The registration quality of the FPR algorithm hardly improves with the considered increase of the overlap size. This is explained by the incorrect interpretation of the surface information, as the approximation of Eq. (10) is less qualitative for the considered large sensor misalignment. In summary the high quality registration performance of the IFPR algorithm with RMS registration errors down to 10 nm in the presence of large sensor misalignment up to 100 µm and 3 mrad is successfully demonstrated. With a computation time of less than 500 ms on a personal computer the evaluation of high-end optical systems in time critical applications is possible.

Conclusions
In this paper, an algorithm for the precise registration of SHS measurements is proposed. The benefit of the proposed algorithm is that it can cope with large sensor misalignment, where other algorithms lack in performance. The wavefront to be registered, can be a plane as well as a divergent wavefront, making the algorithm applicable for a wide range of tasks including the evaluation of freeform optics. The algorithm registers the wavefronts in an iterative manner and mathematics are discussed. A simulation-based analysis of the algorithm performance is carried out, analysing influencing factors such as sensor misalignment, measurement errors and available surface information. In this course the results of the IFPR algorithm are compared to the FPR and the established ICP algorithm. For translational misalignment of up to 100 µm and rotational misalignment of up to 3 mrad, the proposed algorithm reconstructs a plane wavefront with registration errors a factor 4 and a divergent wavefront with registration errors a factor 12 smaller than the registration errors of the FPR and ICP algorithm. The considered wavefronts are reconstructed in a few iterations (3 to 4) by the proposed algorithm. The computation time of the algorithm on a personal computer is less than 500 ms making the algorithm suitable for time critical applications. For an overlap size between the measurements of 40 % of the sensor aperture area, the proposed algorithm achieves RMS registration errors smaller 10 nm enabling a qualitative assessment of high-end optical systems.
Future work concerns applications of the algorithm as well as further analysis of the algorithm with respect to the shape of the wavefront.