Computational Method for Wavefront Sensing Based on Transport-of-Intensity Equation

Recently the transport-of-intensity equation as a phase imaging method turned out as an effective microscopy method that does not require the use of high-resolution optical systems and a priori information about the object. In this paper we propose a mathematical model that adapts the transport-of-intensity equation for the purpose of wavefront sensing of the given light wave. The analysis of the influence of the longitudinal displacement z and the step between intensity distributions measurements on the error in determining the wavefront radius of curvature of a spherical wave is carried out. The proposed method is compared with the traditional Shack–Hartmann method and the method based on computer-generated Fourier holograms. Numerical simulation showed that the proposed method allows measurement of the wavefront radius of curvature with radius of 40 mm and with accuracy of ~200 μm.


Introduction
To date, the majority of devices and systems for determining the surface relief, threedimensional imaging, measuring and correcting aberrations, medical imaging and digital microscopy are based on methods for analyzing light fields by measuring intensity [1][2][3][4][5]. It is necessary to obtain all information about the light complex amplitude to characterize the light field, which is usually limited by the possibility of detecting only the intensity (amplitude) of the light beam. One of the most common ways to detect a phase argument of the complex amplitude is to convert the phase distribution to an amplitude distribution. The first to demonstrate such a possibility was L. Foucault, back in 1857, in his method known as the "Foucault knife-edge test" [6].
Wavefront sensors use the principles of geometric optics to transform the phase distribution into the intensity distribution [7][8][9]. Today the most commonly used wavefront sensor is a Shack-Hartmann sensor, which combines microlens array with a pixel detector. When the detector is placed in the focal plane of the microlens array, the position of the point focused by each microlens can be matched to a certain region of the wavefront. This method allows one to create a compact and simple scheme; therefore, Shack-Hartmann sensors are used in many fields of science [10][11][12][13][14]. However, processing the information recorded by the camera requires mathematical calculations to reconstruct the phase, which severely limits the performance and makes them applicable for certain tasks like atmospheric turbulence measurements only in pair with the state-of-art computers. It is also important that the parameters of the microlens array limit the spatial resolution, dynamic range and sensitivity of the sensor. Despite the fact that modern science has proposed several solutions to the abovementioned problems [15], technological processes still impose significant restrictions

Reconstruction of the Phase of a Coherent Field via Transport-of-Intensity Equation
Modern methods of phase reconstruction, which are the development of ideas formulated in works 40 years ago, allow finding of a complex field on the surface of an object from the measured intensity distributions on the detector. In this case, the diffraction integral used in practice describes the field at the detector through the field distribution on the object surface. As a rule, the calculation of the diffraction integral is reduced to the calculation of direct and inverse Fourier transforms. Let us show the role of the Fourier transform using the example of the transport-ofintensity equation, which describes the propagation of coherent radiation from an object to a detector, and is an exact solution of the wave equation.
It is known that a monochromatic coherent electromagnetic wave propagating in free space, in the absence of other sources, satisfies the Helmholtz equation where ∇ = {∂/∂x; ∂/∂y; ∂/∂z} is the Nabla operator for the three-dimensional space, r ⊥ = (x, y) vector radius in plane (x, y), orthogonal to the direction z, k is the wavenumber, U(r ⊥ , z) is the complex amplitude of the field of the light wave, which is expressed as where I(r ⊥ , z) is the intensity of the wave, and φ(r ⊥ , z) is its phase. Using the complex field amplitude (2), we can obtain Equation (1) in the paraxial approximation where ∇ ⊥ = {∂/∂x; ∂/∂y} is the Nabla operator for the two-dimensional space.
Separating the real part of Equation (3), we obtain Equation (4) is called the transport-of-intensity equation (TIE). It connects the intensity I and its longitudinal derivative ∂I/∂z with the phase φ of the light wave. Note that the TIE was originally derived from the wave equation [27], and later derived from Poynting's theorem [28] and other fundamental relations [29].
Teague [27] proposed one of the approaches to solving the resulting equation. Its essence lies in the application of the Helmholtz expansion theorem [30], which allows representing the vector field I(r ⊥ , z)∇ ⊥ φ(r ⊥ , z) in the form where ψ(r ⊥ , z) is the scalar potential, A(r ⊥ , z) is the vector potential. Neglecting the second term [∇ × A(r ⊥ , z)], expression (4) can be simplified into two standard Poisson equations.
We can present Equation (6) with respect to the phase φ(r ⊥ , z) [28,31] as One of the most popular methods for solving Equation (7) is the method based on the expression of differential operators in terms of the Fourier transform [31]. This solution is simple and effective. Let us consider the main idea of the method in more detail. Using the properties of the Fourier transform, one can express the sum of the partial derivatives of the m-th and n-th powers [32]. where F иF −1 are operators of direct and inverse Fourier transforms, respectively; k x,y = 2πν x,y are the frequency coefficients, ν x,y are the frequency grids. Then the differential operators in Equation (8) can be written as It should be noted that, for the application of relations (10) and (11), in view of the periodic nature of the Fourier transform, it is necessary to fulfill the periodic boundary conditions, under which the value of the function at the boundaries must be cyclically repeated [33].

