- Email: [email protected]
Wavefront reconstruction using iterative discrete Fourier transforms with Fried geometry
ARTICLE IN PRESS
Optik
Optics
Optik 117 (2006) 77–81 www.elsevier.de/ijleo
Wavefront reconstruction using iterative discrete Fourier transforms with Fried geometry Huanqing Guo, Zhaoqi Wang Institute of Modern Optics, Nankai University, Tianjin 300071, PR China Received 29 March 2005; accepted 15 June 2005
Abstract A new kind of wavefront reconstruction algorithm is proposed in this paper, which can estimate the wavefront from its orthogonal gradients at different sample points. We model the gradients measured by wavefront sensor with the Fried geometry. Iterative discrete Fourier transforms and the correspondent Fried inverse filter in spatial frequency domain are introduced. No special boundary control is required to process the boundary condition in this algorithm. The stimulation and the experimental results show that the new algorithm has strong capability to fit the wavefront and it is more available if the sample points are few. r 2005 Elsevier GmbH. All rights reserved. Keywords: Wavefront reconstruction; Fourier transform; Fried geometry
1. Introduction The use of the discrete Fourier transform (DFT) for wavefront reconstruction from slope map was firstly proposed by Freischlad and Koliopoulos [1]. They applied DFT method to square or rectangular domain, and the method to solve the inverse filters in spatial frequency domain for different geometries were derived in their paper. Their work had been improved and used in interferogram analysis and in adaptive optics systems of astronomy [2–5]. Poyneer et al. suggested the boundary and extension method to fill in missing slope information for large circular apertures with use of DFT method [2]. As for the domain with any arbitrary shape, Roddier did good work in wavefront reconstruction using iterative DFT. Roddier’s algorithm came from Gerchberg’s extrapolation algorithm [6] beyond the Corresponding author. Fax: +8622 23508332.
E-mail address: [email protected] (Z. Wang). 0030-4026/$ - see front matter r 2005 Elsevier GmbH. All rights reserved. doi:10.1016/j.ijleo.2005.06.003
boundaries [3,4]. We name Roddier and Gerchberg’s work with R–G algorithm. In this paper, we introduced both Fried geometry [7] and DFT method to reconstruct the wavefront. And we combined inverse filter and iterative Fourier transforms to give more exact and easy estimate of the wavefront. Our work can be looked as the modification and the improvement of R–G algorithm.
2. Fourier iterative algorithm with Fried geometry A flow chart of our modified new iterative algorithm is shown in Fig. 1. In this chart, the input data is the gradient matrices which can be measured by wavefront sensor, such as the Hartmann–Shack (H–S) sensor. The output of the procedure shown in Fig. 1 is the estimated wavefront with acceptable error. In order to explain the work flow chart, we take the H–S wavefront sensor as
ARTICLE IN PRESS 78
H. Guo, Z. Wang / Optik 117 (2006) 77–81 W(m-1,n-1)
Gradient matrixes with zero value out of the pupil
DFT
Put original gradient matrixes back inside the pupil
Filter for Fried geometry in spatial frequency domain
Recalculated gradient matrixes with Fried geometry
W(m-1,n)
W(m,n-1)
W(m+1,n-1)
W(m,n)
W(m+1,n)
Sx[m,n]
Inverse DFT Error acceptable W(m-1,n+1)
Estimated wavefront
W(m,n+1)
Sy[m,n]
W(m+1,n+1)
Fig. 3. The Fried sensor geometry.
Fig. 1. Flow chart of the new iterative discrete Fourier transform algorithm used to reconstruct a wavefront from the measured wavefront gradients.
have the difference in x- and y-directions S x ¼ rx W ðm; nÞ ¼ 1=2½W ðm þ 1; n þ 1Þ þ W ðm þ 1; nÞ 1=2½W ðm; n þ 1Þ þ W ðm; nÞ,
ð2Þ
S y ¼ ry W ðm; nÞ ¼ 1=2½W ðm þ 1; n þ 1Þ þ W ðm; n þ 1Þ 1=2½W ðm þ 1; nÞ þ W ðm; nÞ. ð3Þ
Fig. 2. Hartmann–Shack grids.
an example. Figs. 2 and 3 show the H–S grids and estimated points of wavefront W(m,n) and the gradients Sx, Sy measured by H–S sensor in Fried sensor geometry, respectively. We want to determine the wavefront W(m,n) over the measurement domain within the pupil from its noisy discrete gradient S, given by Sðm; nÞ ¼ rW ðm; nÞ þ nðm; nÞ with S¼
Sx Sy
! ;
n¼
nx ny
(1)
! ,
where n represents the additive noise to the measured gradient and it will not be discussed in this paper. The first coordinate is defined as the y-coordinate which is the columns in a numeric matrix. And the rows of the numeric matrix are along the x-coordinate. For the lowright sensor locations in Fig. 3 of Fried geometry, we
In the work flow chart that Fig. 1 shows, the gradient matrices are recalculated from the estimated wavefront points with Eqs. (2) and (3). The DFT is applied to both the equations with use of the shift property of the Fourier transform, Eqs. (2) and (3) become 1 i2ðk þ lÞp i2kp ~ lÞ exp DFTfS x g ¼ wðk; þ exp 2 N N i2lp exp 1 , ð4Þ N 1 i2ðk þ lÞp i2lp ~ lÞ exp þ exp DFTfS y g ¼ wðk; 2 N N i2kp exp 1 , ð5Þ N where the number of sample grids is N N discrete signal, k, l are the coordinate in spatial frequency domain, and i2 ¼ 1. We can derive the inverse filter by combine Eqs. (4) and (5). The following formula is the ~ lÞ for Fried geometry: solved wðk; 8 0; k ¼ l ¼ 0; k ¼ l ¼ N=2 > >
i2kp > > < exp N 1 exp i2lp N þ 1 DFTfS x g
~ lÞ ¼ wðk; 1 exp i2kp þ 1 DFTfS y g þ exp i2lp > N N > > >
: 2 kp 2 lp 2 kp 1 else: 8 sin N cos2 lp N þ sin N cos N (6) ~ lÞ will Taking the inverse Fourier transform of wðk; ~ produce the estimated wavefront wðm; nÞ. There are a key assumption when we calculate DFT{Sx} and DFT{Sy} in Eq. (6). It is necessary to
ARTICLE IN PRESS H. Guo, Z. Wang / Optik 117 (2006) 77–81
have difference data consistent with a periodic reconstructed wavefront of period N is x- and y-directions. From Fig. 3, we can conclude that only ðN 1Þ ðN 1Þ array of x and y differences are practically measured by the sensor. In order to calculate DFT during the iterative procedure, we must process reasonably the missing information to obtain the N N array of differences. Furthermore, when the aperture is circular or even arbitrary, the estimated wavefront cross the boundary of the aperture will be biased if no boundary control algorithm is used [2]. Eq. (6) gives the formula of the filter for Fried geometry in spatial frequency domain. This formula will enter the loop shown in Fig. 1. For certain size of matrices, the exponential and the trigonometric operation in Eq. (8) can be calculated before starting the DFT iterative operations. In our new iterative DFT algorithm, boundary control can be done during the course of iteration and
79
has the fewest iterative times to arrive at the minimum reconstruction error. After this, in order to continue the iteration, we put back the measurement data inside the aperture and do the iteration sequentially, until the wavefront estimate is error acceptable. The measurement of current error in the iterative reconstruction process is to calculate the root mean square (RMS) difference. If the simulation is considered, the RMS difference is between the real wavefront and the estimated wavefront 2RR
3 ~ ðx; yÞ2 dx dy 1=2 ½W real ðx; yÞ W 6 7 RR RMS ¼ 4 A 5 , dx dy
(8)
A
where A represents the area of the circular aperture. If real data is considered, the RMS difference is between the measured gradients and the estimated gradients
2RR n 31=2
2
2 o S xreal ðx; yÞ S~ x ðx; yÞ þ Syreal ðx; yÞ S~ y ðx; yÞ dx dy 6A 7 7 . RR RMS ¼ 6 4 5 2 dx dy
(9)
A
it is not necessary to derive additional boundary method. The iterative procedure produces exact approximation cross the boundary of the estimated wavefront. In Fig. 2, the black bold square is ðN 1Þ ðN 1Þ H–S sensor array. The circular boundary tangential with the square is used to indicate those lenslets (also called ðN 1Þ ðN 1Þ) whose outputs are used in the estimation of the wavefront. The biggest square outside the black bold one adds some complementary points and its size is M M ðM4NÞ. In order to accelerate the speed of the computation, M is often set to be the integral power of 2. Before the first iteration, all the measurement gradients are within the circular aperture and the outside domain of it is padded with zero. After the first iteration, we recalculated two ðM 1Þ ðM 1Þ gradient matrices S½ðM 1Þ ðM 1Þ from the first ~ ðm; nÞ by use of Eqs. (2) and (3). estimated wavefront W And we set ~ ðM; nÞ, S x ðM; nÞ ¼ S y ðM; nÞ ¼ W ~ ðn; MÞ S x ðn; MÞ ¼ S y ðn; MÞ ¼ W
ð7Þ
for 1pnpM. Now the new gradient matrices have the same size M M with the original two. In fact, other values can be used to complement the missing row and column, such as zero or constant. At present, we are not so sure about which is the best value for the missing gradients data. The complement method with Eq. (7)
The units of both Eqs. (8) and (9) are the unit of length, i.e. micron or wavelength.
3. Simulation results We develop a simulation that the wavefront is given by the 12th fringe Zernike polynomial (ZP), say, spherical aberration with its coefficient C varied from 0.1 to 6 mm. Its expression is Z ¼ Cð6x4 þ 6y4 þ 12x2 y2 6x2 6y2 þ 1Þ. The circular pupil inside a square sample area is considered and there is no noise added in our work. Fig. 4(a) and (b) show the simulation results. In order to compare with these results, the reconstruction results of R–G algorithm are shown in Fig. 4 too. According to Marechal limit, the RMS fitting error will be acceptable if it is lower than l=14. The wavelength l ¼ 633 nm in our simulation. This limit is also drawn with blue dashed line in Fig. 4. In the case that Fig. 4(a) shows, the RMS errors of both algorithms decrease while the sample points increase. And the new algorithm has lower RMS error than R–G algorithm on every sample level. Even the fitting error of the new algorithm at 30 30 is less than that of R–G algorithm at 100 100. If the sample grids
ARTICLE IN PRESS 80
H. Guo, Z. Wang / Optik 117 (2006) 77–81
Fig. 5. The convergence of the wavefront’s RMS fitting error comes up with the iterative times of Fourier transforms.
to its minimal value within 10 iterative times. The new algorithm always has greater RMS fitting error after the first iteration. If a few more iterations are carried out, new algorithm will give better results than R–G algorithm. For example, it requires about 15 iterations to get to the equal fitting error of R–G algorithm if the sample grids is 100 100. In fact, we also found the similar iterative cases in other spherical aberration coefficients as shown in Fig. 5.
Fig. 4. Simulation results with the new iterative algorithm (blue color symbols) and R–G algorithm (red color symbols) for spherical aberration 0.1–6 mm: (a) 200 iterative times and (b) 10 iterative times.
is 10 10 or lower, R–G algorithm produces obvious fitting error if the coefficient is more than 0.2 mm. And the new algorithm produces the same error when the coefficient gets to 1.4 mm. In addition, the new algorithm gives results that can be totally acceptable when the sample grids is 30 30. At the same time, the results of R–G algorithm can be acceptable only if the coefficient is less than 2.5 mm. The iterative DFT times of Fig. 4(a) is 200, which perhaps will be not suitable for real-time system. Then we just restrict the iterative times to 10 in Fig. 4(b). In this figure, the new algorithm is no better than R–G algorithm if the sample grids is 100 100. Even the RMS error under this sample is bigger than that of 60 60 and almost equal to 30 30 using the new algorithm. It is clear that iterative times will influence the results of the algorithm. The iterative times of different sample density for 1 mm spherical aberration coefficient is shown in Fig. 5. In this figure, we can see that the R–G algorithm converges more quickly than the new algorithm and gets
4. Results of real data tests In order to observe the behavior of the new algorithm, we use the Hartmann–Shack (H–S) wavefront sensor to measure the wavefront aberration of the eye. With use of the experimental setup that paper [8] showed, one human eye was measured and 75 facula array pictures were gotten within 3 s. We analyzed 25 pictures which were obtained in the 2nd second of the measurement. The sample points were 18 18 and the pupil size was 5.6 mm.15, 35, and 65 ZPs were also used, respectively, as the basis functions to fit the measured wavefront. The iterative times of both the R–G algorithm and the new algorithm were confined to 30 times. The mean RMS errors of five processing methods are shown in Fig. 6. It is clear that the RMS difference between the gradients measured and gradients estimated has minimum value using the new iterative DFT algorithm. It decreases about 33% compare with 65 ZPs fitting and about 25% compare with R–G algorithm.
5. Conclusions The new iterative DFT algorithm developed here has the following features. The first feature inherits from the Roddier–Gerchberg algorithm, which can used to
ARTICLE IN PRESS H. Guo, Z. Wang / Optik 117 (2006) 77–81
81
No. 60438030), the Chinese Ministry of Education’s Nankai University and Tianjin University Cooperation Foundation and the Key Research Foundation of Scientific and Technical Committee of Tianjin City of China (No. 033183711).
References Fig. 6. RMS error for measuring the wavefront of the eye with H–S sensor.
analyze the wavefront in circular or the arbitrary pupil. Secondly, no special boundary control algorithm is needed in this algorithm for circular or arbitrary pupil. The production of periodic signal required by Fourier transform is solved with Gershberg’s extrapolation method, which will give accurate approximation of boundary differences during the iterative procedure. Thirdly, without the noise, the new algorithm can reconstruct the wavefront more accurately than its previous algorithm, especially in the case of low number of the sample points. If the sample density is high, more iterative times will be required. At last, the measured orthogonal slopes centered in-between the estimated points on the wavefront, which is determined by Fried geometry.
Acknowledgments This work was supported by the National Nature Science Foundation of China (Key Research Project
[1] K.R. Freischlad, C.L. Koliopoulos, Modal estimation of a wave front from difference measurements using the discrete Fourier transform, J. Opt. Soc. Am. A 3 (11) (1986) 1852–1861. [2] L.A. Poyneer, D.T. Gavel, J.M. Brase, Fast wave-front reconstruction in large adaptive optics systems with use of the Fourier transform, J. Opt. Soc. Am. A 19 (10) (2002) 2100–2111. [3] C. Roddier, F. Roddier, Interferogram analysis using Fourier transform techniques, Appl. Opt. 26 (9) (1987) 1668–1673. [4] F. Roddier, C. Roddier, Wavefront reconstruction using iterative Fourier transforms, Appl. Opt. 30 (11) (1991) 1325–1327. [5] A. Dubra, C. Paterson, C. Daninty, Wave-front reconstruction from shear phase maps by use of the discrete Fourier transform, Appl. Opt. 43 (5) (2004) 1108–1113. [6] R.W. Gerchberg, Super-resolution through error energy reduction, Opt. Acta 21 (1974) 709–720. [7] D.L. Fried, Least-Square fitting a wave-front distortion estimate to an array of phase-difference measurements, J. Opt. Soc. Am. 67 (1977) 370–375. [8] H. Guo, Z. Wang, Q. Zhao, W. Quan, Y. Wang, Individual eye model based on wavefront aberration, Optik 116 (2) (2005) 80–85.
Optik
Optics
Optik 117 (2006) 77–81 www.elsevier.de/ijleo
Wavefront reconstruction using iterative discrete Fourier transforms with Fried geometry Huanqing Guo, Zhaoqi Wang Institute of Modern Optics, Nankai University, Tianjin 300071, PR China Received 29 March 2005; accepted 15 June 2005
Abstract A new kind of wavefront reconstruction algorithm is proposed in this paper, which can estimate the wavefront from its orthogonal gradients at different sample points. We model the gradients measured by wavefront sensor with the Fried geometry. Iterative discrete Fourier transforms and the correspondent Fried inverse filter in spatial frequency domain are introduced. No special boundary control is required to process the boundary condition in this algorithm. The stimulation and the experimental results show that the new algorithm has strong capability to fit the wavefront and it is more available if the sample points are few. r 2005 Elsevier GmbH. All rights reserved. Keywords: Wavefront reconstruction; Fourier transform; Fried geometry
1. Introduction The use of the discrete Fourier transform (DFT) for wavefront reconstruction from slope map was firstly proposed by Freischlad and Koliopoulos [1]. They applied DFT method to square or rectangular domain, and the method to solve the inverse filters in spatial frequency domain for different geometries were derived in their paper. Their work had been improved and used in interferogram analysis and in adaptive optics systems of astronomy [2–5]. Poyneer et al. suggested the boundary and extension method to fill in missing slope information for large circular apertures with use of DFT method [2]. As for the domain with any arbitrary shape, Roddier did good work in wavefront reconstruction using iterative DFT. Roddier’s algorithm came from Gerchberg’s extrapolation algorithm [6] beyond the Corresponding author. Fax: +8622 23508332.
E-mail address: [email protected] (Z. Wang). 0030-4026/$ - see front matter r 2005 Elsevier GmbH. All rights reserved. doi:10.1016/j.ijleo.2005.06.003
boundaries [3,4]. We name Roddier and Gerchberg’s work with R–G algorithm. In this paper, we introduced both Fried geometry [7] and DFT method to reconstruct the wavefront. And we combined inverse filter and iterative Fourier transforms to give more exact and easy estimate of the wavefront. Our work can be looked as the modification and the improvement of R–G algorithm.
2. Fourier iterative algorithm with Fried geometry A flow chart of our modified new iterative algorithm is shown in Fig. 1. In this chart, the input data is the gradient matrices which can be measured by wavefront sensor, such as the Hartmann–Shack (H–S) sensor. The output of the procedure shown in Fig. 1 is the estimated wavefront with acceptable error. In order to explain the work flow chart, we take the H–S wavefront sensor as
ARTICLE IN PRESS 78
H. Guo, Z. Wang / Optik 117 (2006) 77–81 W(m-1,n-1)
Gradient matrixes with zero value out of the pupil
DFT
Put original gradient matrixes back inside the pupil
Filter for Fried geometry in spatial frequency domain
Recalculated gradient matrixes with Fried geometry
W(m-1,n)
W(m,n-1)
W(m+1,n-1)
W(m,n)
W(m+1,n)
Sx[m,n]
Inverse DFT Error acceptable W(m-1,n+1)
Estimated wavefront
W(m,n+1)
Sy[m,n]
W(m+1,n+1)
Fig. 3. The Fried sensor geometry.
Fig. 1. Flow chart of the new iterative discrete Fourier transform algorithm used to reconstruct a wavefront from the measured wavefront gradients.
have the difference in x- and y-directions S x ¼ rx W ðm; nÞ ¼ 1=2½W ðm þ 1; n þ 1Þ þ W ðm þ 1; nÞ 1=2½W ðm; n þ 1Þ þ W ðm; nÞ,
ð2Þ
S y ¼ ry W ðm; nÞ ¼ 1=2½W ðm þ 1; n þ 1Þ þ W ðm; n þ 1Þ 1=2½W ðm þ 1; nÞ þ W ðm; nÞ. ð3Þ
Fig. 2. Hartmann–Shack grids.
an example. Figs. 2 and 3 show the H–S grids and estimated points of wavefront W(m,n) and the gradients Sx, Sy measured by H–S sensor in Fried sensor geometry, respectively. We want to determine the wavefront W(m,n) over the measurement domain within the pupil from its noisy discrete gradient S, given by Sðm; nÞ ¼ rW ðm; nÞ þ nðm; nÞ with S¼
Sx Sy
! ;
n¼
nx ny
(1)
! ,
where n represents the additive noise to the measured gradient and it will not be discussed in this paper. The first coordinate is defined as the y-coordinate which is the columns in a numeric matrix. And the rows of the numeric matrix are along the x-coordinate. For the lowright sensor locations in Fig. 3 of Fried geometry, we
In the work flow chart that Fig. 1 shows, the gradient matrices are recalculated from the estimated wavefront points with Eqs. (2) and (3). The DFT is applied to both the equations with use of the shift property of the Fourier transform, Eqs. (2) and (3) become 1 i2ðk þ lÞp i2kp ~ lÞ exp DFTfS x g ¼ wðk; þ exp 2 N N i2lp exp 1 , ð4Þ N 1 i2ðk þ lÞp i2lp ~ lÞ exp þ exp DFTfS y g ¼ wðk; 2 N N i2kp exp 1 , ð5Þ N where the number of sample grids is N N discrete signal, k, l are the coordinate in spatial frequency domain, and i2 ¼ 1. We can derive the inverse filter by combine Eqs. (4) and (5). The following formula is the ~ lÞ for Fried geometry: solved wðk; 8 0; k ¼ l ¼ 0; k ¼ l ¼ N=2 > >
i2kp > > < exp N 1 exp i2lp N þ 1 DFTfS x g
~ lÞ ¼ wðk; 1 exp i2kp þ 1 DFTfS y g þ exp i2lp > N N > > >
: 2 kp 2 lp 2 kp 1 else: 8 sin N cos2 lp N þ sin N cos N (6) ~ lÞ will Taking the inverse Fourier transform of wðk; ~ produce the estimated wavefront wðm; nÞ. There are a key assumption when we calculate DFT{Sx} and DFT{Sy} in Eq. (6). It is necessary to
ARTICLE IN PRESS H. Guo, Z. Wang / Optik 117 (2006) 77–81
have difference data consistent with a periodic reconstructed wavefront of period N is x- and y-directions. From Fig. 3, we can conclude that only ðN 1Þ ðN 1Þ array of x and y differences are practically measured by the sensor. In order to calculate DFT during the iterative procedure, we must process reasonably the missing information to obtain the N N array of differences. Furthermore, when the aperture is circular or even arbitrary, the estimated wavefront cross the boundary of the aperture will be biased if no boundary control algorithm is used [2]. Eq. (6) gives the formula of the filter for Fried geometry in spatial frequency domain. This formula will enter the loop shown in Fig. 1. For certain size of matrices, the exponential and the trigonometric operation in Eq. (8) can be calculated before starting the DFT iterative operations. In our new iterative DFT algorithm, boundary control can be done during the course of iteration and
79
has the fewest iterative times to arrive at the minimum reconstruction error. After this, in order to continue the iteration, we put back the measurement data inside the aperture and do the iteration sequentially, until the wavefront estimate is error acceptable. The measurement of current error in the iterative reconstruction process is to calculate the root mean square (RMS) difference. If the simulation is considered, the RMS difference is between the real wavefront and the estimated wavefront 2RR
3 ~ ðx; yÞ2 dx dy 1=2 ½W real ðx; yÞ W 6 7 RR RMS ¼ 4 A 5 , dx dy
(8)
A
where A represents the area of the circular aperture. If real data is considered, the RMS difference is between the measured gradients and the estimated gradients
2RR n 31=2
2
2 o S xreal ðx; yÞ S~ x ðx; yÞ þ Syreal ðx; yÞ S~ y ðx; yÞ dx dy 6A 7 7 . RR RMS ¼ 6 4 5 2 dx dy
(9)
A
it is not necessary to derive additional boundary method. The iterative procedure produces exact approximation cross the boundary of the estimated wavefront. In Fig. 2, the black bold square is ðN 1Þ ðN 1Þ H–S sensor array. The circular boundary tangential with the square is used to indicate those lenslets (also called ðN 1Þ ðN 1Þ) whose outputs are used in the estimation of the wavefront. The biggest square outside the black bold one adds some complementary points and its size is M M ðM4NÞ. In order to accelerate the speed of the computation, M is often set to be the integral power of 2. Before the first iteration, all the measurement gradients are within the circular aperture and the outside domain of it is padded with zero. After the first iteration, we recalculated two ðM 1Þ ðM 1Þ gradient matrices S½ðM 1Þ ðM 1Þ from the first ~ ðm; nÞ by use of Eqs. (2) and (3). estimated wavefront W And we set ~ ðM; nÞ, S x ðM; nÞ ¼ S y ðM; nÞ ¼ W ~ ðn; MÞ S x ðn; MÞ ¼ S y ðn; MÞ ¼ W
ð7Þ
for 1pnpM. Now the new gradient matrices have the same size M M with the original two. In fact, other values can be used to complement the missing row and column, such as zero or constant. At present, we are not so sure about which is the best value for the missing gradients data. The complement method with Eq. (7)
The units of both Eqs. (8) and (9) are the unit of length, i.e. micron or wavelength.
3. Simulation results We develop a simulation that the wavefront is given by the 12th fringe Zernike polynomial (ZP), say, spherical aberration with its coefficient C varied from 0.1 to 6 mm. Its expression is Z ¼ Cð6x4 þ 6y4 þ 12x2 y2 6x2 6y2 þ 1Þ. The circular pupil inside a square sample area is considered and there is no noise added in our work. Fig. 4(a) and (b) show the simulation results. In order to compare with these results, the reconstruction results of R–G algorithm are shown in Fig. 4 too. According to Marechal limit, the RMS fitting error will be acceptable if it is lower than l=14. The wavelength l ¼ 633 nm in our simulation. This limit is also drawn with blue dashed line in Fig. 4. In the case that Fig. 4(a) shows, the RMS errors of both algorithms decrease while the sample points increase. And the new algorithm has lower RMS error than R–G algorithm on every sample level. Even the fitting error of the new algorithm at 30 30 is less than that of R–G algorithm at 100 100. If the sample grids
ARTICLE IN PRESS 80
H. Guo, Z. Wang / Optik 117 (2006) 77–81
Fig. 5. The convergence of the wavefront’s RMS fitting error comes up with the iterative times of Fourier transforms.
to its minimal value within 10 iterative times. The new algorithm always has greater RMS fitting error after the first iteration. If a few more iterations are carried out, new algorithm will give better results than R–G algorithm. For example, it requires about 15 iterations to get to the equal fitting error of R–G algorithm if the sample grids is 100 100. In fact, we also found the similar iterative cases in other spherical aberration coefficients as shown in Fig. 5.
Fig. 4. Simulation results with the new iterative algorithm (blue color symbols) and R–G algorithm (red color symbols) for spherical aberration 0.1–6 mm: (a) 200 iterative times and (b) 10 iterative times.
is 10 10 or lower, R–G algorithm produces obvious fitting error if the coefficient is more than 0.2 mm. And the new algorithm produces the same error when the coefficient gets to 1.4 mm. In addition, the new algorithm gives results that can be totally acceptable when the sample grids is 30 30. At the same time, the results of R–G algorithm can be acceptable only if the coefficient is less than 2.5 mm. The iterative DFT times of Fig. 4(a) is 200, which perhaps will be not suitable for real-time system. Then we just restrict the iterative times to 10 in Fig. 4(b). In this figure, the new algorithm is no better than R–G algorithm if the sample grids is 100 100. Even the RMS error under this sample is bigger than that of 60 60 and almost equal to 30 30 using the new algorithm. It is clear that iterative times will influence the results of the algorithm. The iterative times of different sample density for 1 mm spherical aberration coefficient is shown in Fig. 5. In this figure, we can see that the R–G algorithm converges more quickly than the new algorithm and gets
4. Results of real data tests In order to observe the behavior of the new algorithm, we use the Hartmann–Shack (H–S) wavefront sensor to measure the wavefront aberration of the eye. With use of the experimental setup that paper [8] showed, one human eye was measured and 75 facula array pictures were gotten within 3 s. We analyzed 25 pictures which were obtained in the 2nd second of the measurement. The sample points were 18 18 and the pupil size was 5.6 mm.15, 35, and 65 ZPs were also used, respectively, as the basis functions to fit the measured wavefront. The iterative times of both the R–G algorithm and the new algorithm were confined to 30 times. The mean RMS errors of five processing methods are shown in Fig. 6. It is clear that the RMS difference between the gradients measured and gradients estimated has minimum value using the new iterative DFT algorithm. It decreases about 33% compare with 65 ZPs fitting and about 25% compare with R–G algorithm.
5. Conclusions The new iterative DFT algorithm developed here has the following features. The first feature inherits from the Roddier–Gerchberg algorithm, which can used to
ARTICLE IN PRESS H. Guo, Z. Wang / Optik 117 (2006) 77–81
81
No. 60438030), the Chinese Ministry of Education’s Nankai University and Tianjin University Cooperation Foundation and the Key Research Foundation of Scientific and Technical Committee of Tianjin City of China (No. 033183711).
References Fig. 6. RMS error for measuring the wavefront of the eye with H–S sensor.
analyze the wavefront in circular or the arbitrary pupil. Secondly, no special boundary control algorithm is needed in this algorithm for circular or arbitrary pupil. The production of periodic signal required by Fourier transform is solved with Gershberg’s extrapolation method, which will give accurate approximation of boundary differences during the iterative procedure. Thirdly, without the noise, the new algorithm can reconstruct the wavefront more accurately than its previous algorithm, especially in the case of low number of the sample points. If the sample density is high, more iterative times will be required. At last, the measured orthogonal slopes centered in-between the estimated points on the wavefront, which is determined by Fried geometry.
Acknowledgments This work was supported by the National Nature Science Foundation of China (Key Research Project
[1] K.R. Freischlad, C.L. Koliopoulos, Modal estimation of a wave front from difference measurements using the discrete Fourier transform, J. Opt. Soc. Am. A 3 (11) (1986) 1852–1861. [2] L.A. Poyneer, D.T. Gavel, J.M. Brase, Fast wave-front reconstruction in large adaptive optics systems with use of the Fourier transform, J. Opt. Soc. Am. A 19 (10) (2002) 2100–2111. [3] C. Roddier, F. Roddier, Interferogram analysis using Fourier transform techniques, Appl. Opt. 26 (9) (1987) 1668–1673. [4] F. Roddier, C. Roddier, Wavefront reconstruction using iterative Fourier transforms, Appl. Opt. 30 (11) (1991) 1325–1327. [5] A. Dubra, C. Paterson, C. Daninty, Wave-front reconstruction from shear phase maps by use of the discrete Fourier transform, Appl. Opt. 43 (5) (2004) 1108–1113. [6] R.W. Gerchberg, Super-resolution through error energy reduction, Opt. Acta 21 (1974) 709–720. [7] D.L. Fried, Least-Square fitting a wave-front distortion estimate to an array of phase-difference measurements, J. Opt. Soc. Am. 67 (1977) 370–375. [8] H. Guo, Z. Wang, Q. Zhao, W. Quan, Y. Wang, Individual eye model based on wavefront aberration, Optik 116 (2) (2005) 80–85.