A theoretical model for reconstruction errors in the Hartmann-Shack wave-front sensor was constructed. The measured pattern of a regular grid of spots can be analyzed by their individual centroiding or by a global Fourier demodulation. We investigate pixelization errors in the camera and Poisson errors in the camera pixels. We show that by the Fourier demodulation technique it is possible to overcome pixelization errors which occur in the traditional centroid technique. By supporting simulations we show that the two methods coincide for an infinite number of pixels per spot.
© 2006 Optical Society of America
In the Hartmann-Shack sensor , photons pass through an array of lenslets, focusing corresponding sections of the incoming wave-front on an array of pixels, usually a CCD (coupled charge device) camera. This picture of the foci of the dissected wave-front is then sampled to yield the wave-front slopes. There are two methods for finding the sampled wave-front slopes. The first and most common method is the centroid approach [2, 3, 8]. In this method we regard each focus as an individual light spot. We then measure the displacement of the actual incoming light spot from the physical center of the lenslet. The physical center represents the place the ray would have been focused at, had it been perpendicular to the lens. This displacement provides us with the local slope: the necessary information for the reconstruction of the phase of the incoming wave, or for closing an adaptive optics phase correction loop.
The second method for calculating the spot pattern uses a two-dimensional Fourier demodulation . In this technique we look at the picture as a whole, rather than investigating each focus spot by itself. When the grid of spots is cartesian, we regard the small displacements from regularity as small phase aberrations, which are measurable directly in the Fourier plane. Wave-front reconstruction can be performed directly in that plane [5,6].
2. Centroid error
2.1. Pixelization error
When using a CCD or another pixelated array, we cannot regard the world as continuous, but limit ourselves to a sampled image of the incoming beam. Thus we suffer from aliasing effects - the center of the discrete pattern does not coincide with the center of the original (continuous power distributed) beam. In order to construct our model we introduce some fundamental assumptions:
- The intensity spread (lenslet diffraction spot) is a two dimensional Gaussian random vector, with independent elements I = (2πσ2)-1 exp[- (x 2 + y 2)/2σ2].
- The total area assigned to one spot (lenslet shadow) of light by the camera is D × D, divided into N × N pixels (the whole detector holds M × M pixels). In addition, let us define one pixel’s size as d = D/N. We also define D = 2R, ρ =R/σ, and μ = d/σ, the lenslet and pixel sizes relative to the spot size.
- The axis origin is in the center of the sub-aperture of the respective lenslet.
- The signal in every pixel is the accumulation of the intensity in the area of that pixel. Let us define the signal in the element (n, m) of the CCD array as
where T = [Q(-ρ)-Q(ρ)]-2 is the normalization factor and Q(x) is defined through the Gaussian error function (erfc)
We can now continue and write the pixel signal as
and the separability of Fnm above is permitted due to the separability of I(x,y). The real position of the center of mass is now calculated using the discrete power distributions. Let us denote the measured position as x̂, the measured pixel as n,m, and the real (continuous) position as x. Then we can write the expectation values E[ ] = d〈(n,m)〉 = d(N/2,N/2) , where 〈(n,m)〉 is the average of the position vector (n,m) in pixel units.
Hence 〈(n,m)〉 = (ΣnΣm nFnm,ΣnΣm mFnm) , and the summation is over the whole set of pixels in the CCD array. We now continue with the second moment calculation. Without loosing generality, let us look at one coordinate only (due to symmetry 〈m 2〉 = 〈n 2〉).
Let us now continue and investigate the error of the discrete result due to the pixelization. First let us define the error here as the distance between the discrete centroid and the continuous centroid, or
2.2. Photons quantization error
Due to the dual characteristics of light, we can regard the incoming power as the summation of a discrete number of incident particles . Let us now look at a single picture shot for one of the pixels, during one intensity integration time. Within that interval of time, K photons hit this pixel, each independent of the others. Therefore we have K independent spatially sampled integrated Gaussian random variables.
Again we would like to calculate the average error of the centroid in this case. Let us denote by Xi the displacement of the ith photon from the origin, and X (see Fig. 1) as the average displacement (the centroid), namely X = K -1 ΣK i=1 Xi. Since E[Xi] = 0, we also have E[X] = 0 and
Again we can calculate the total error in the two axes by the root summation of the squares of each, hence
where X̂ represents a vector in the XY plane, and X a scalar (one of the axes, arbitrarily). We will now apply the smoothing theorem on σ2 xi , while recalling that the arrival of photons per spot is a Poisson random variable with an average number of photons per spot λ.
To simplify this function let us differentiate and integrate the following quantity:
3. Fourier demodulation error
3.1. Pixelization error
We will now investigate the error in the Fourier demodulation method, where the regular, cartesian spot pattern is treated as a whole, rather than each single spot individually [4, 5, 6]. As explained in , we model the intensity pattern of the light on the focal plane as a function whose period is k, represented by the first harmony:
where F is the focal length of the lenslets and ϕx,ϕy are the wave-front gradient components. V(x,y) is the aperture function holding the Hartmann-Shack pattern. As before, we will look at the power integral due to the finite size of one CCD pixel. Now, however, we will integrate over the whole pattern and not over one spot only. Let as denote the number of pixels of the whole CCD array as M × M, and its physical size as L × L = Md × Md. For convenience we also define H = (L - d)/2, and W as the normalization factor. The pixelated light intensity is now
We express this term as a sampled convolution (the subscript D represents the discreteness):
where we make use of the Λ or Rect function
Let us substitute x ⇒ x-H,y ⇒ y-H. Hence
We will now demodulate the pattern in order to derive the x gradient component ϕx .
3.1.1. Fourier transform
We start by calculating the Fourier transform ĨD from the discrete pattern ID. Since ID is a product of 3 elements, its Fourier transform is the convolution of the Fourier transform of each,
The first of the three convolved components is the transform of the lenslet sub-aperture
where sinc(x) = sin(x)/x. Its frequency support is (-1/L,1/L). The second component is the signal filtered with a low pass filter
Its frequency support is (-1/D - q 0/2π,1/D-q 0/2π). The last component is a bed of nails which is the result of the sampling. It has δ elements spread 1/d from each other,
When the number of lenslets is large (M ≫ N), we can assume A ≈ δ(u,v). Also, examining the convolution B∗C, we notice that the support of B is smaller than the distance between the δ functions of C for N > 2, or 1/d > 2/D + q 0/π. Thus, there are no aliasing effects due to the sampling. Here we have assumed a rectangular aperture, whereas the more common case is a circle. Since we are only interested in the width of the aperture relative to the rest of the functions in the Fourier space, it does not make a difference which of the two we choose, since both the sinc and the airy functions are much narrower than the rest of the functions in our system.
We will now shift the whole function so that the first x side lobe of ÎD will be at the origin, and then apply a low-pass filter with q 0/2π cutoff frequency. Let us denote the original side lobe of I as Cx = V(r) exp(-iFϕx), whose transform is Cx . The shifted, smoothed result is
Applying the inverse Fourier transform yields:
The analyzed array size M holds (or is extended mathematically to hold) an integer number of lenslets,M = kN. Thus exp[-iπ(M- 1)/N] = exp(iπ/N), and
Using the adiabatic property of Cx , F∇2ϕ(x) ≪ ∇(2πx/D), the wave length of the gradient variations is λϕx ≫ D > d. Thus we can take Cx out of the integral and write:
3.1.3. Phase extraction
We will now calculate the phase of the wave-front directly. If we disregard the pixelization effect, the phase (the argument) of the complex K(x,y) should be the phase of the wave-front we are looking for, e.g. ϕx(x,y). Thus the error of the pixelization is ε = |Fϕx - ∠K|. However, as we have shown, the error in this case is not statistical, but deterministic. This allows us to fully recover ϕx from K with no pixelization errors at all. This should not surprise us, since we are sampling beyond Nyquist’s frequency (our adiabatic assumption). When ε = 0 the correction for ϕx can be written as
3.2. Photon quantization error
In the Fourier model (Eq. 14) a sub-aperture spot is approximated by one period of a cosine function. Thus the probability function for each photon per one spot can be written as P(x) = (kx/2n)[1+cos(kx ∙x)]Λ(-π/kx, π/kx). Let us assume that the ith spot receives Ki photons, and the average sampled intensity is a cosine function. Under this model, shifts in the spot position are equivalent to a change in the cosine phase, to which we shall refer as θ.
In order to determine the phase error, we will investigate the errors due to the different displacements from the cosine equilibrium, and then see how they affect the phase. Let us denote by Xi the displacement of the ith spot from the center, and by xi,j the displacement of the jth photon in the ith spot from the center. As in Eq. 7 we get E |Xi] = 0
where σxi (the displacement error) was calculated by:
We express the phase error as a function of the position fluctuation. Since the two axes are assumed to be independent, we simply need to multiply the error by √2 to make the expression two dimensional. Using simple geometry, the phase of the cosine function can be expressed as θi = √2kxXi. Thus, its mean is E[θi] = 0, and its variance (the phase error) can be calculated:
Let us now denote the average number of photons per sub-aperture as λ, and apply the smoothing theorem on σθi:
Using Eq. 12 we can write the total error in X due to the photon quantization as
(The displacement error is achieved simply by a dividing the above σθ by kx).
4. Comparison between the different methods
To compare the two demodulation techniques we need to match their probability characteristics. We will set the mean value and the variance of the two methods to the same values. The mean value had been considered to be zero throughout, and thus no further adjustments are required. For the variance, we will use the one calculated for the cosine (from the Fourier method, Eq. 28) as the variance of the Gaussian (from the centroid method). This will allow us to compare the two position errors derived before (Eqs. 13 and 32). The assumption that the beams have a wider profile can be justified when the source is large, the quality of the lenslets is low, or when integrating over a fluctuating wave front.
Fig. 2 clearly shows the fundamental result expected by the cancellation of the pixelization error for the Fourier method: as the number of pixels per spot decreases, the difference between the two errors drops dramatically, in favor of the Fourier demodulation method. For infinite number of pixels per spot, as expected, the two methods coincide and result with the same error.
4.1. Readout noise
Another source of noise in CCD arrays is the read out noise, affecting the number of electrons read in the output amplifier. In order to provide a feeling for its effect, we used a very simple model, where we take it as an addition of “fake photons” to every pixel. This avoids the necessity to include occasional negative noise, but a gaussian model yields similar results. The read electron distribution is regarded as an addition of a Poisson random variable λRoN with an average of RN photon hits per one pixel (thus N 2 RN per one spot). Let us denote the average displacement due to the read out noise only by XRoN. We assume that the displacement of one electron due to read out noise can be modeled as a uniform random variable (since it does not depend on the location of the pixel in the CCD). Thus the displacement is a uniform random variable, while the number of photons is spread by a Poisson random variable. Let us denote by Xγ and λγ the real photons displacement and amount respectively. To ease our calculations, we define XTOT as the weighted average of both real and fake photons in one spot. We calculate XγwR - the error in Xγ when taking into account the readout noise as well, using XTOT :
The problem of the statistics of division of Poisson-dependent processes was tackled early on in the framework of stellar interferometry (,  & ) and imaging ( & ). To estimate the fringe visibility, differences of photon events are divided by their sums, but the noise in the visibility is difficult to analyse. We adopt the common approach in these estimators, where the sample-to-sample variation in the denominator is replaced by its ensemble average, which is valid when λ > 1. (Note the factor of 2 for the variance in two dimensions). For each demodulation technique, σXγ should be replaced with its own displacement variance calculated above (Eq. 13 and Eq. 32).XRoN is calculated by substituting uniform variable variance instead of σxi. in Eq. 12. Hence σ2 XRoN = 2D 2/12exp(- N 2 RN) [Ei(N 2 RN) - lnN 2 RN]
The final errors according to this simple model are shown in Fig. 2. A more complex model is now under construction.
5. Summary and conclusions
We have shown that unlike the centroid method, the Fourier method does not suffer from pix-elization error. The basic principle behind this phenomenon is the fact that our sampling is performed beyond the Nyquist frequency, and so no aliasing effects occur. We wish to emphasize that the Nyquist criterion was preserved under the assumption of the cosine shape of the incoming photons while using the Fourier demodulation methods. In other words, spreading the photons over more pixels results in better detection efficiency. Hence lower quality optics may be employed for the lenslet array. For optimal results, the point spread function of the foci should be designed to be sinusoidal. In a more cumbersome computation, additional harmonics could be used in the demodulation. On the other hand, having more pixels can contribute to more read out noise, so the two should be balanced according to the application.
Parts of this project were funded in conjunction with the Israeli Ministry of Science and the European Research and Training Network SHARP EYE.
References and links
01. I. Ghozeil, in: Optical Shop Testing, D. Malacara , ed. (Wiley, New York, 1978) Chap. 10.
02. R. K. Tyson, Principles Of Adaptive Optics, Second ed. (Academic, New York, 1998).
03. R. K. Tyson, ed., Adaptive Optics Engineering Handbook (Marcel Decker, New York, 2000).
04. Y. Carmon and E. N. Ribak: “Phase retrieval by demodulation of a Hartmann-Shack sensor,” Opt. Comm. 215, 285–288 (2003). [CrossRef]
05. A. Talmi and E. N. Ribak: “Direct demodulation of Hartmann-Shack patterns,” J. Opt. Soc. Am. A 21, 632–639 (2004). [CrossRef]
06. Y. Carmon and E. N. Ribak: “Fast Fourier demodulation,” Appl. Phys. Lett. 84, 4656–4657 (2004). [CrossRef]
07. J. W. Goodman, Statistical Optics (Wiley-Interscience, New York, 1985).
08. J. W. Hardy, Adaptive Optics for Astronomical Telescopes, (Oxford University Press, 1998).
09. W. J. Tango and R. Q. Twiss, “Michelson stellar interferometry,” in Progress in Optics XVII, E. Wolf, ed. (North Holland, 1980) 239–277.
11. M. M. Colavita: “Fringe Visibility Estimators for the Palomar Testbed Interferometer,” Pub. Astron. Soc. Pac. , 111, 111–117 (1999). [CrossRef]
12. J. A. Zadnik and J. W. Beletic: “Effect of CCD readout noise in astronomical speckle imaging,” Appl. Opt. 37361–368 (1998). [CrossRef]
13. S. B. Howell, B. Koehn, E. Bowell, and M. Hoffman: “Detection and Measurement of Poorly Sampled Point Sources Imaged With 2-D Array,” Astron. J. 112, 1302–3611 (1996). [CrossRef]