Proposed Method
This section describes the mathematical model for determining the RoC of a spherical wave with a Gaussian intensity ( Figure 1). We used the spherical wave as an object, which can be characterized by the RoC, while the gaussian intensity distribution is the most common distribution when it comes to coherent beams, such as laser beams. The novelty of the proposed approach is in the combination of two methods: TIE to reconstruct the phase of the complex field amplitude from a set of intensities and a geometric method (to be explained below) to calculate the RoC from a known complex field amplitude.
Then the differential operators in Equation (8) can be written as , , , , It should be noted that, for the application of relations (10) and (11), periodic nature of the Fourier transform, it is necessary to fulfill the perio conditions, under which the value of the function at the boundaries mus repeated [33].

Proposed Method
This section describes the mathematical model for determining the RoC wave with a Gaussian intensity ( Figure 1). We used the spherical wave as an can be characterized by the RoC, while the gaussian intensity distribution is mon distribution when it comes to coherent beams, such as laser beams. T the proposed approach is in the combination of two methods: TIE to reconst of the complex field amplitude from a set of intensities and a geometric explained below) to calculate the RoC from a known complex field amplitu    Obtaining a set of complex amplitudes U(r ⊥ , z i ), . . . , U(r ⊥ , z i+1 ) by propagating the original field through space using the angular spectrum method [34], or acquiring a set of intensities I(r ⊥ , z i ) from measurements in a physical experiment (dashed block); 2.
Calculation of the theoretical RoC of the wavefront R t (z) using the geometric method; 3.
Calculation of RoC of the wavefront R TIE (z) using the initial intensities (step 1) and phases obtained in step 3 by the geometric method; 5.
Comparative analysis of R t (z) and R TIE (z). Next, the algorithm is explained more thoroughly. Note that to solve TIE (Equation (8)), it is necessary to know the longitudinal derivative with respect to the intensity ∂I/∂z, which can be approximated by the finite difference method (Equation (12)), if at least two intensity distributions and the distance ∂z between them are known ( Figure 1).
Therefore, the first step of the proposed algorithm ( Figure 2) is to obtain a set of complex amplitudes ( , ) by propagating the initial field ( , 0) U r z ⊥ = through space using angular spectrum method [34]. The propagation of the initial field is done only for simulation purposes: To attain intensities used later for TIE and for comparison in the end. For a physical experiment, acquiring a set of intensities in different planes will suffice.
At the second step, the theoretical radius of curvature Rt is calculated using the geometric method. Let's consider it in more detail. It is known [35] that the radius of the circle R is related to the sag s and the chord l by the following relation The chord l in our case is the width of the Gaussian function in level Imax(x)/e 2 ( Figure  3a). The search for the sag s assumes that the complex field amplitude is zero where I<Imax/e 2 . Therefore, it is necessary to form an aperture with a diameter equal to the chord length l and superimpose it on the wrapped phase ϕw (Figure 3b). The sag s is found as the phase maximum after its unwrapping operation [36]. Thus, using Equation (13), one can determine the RoC of the wavefront. Next, the algorithm is explained more thoroughly. Note that to solve TIE (Equation (8)), it is necessary to know the longitudinal derivative with respect to the intensity ∂I/∂z, which can be approximated by the finite difference method (Equation (12)), if at least two intensity distributions and the distance ∂z between them are known ( Figure 1).

∂I(z) ∂z
Therefore, the first step of the proposed algorithm ( Figure 2) is to obtain a set of complex amplitudes U(r ⊥ , z i ), . . . , U(r ⊥ , z i+1 ) by propagating the initial field U(r ⊥ , z = 0) through space using angular spectrum method [34]. The propagation of the initial field is done only for simulation purposes: To attain intensities used later for TIE and for comparison in the end. For a physical experiment, acquiring a set of intensities in different planes will suffice.
At the second step, the theoretical radius of curvature R t is calculated using the geometric method. Let's consider it in more detail. It is known [35] that the radius of the circle R is related to the sag s and the chord l by the following relation The chord l in our case is the width of the Gaussian function in level I max (x)/e 2 (Figure 3a). The search for the sag s assumes that the complex field amplitude is zero where I < I max /e 2 . Therefore, it is necessary to form an aperture with a diameter equal to the chord length l and superimpose it on the wrapped phase φ w (Figure 3b). The sag s is found as the phase maximum after its unwrapping operation [36]. Thus, using Equation (13), one can determine the RoC of the wavefront. Note that the calculation of the theoretical RoC Rt was done only for simulation purposes: To compare with the proposed method and to figure out its limitations and applicability. In a physical experiment there will be no calculation of the theoretical RoC. In this case other methods will be used to evaluate the performance of the proposed method (e.g., interferometry, Shack-Hartmann sensors and others).
Finally, in the third step of the algorithm, we find the derivative using Equation (12), solve TIE (Equation (8)) and get the already unwrapped phase ϕuw,TIE (one of the advantages of TIE). Next, we calculate the chord l and the value of the sag sTIE and we calculate RTIE using Equation (13).
At the last step, a comparison of Rt and RTIE is carried out, the results of which will be discussed in the next section.

The Limits of Applicability of Methods for Measuring the Curvature of the Wavefront
The proposed method is compared with the traditional method based on the Shack-Hartmann sensor and the holographic method using numerical simulation. The simulation was carried out in the Python programming language (http://www.python.org, accessed on 05 May 2021), which in recent years has become a powerful tool for fundamental and applied research in astronomy [37] and other related fields [38]. However, such comparison could be implemented in any other programming language (C++, MATLAB, etc.). For comparison, we took parameters of commercially available devices, such as Shack-Hartmann wavefront sensor and spatial light modulator with Full HD resolution.

Shack-Hartmann Sensor
Wavefront sensors are devices that measure the deviation of an optical wavefront from a plane or sphere. Let us calculate the minimum and maximum RoC that such a sensor can detect using the scheme shown in Figure 4.
The minimum radius Rmin can be calculated from the condition of the maximum angle θmax, by which one of the elements of the microlens array can deflect the incident beam, provided that the spot focused by the microlens does not fall on the area which is matched with another microlens. This angle is determined by the Equation (14): Note that the calculation of the theoretical RoC R t was done only for simulation purposes: To compare with the proposed method and to figure out its limitations and applicability. In a physical experiment there will be no calculation of the theoretical RoC. In this case other methods will be used to evaluate the performance of the proposed method (e.g., interferometry, Shack-Hartmann sensors and others).
Finally, in the third step of the algorithm, we find the derivative using Equation (12), solve TIE (Equation (8)) and get the already unwrapped phase φ uw,TIE (one of the advantages of TIE). Next, we calculate the chord l and the value of the sag s TIE and we calculate R TIE using Equation (13).
At the last step, a comparison of R t and R TIE is carried out, the results of which will be discussed in the next section.

The Limits of Applicability of Methods for Measuring the Curvature of the Wavefront
The proposed method is compared with the traditional method based on the Shack-Hartmann sensor and the holographic method using numerical simulation. The simulation was carried out in the Python programming language (http://www.python.org, accessed on 5 May 2021), which in recent years has become a powerful tool for fundamental and applied research in astronomy [37] and other related fields [38]. However, such comparison could be implemented in any other programming language (C++, MATLAB, etc.). For comparison, we took parameters of commercially available devices, such as Shack-Hartmann wavefront sensor and spatial light modulator with Full HD resolution.

Shack-Hartmann Sensor
Wavefront sensors are devices that measure the deviation of an optical wavefront from a plane or sphere. Let us calculate the minimum and maximum RoC that such a sensor can detect using the scheme shown in Figure 4.
The minimum radius R min can be calculated from the condition of the maximum angle θ max , by which one of the elements of the microlens array can deflect the incident beam, provided that the spot focused by the microlens does not fall on the area which is matched with another microlens. This angle is determined by the Equation (14): where d is the microlens diameter, f is its equivalent focal length. For example, for a Thorlabs WFS 300-14AR sensor with d = 300 µm and f = 14.2 mm, the θ max value is defined as 10.56 mrad.
where d is the microlens diameter, f is its equivalent focal length. For example, for a Thorlabs WFS 300-14AR sensor with d = 300 μm and f = 14.2 mm, the θmax value is defined as 10.56 mrad. The maximum radius Rmax can be calculated from the condition of the minimum detectable angle θmin, it is determined by the Equation (15) min min , y f Δ θ = (15) where ∆ymin is the minimum displacement of the light spot in the plane of the CCD that can be detected [39]. For CCD, it can be calculated as ∆ymin = a/20, where a = 4.65 μm is the pixel size of the WFS 300-14AR sensor. Then we get that θmin = 16.37 μrad. Based on the maximum and minimum angles, it is easy to calculate the minimum and maximum RoC that the sensor can register. The maximum radius Rmax can be calculated as follows where rdet is the radius of the aperture in which the wavefront curvature is measured. The maximum radius of the aperture that can be inscribed in the WFS 300-14AR sensor aperture is 4.76 mm, then Rmax = 360 m. The minimum radius Rmin is calculated from the condition θ ≤ θmax, where θ is the angle by which the wavefront is deflected by the microlens located at the edge of the array. For the considered WFS 300-14AR sensor Rmin = 340.6 mm. At the same time, proceeding from the accuracy of wavefront measurements (λ/50 for the selected WFS), we can determine the accuracy of the RoC measurement. The change of the sag s (Equation (13)) of the local wavefront caused by increasing/decreasing of defocus value by λ/50 will result in changing of the RoC by 517 μm, while working around Rmin.
Shack-Hartmann wavefront sensors output characteristic is linear on its entire dynamic range and its dynamic range is determined mostly by the microlens array used in the development. However, it was shown that the dynamic range can be expanded using different techniques. Namely the defocus aberration--which can be directly translated to RoC of wavefront--was measured with a maximum of 11% error over a range eight times larger than the microlens-bound definition [40]. However, such an error in estimation of the defocus value will translate to even bigger error in determining of the wavefront RoC (∆R = 5mm for the RoC R = 50 mm for the parameters of the wavefront sensor given above). The maximum radius R max can be calculated from the condition of the minimum detectable angle θ min , it is determined by the Equation (15) where ∆y min is the minimum displacement of the light spot in the plane of the CCD that can be detected [39]. For CCD, it can be calculated as ∆y min = a/20, where a = 4.65 µm is the pixel size of the WFS 300-14AR sensor. Then we get that θ min = 16.37 µrad.
Based on the maximum and minimum angles, it is easy to calculate the minimum and maximum RoC that the sensor can register. The maximum radius R max can be calculated as follows where r det is the radius of the aperture in which the wavefront curvature is measured. The maximum radius of the aperture that can be inscribed in the WFS 300-14AR sensor aperture is 4.76 mm, then R max = 360 m. The minimum radius R min is calculated from the condition θ ≤ θ max , where θ is the angle by which the wavefront is deflected by the microlens located at the edge of the array. For the considered WFS 300-14AR sensor R min = 340.6 mm. At the same time, proceeding from the accuracy of wavefront measurements (λ/50 for the selected WFS), we can determine the accuracy of the RoC measurement. The change of the sag s (Equation (13)) of the local wavefront caused by increasing/decreasing of defocus value by λ/50 will result in changing of the RoC by 517 µm, while working around R min .
Shack-Hartmann wavefront sensors output characteristic is linear on its entire dynamic range and its dynamic range is determined mostly by the microlens array used in the development. However, it was shown that the dynamic range can be expanded using different techniques. Namely the defocus aberration--which can be directly translated to RoC of wavefront--was measured with a maximum of 11% error over a range eight times larger than the microlens-bound definition [40]. However, such an error in estimation of the defocus value will translate to even bigger error in determining of the wavefront RoC (∆R = 5 mm for the RoC R = 50 mm for the parameters of the wavefront sensor given above).

Holographic Method Based on a Spatial Light Modulator
There is a compact device for measuring the wavefront based on highly efficient computer-generated holograms (CGH) in addition to traditional wavefront sensors [41,42]. The dimensions and resolution of the modulator used to output the CGH in such devices would determine the maximum spatial frequency of the CGH. Based on this, it is possible to determine the dynamic range of the measured values of the RoC for the method based on the spatial light modulator. According to [43], the resolution of the object recorded on the hologram should be four times less than the resolution of the final hologram. However, the expression is valid for the hologram recording in a physical medium. In the case when the object is a CGH, it is enough to display it without distortion on the modulator. This can be achieved if at least 4 pixels of the hologram fall on one period of the interference pattern (structure) [44].
The CGH itself consists of two elements: the blazed diffraction grating and the phase function imposed on it. When measuring the RoC, the Zernike polynomial corresponding for defocus (Z 2 2 ), acts as a phase structure. In the limiting case, when work is carried out in zero order, the blazed diffraction grating can be neglected. The minimum RoC that can be registered without distortion can be calculated using the phase structure. With a modulator resolution of 1920 × 1080 pixels R min = 170 mm. With a further decrease in the radius on the final CGH, there are areas with a period of T ≤ 3 pixels. The use of such CGHs will lead to an increase in the measurement error of RoC.
In general, the holographic method output characteristic is linear for a certain interval, which is determined by the values of the aberration modes, used during hologram synthesis [45,46]. Increasing the dynamic range by using the higher value of the aberration modes will decrease the accuracy of the method. The accuracy for wavefront aberration modes measurement vary from λ/10 [47] to λ/50 [18,21], which for the given parameters of the modulator (aperture 6912 mm for the pixel size of 6.4 µm) translates to RoC measurement accuracy ∆R = 201 µm and ∆R = 41 µm, respectively.

Method Based on the Transport-of-Intensity Equation
This section presents the results of three numerical simulations carried out using the mathematical model described in Section 3. In all numerical simulations, the following parameters were used: radiation wavelength λ = 650 nm; field size w × h = 1500 × 1500 pixels; sampling step 5.04 µm. The intensity and phase in the z = 0 plane were specified by the following relations φ(x, y, 0) = x 2 + y 2 + R 2 0 , where l 0 = 750 pixels, R 0 value depends on the simulation. In the first numerical simulation, the linearity of the proposed method was investigated for RoC R 0 = 100 mm, the range of distances z = 0, ∂z, 2∂z, . . . , R 0 , where ∂z = 1 mm. The obtained dependence ( Figure 5) shows that the nonlinearity increases in direct proportion to the distance from the point z = R 0 , that is, the error in determining RoC is proportional to the value of the chord l (the width of the Gaussian function). During the second numerical simulation, the influence of the step value ∂z was investigated for the following parameters: R0 = 100 mm; z = 80, 80 + ∂z, 80 + 2∂z, …, R0; ∂z = 0.1, 1, 5, 10 mm. Earlier [48] it was shown that the value of the step ∂z should be large enough to cover the influence of noise and, at the same time, small enough so as not to violate the linearity of the derivative approximation (Equation (12)). In our simulation, During the second numerical simulation, the influence of the step value ∂z was investigated for the following parameters: R 0 = 100 mm; z = 80, 80 + ∂z, 80 + 2∂z, . . . , R 0 ; ∂z = 0.1, 1, 5, 10 mm. Earlier [48] it was shown that the value of the step ∂z should be large enough to cover the influence of noise and, at the same time, small enough so as not to violate the linearity of the derivative approximation (Equation (12)). In our simulation, there is no noise, so the step value is limited only from above, which is observed in the obtained dependence ( Figure 6). More thorough study on the accuracy in the displacement of the between planes where intensity distributions were acquired was carried out in [49]. During the second numerical simulation, the influence of the step value ∂z was investigated for the following parameters: R0 = 100 mm; z = 80, 80 + ∂z, 80 + 2∂z, …, R0; ∂z = 0.1, 1, 5, 10 mm. Earlier [48] it was shown that the value of the step ∂z should be large enough to cover the influence of noise and, at the same time, small enough so as not to violate the linearity of the derivative approximation (Equation (12)). In our simulation, there is no noise, so the step value is limited only from above, which is observed in the obtained dependence ( Figure 6). More thorough study on the accuracy in the displacement of the between planes where intensity distributions were acquired was carried out in [49]. Note that it is difficult to find the exact value of the chord l near R0 due to the effect of sampling. So, at z = 97 mm the chord l = 5 px, and the values in these 5 pixels differ from each other by 4•10 2 (for Imax(z = 0) = 1.0). That is, to improve the accuracy in the region z ≅ R0, interpolation is necessary, which will allow obtaining a more accurate value of the chord l. Note that it is difficult to find the exact value of the chord l near R 0 due to the effect of sampling. So, at z = 97 mm the chord l = 5 px, and the values in these 5 pixels differ from each other by 4·10 2 (for I max (z = 0) = 1.0). That is, to improve the accuracy in the region z ∼ = R 0 , interpolation is necessary, which will allow obtaining a more accurate value of the chord l.
In the last numerical simulation, R min was determined for the comparison with the methods presented above. The distance (R 0 -z) was fixed: (R 0 -z) = 15.5 mm-that is, 15.5 mm in front of R 0 ; step ∂z = 1 mm; R 0 = 20 . . . 100 mm. The resulting dependence (Figure 7) demonstrates the possibility of using this mathematical model at R min ∼ = 40 mm. In the last numerical simulation, Rmin was determined for the comparison with the methods presented above. The distance (R0-z) was fixed: (R0-z) = 15.5 mm-that is, 15.5 mm in front of R0; step ∂z = 1 mm; R0 = 20 … 100 mm. The resulting dependence ( Figure  7) demonstrates the possibility of using this mathematical model at Rmin ≅ 40 mm.

Results and Discussion
The main data from Section 4 was assembled in the Table 1. All of the errors were calculated based on the same radiation wavelength λ = 650 nm. It can be seen, that while the proposed method does not achieve the highest accuracy of wavefront aberrations measurement, it still allows one to measure wavefront RoC with curvature of less than 60mm with adequate accuracy in comparison to other methods. So, we can say that the

Results and Discussion
The main data from Section 4 was assembled in the Table 1. All of the errors were calculated based on the same radiation wavelength λ = 650 nm. It can be seen, that while the proposed method does not achieve the highest accuracy of wavefront aberrations measurement, it still allows one to measure wavefront RoC with curvature of less than 60mm with adequate accuracy in comparison to other methods. So, we can say that the proposed method can expand the lowest range of wavefront RoC measurement, where other methods start showing high errors. In addition, all of this can be achieved with just a camera and software without the need of using focusing optics, which itself can introduce aberrations into the beam. Although it is worth noting that, for all the methods, the accuracy increases with an increase of RoC. That happens because the same change of the sag s translates to smaller relative change in wavefront RoC. Figure 8 shows the graphical comparison of the dynamic range limits of the three methods.  Certainly, there are devices with higher resolution (up to 4 and 8 K), that will allow to decrease the Rmin, while increasing the accuracy. Likewise, the same changes could be done to the proposed method. Thus, even using more common commercially available devices one can achieve better results while implementing the proposed method. In this work we showed the principle of the proposed method for the case of the wavefront RoC measurement with symmetric wavefront profile (defocus aberration mode). However, for the TIE it does not matter if the wavefront profile is symmetric or asymmetric as long as the wavefront function satisfies the boundary conditions. Naturally the asymmetry will impact the RoC calculation process, but the algorithm can be adapted by introducing the polynomials into the calculation process [50].

Conclusions
In this study we proposed a method for wavefront sensing of the light radiation based on the transport-of-intensity equation. The proposed method was studied using numerical simulation and compared to the common devices, such as Shack-Hartmann wavefront sensors and Holographic wavefront sensors. The comparison was done with the radius of curvature of the wavefront of spherical wave. The resulting error in measuring the wavefront radius of curvature R0 = 40 mm at λ = 650 nm was 0.20 mm. The pro- Certainly, there are devices with higher resolution (up to 4 and 8 K), that will allow to decrease the R min , while increasing the accuracy. Likewise, the same changes could be done to the proposed method. Thus, even using more common commercially available devices one can achieve better results while implementing the proposed method. In this work we showed the principle of the proposed method for the case of the wavefront RoC measurement with symmetric wavefront profile (defocus aberration mode). However, for the TIE it does not matter if the wavefront profile is symmetric or asymmetric as long as the wavefront function satisfies the boundary conditions. Naturally the asymmetry will impact the RoC calculation process, but the algorithm can be adapted by introducing the polynomials into the calculation process [50].

Conclusions
In this study we proposed a method for wavefront sensing of the light radiation based on the transport-of-intensity equation. The proposed method was studied using numerical simulation and compared to the common devices, such as Shack-Hartmann wavefront sensors and Holographic wavefront sensors. The comparison was done with the radius of curvature of the wavefront of spherical wave. The resulting error in measuring the wavefront radius of curvature R 0 = 40 mm at λ = 650 nm was 0.20 mm. The proposed method enables to expand the dynamic range of the radius of curvature of the wavefront measurement.
The novelty of the proposed method lies in the fact that the transport-of-intensity equation was for the first time used to directly measure the wavefront of the light wave (i.e., wavefront sensing), unlike its previous uses in microscopy imaging. The estimates obtained in the article could be used in the development of optoelectronic systems designed to register, process and display information about optical radiation.

Data Availability Statement:
The data is available from the corresponding author upon request.