2038
J. Opt. Soc. Am. A / Vol. 27, No. 9 / September 2010
Jia et al.
Minimum variance unbiased subpixel centroid estimation of point image limited by photon shot noise Hui Jia,* Jiankun Yang, and Xiujian Li Tech-Physical Research Center of Science College, National University of Defense Technology, Changsha, Hunan 410073, China *Corresponding author:
[email protected] Received January 25, 2010; revised June 24, 2010; accepted June 30, 2010; posted July 6, 2010 (Doc. ID 123287); published August 18, 2010 An unbiased subpixel centroid estimation algorithm of point image is proposed through the compensation of the systematic error of the center of mass method. The Cramér–Rao lower bound on centroid estimation variances is derived under the photon shot noise condition and is utilized to evaluate the proposed algorithm. Numerical analysis shows that the proposed centroid estimator attains the required lower bound; thus the proposed algorithm can be asserted as a minimum variance estimator. Simulation results indicate that the centroid accuracy is maximized when the Gaussian width of the signal spot is 0.2–0.3 pixel and the estimator can attain subpixel accuracy close to 1/100 pixel when 1000 photons are detected. © 2010 Optical Society of America OCIS codes: 040.0040, 120.1880, 150.1488, 100.6640.
1. INTRODUCTION Centroid estimation of the point source target is widely applied in many fields, such as star trackers that estimate the attitude of spacecraft utilizing the location of stars [1–4], and Shack–Hartmann wavefront sensors that detect the wavefront through the position of spots generated by lens arrays [5–7]. Other applications include biological particle tracking using optical microscopy [8] and edge detection of machine vision [9,10]. In these applications, the imaging systems are always slightly defocused in order to expand the target spot to cover several pixels; thus the accuracy can attain subpixel accuracy through a centroid estimation algorithm [11]. It is meaningful to find the minimum variance unbiased estimation (MVUE) algorithms which can attain the accuracy performance bound that is mainly influenced by the systematic error and the noise of the imaging sensor. The cause of systematic error is the calculation bias of the estimation algorithm and the insufficient sampling of the spot image. In addition it is difficult to find an unbiased estimation algorithm, because the detected signal of each pixel of a charge-coupled device (CCD) cannot be expressed analytically in most situations. Fisher and Naidu studied the systematic error of different algorithms under the condition of first-order approximation of a Gaussian expression [10]. They evaluated the Gaussian estimator, the center of mass (CoM) algorithm, the linear interpolation algorithm, the parabolic estimator, and various other methods. They concluded that only the Gaussian estimator algorithm is unbiased. However, even though the signal intensity profile through an optical image system can be assumed to be a Gaussian distribution, the distribution of signal detected by each pixel will not be Gaussian after each pixel integrates light over its field of view. Thus 1084-7529/10/092038-8/$15.00
their assumption cannot be realized in a practical system. Brian and Alexander gave a general expression of systematical error for the CoM method through Fourier spectrum analysis [12]. Numerical simulation of the systematical error of different weighted CoM algorithms is presented by H. Li et al. [5]. The systematic errors under different Gaussian widths, CCD fill ratios, and detector sizes were studied numerically and experimentally in [13,14]. This paper furthers the Fourier spectrum analysis of systematical error of the CoM algorithm under the Gaussian signal intensity profile assumption. The analytical expression of systematical error under different Gaussian widths is given which is consistent with the numerical simulation. An unbiased centroid estimation algorithm is proposed through the compensation of the systematic error of the CoM method, which can be known as the compensated center of mass (CCoM) algorithm. The noise of an imaging sensor such as a CCD and a complementary metal-oxide semiconductor is another main factor restricting the accuracy performance of the centroid estimation. The noises mainly include photon shot noise (PSN), background noise, dark current noise, and readout noise. This paper concentrates on PSN, because PSN is the ultimate limitation on the accuracy performance of the centroid estimation. Other noises can be reduced through the developments of technology. The background noise can be suppressed by baffle [15]. The dark current noise can be reduced by lower temperature [16]. The readout noise of the recently developed electron multiplying CCD (EMCCD) can be reduced to halfelectron per frame [17]. It is therefore reasonable to discuss the bound of the accuracy performance of centroid estimation under the shot noise limit, especially when a faint object is imaged. © 2010 Optical Society of America
Jia et al.
Vol. 27, No. 9 / September 2010 / J. Opt. Soc. Am. A
The proposed CCoM algorithm is evaluated by comparing the accuracy with the Cramér–Rao lower bound (CRLB). The CRLB is a lower limit on the variance of any unbiased estimated method in the classic statistic theory. Winick derived the CRLB of the variance of centroid estimation in two dimensions under the photo noise assumption [1]. The problem of calculating the position of a Gaussian pulse in Poisson noise when there is no background noise was considered by Rye and Hardesty [18]. Chen and Rao calculated the estimation ratio of the CoM algorithm through the comparison of the CRLB [6], which indicated that the (CoM) algorithm cannot attain the CRLB well. The CRLB of the variance of centroid estimation under the PSN limit is derived directly from statistic theory in this paper. Numerical analysis shows that the proposed centroid estimator attains the lower bound; thus the proposed algorithm can be asserted as a minimum variance estimator. We propose an unbiased subpixel accuracy centroid estimation algorithm which compensates the systematical error of the CoM method. The systematical error is analyzed and our algorithm is proposed in Section 2. The CRLB under the PSN assumption is derived in Section 3. Section 4 gives the accuracy analysis of our method under the PSN, and the comparison with the CRLB is presented which shows that the proposed method is a kind of MVUE. Practical application considerations are discussed in Section 5. Finally, the conclusion is given in Section 6.
2. ANALYSIS OF SYSTEMATIC ERROR AND COMPENSATED ESTIMATION ALGORITHM A. Systematic Error of CoM Algorithm The CoM algorithm is the simplest and most widely used method to calculate the centroid position of a symmetric image spot. The estimated location of the centroid is given by
兺xU i
xˆ =
i,j
兺U
兺yU
ij
j
,
yˆ =
ij
i,j
ij
i,j
兺U
共1兲
,
I共x,x0兲 =
冋
1
冑2
exp
− 共x − x0兲2 22
册
.
2039
共3兲
For convenience, we measure all distances in units of the pixel pitch, and the fill factor of the CCD is 100%. We assume that each pixel has same photon response, and the center of the ith pixel is xi. Then the average number of photoelectrons generated by this pixel from the image spot is Si = Nphgi共x0兲 =
冕
xi+0.5
I共x,x0兲dx,
共4兲
xi−0.5
where Nph is the average number of photoelectrons generated from the image spot. Then, the systematic error of the CoM algorithm can be given as
兺xU i
␦x = xˆ − x0 =
i
兺U i
兺xg
i
i i
− x0 = i
i
兺g
− x0 .
共5兲
i
i
From Eqs. (3)–(5), the systematic error ␦x is related to the real centroid x0 and the Gaussian width . The numerical simulation results of the relationship between ␦x and x0 under different ’s are shown in Fig. 1. It can be seen that there is an approximately sine relationship between ␦x and x0 under different ’s. We can also conclude that when the systematic error decreases the Gaussian width increases. Extending from that observation it will be beneficial to find the analytical relationship of the systematic error ␦x between the real centroid x0 and Gaussian width . The pixel signal function can be written as f共x兲 = I共x,x0兲 丢 rect共x兲,
共6兲
where rect共x兲 is a rectangle function denoting the pixel response and 丢 denotes convolution operation. Then f共x兲 is multiplied by the sampling function comb共x / T兲,
ij
i,j
where 共xi , yj兲 is the coordinate of each pixel and Uij is the detected signal at pixel 共xi , yj兲. We assumed that the signal power is Gaussian distributed, because a Gaussian function is a reasonable approximation to illustrate the signal blurring effect in an image system. The signal intensity profile can be written as I共x,y,x0,y0兲 =
1 2
2
冋
exp
− 共x − x0兲2 2
2
册 冋 exp
− 共y − y0兲2 22
册
, 共2兲
where 共x0 , y0兲 is the true centroid position and is the standard deviation of the two-dimensional (2-D) Gaussian profile which is also the Gaussian width. We will discuss the one-dimensional (1-D) situation in the following, as the 2-D situation can be analyzed similarly. For the 1-D case, Eq. (2) reduces to
Fig. 1. Systematic error of CoM algorithm under different Gaussian widths.
2040
J. Opt. Soc. Am. A / Vol. 27, No. 9 / September 2010
Jia et al.
共7兲
g共x兲 = f共x兲comb共x/T兲,
where T denotes the sample period which is equal to 1 pixel in our assumption. Then Eq. (5) can be written as [12]
冕 冕
⬁
xg共x兲dx
−⬁
␦x =
− x0 = −
⬁
g共x兲dx
rithm that can gain additional precision is our CCoM algorithm. From Eqs. (5) and (12), we can get this following expression: xˆ = x0 + ␦x = x0 −
1
⬁
兺 共− 1兲 exp关− 2共n兲 兴sin共2x n兲/n. n
2
共15兲
G⬘共0兲 2iG共0兲
共8兲
− x0 ,
Let the estimated centroid be xˆc, and then the estimator of CCoM can be further reformulated to
−⬁
where G共s兲 denotes the Fourier transform of the sampled pixel signal function g共x兲. Brian and Alexander derived the systematic error expression through Fourier analysis approach [12],
兺 F⬘共n/T兲sin共2x n/T兲 0
e
冋
n=1
⬁
Fe共0兲 +
兺 F 共n/T兲2 cos共2x n/T兲 e
0
n=1
册
,
共9兲
where Fe共s兲 denotes the Fourier transform of the zero centered pixel signal function fe共x兲 and T denotes the sample period which is equal to 1 pixel in our assumption. Then, from Eqs. (3) and (6), we can get Fe共s兲 = F兵fe共x兲其 = F兵I共x,0兲 丢 rect共x兲其 = exp关− 2共s兲2兴sinc共s兲.
共10兲
So Fe共0兲 = 1, and Fe共n兲 = 0. Then F⬘e 共n兲 = exp关− 2共n兲2兴cos共n兲/n = 共− 1兲nexp关− 2共n兲2兴/n. 共11兲 Therefore
␦x =
1
兺
⬁ n=1
共− 1兲nexp关− 2共n兲2兴sin共2x0n兲/n. 共12兲
Furthermore when ⬎ 0.2 this means exp关−2共兲2兴 Ⰷ exp关−2共n兲2兴, where 共n ⬎ 1兲 we can get
␦x ⬇
1
xˆc −
1
⬁
兺 共− 1兲 exp关− 2共n兲 兴sin共2xˆ n兲/n =
n=1
兺xU i
n
2
c
i
兺U
i
. i
i
共16兲
⬁
␦x =
0
n=1
exp关− 2共兲2兴sin共2x0兲.
We can get the estimated centroid xˆc through solving Eq. (16) using a numerical iteration approach and sufficient accuracy can be achieved when n = 3. The Gaussian width is determined by the image system, and it can be calibrated through appropriate approaches. The accuracy of the CCoM algorithm without noise is evaluated. Note that the center of the image spot x0 appears equally anywhere of the pixel; the root-meansquare (RMS) systematic error of CoM and the RMS residual error of CCoM are calculated assuming that x0 is distributed equally in the range [⫺0.5, 0.5]. Figure 2 shows the RMS error of the CoM algorithm and the RMS error of the CCoM algorithm under different Gaussian widths with the centroids calculated over 5 pixels window. It can be seen from the figure that the residual error is substantially reduced after the compensation. However, the RMS error increases rapidly when the Gaussian width is too small and too large. The reason is that when the Gaussian width is extremely small, the total power is concentrated on 1 pixel which means that the detected signals of each pixel are not sensitive to the movement of the spot image. This phenomenon is consistent with the result shown in Eq. (14) and can be seen clearly in Fig. 3. When the Gaussian width is extremely large, the spot sig-
共13兲
Therefore as shown in Eq. (13) the systematic error ␦x has a sine relationship with real centroid x0 and decreases quickly when the Gaussian width increases which is consistent with the numerical simulation results. When → 0, exp关−2共n兲2兴 → 1, we can get
␦x →
1
兺
⬁ n=1
共− 1兲nsin共2x0n兲/n → − x0 .
共14兲
From Eqs. (5) and (14), we can see that when the Gaussian width is too small, the estimated centroid xˆ will approach zero. B. CCoM Algorithm If we correct the systematic error presented in Eq. (12), better accuracy performance can be expected. One algo-
Fig. 2. RMS error of CoM and CCoM algorithms under different Gaussian widths.
Jia et al.
Vol. 27, No. 9 / September 2010 / J. Opt. Soc. Am. A
2041
兺 g 共x 兲 = 1. i
共20兲
0
i
Therefore
ln p共c;x0兲 x0
=
兺
cigi⬘共x0兲
i
gi共x0兲
共21兲
,
where gi⬘共x0兲 =
Fig. 3. Signal intensity profile and pixel detected signal under different real spot positions x0 and Gaussian widths (line: signal intensity profile; bar: pixel detected signal).
nal cannot be adequately detected by the finite pixels window (5 pixels in this example). The effects of the systematic error due to the finite pixels window are not considered in this paper because the centroid estimation error due to noise increases as the Gaussian width increases, which will be presented in Section 4, and the centroid estimation error due to noise will then be the main component of centroid estimation error instead of systematic error.
gi共x0兲 x0
共22兲
,
then
2 ln p共c;x0兲
兺
=−
x02
ci关gi⬘共x0兲兴2 gi2共x0兲
i
+
兺
cigi⬙共x0兲 gi共x0兲
i
共23兲
.
From Eq. (23), it follows that E
冋
2 ln p共c;x0兲 x02
册
=−
兺
关gi⬘共x0兲兴2 gi2共x0兲
i
E关ci兴 +
gi⬙共x0兲
兺 g 共x 兲 E关c 兴. i
i
i
0
共24兲 Because ci is Poisson distributed with mean Nphgi共x0兲, we can get E关ci兴 = Nphgi共x0兲.
3. CRLB OF CENTROID ESTIMATION LIMITED BY PHOTON SHOT NOISE
Substituting Eq. (25) into Eq. (24) yields
The photocurrent would be unstable due to the fluctuation of photon detection when the CCD is illuminated by the image spot. This fluctuation is known as PSN, and can be described through Poisson process [18,19]. In this section, we will discuss the CRLB of centroid estimation limited by PSN. The probability density of the output signal ci of the ith pixel can be written as p共ci ;x0兲 = exp关− Nphgi共x0兲兴
关Nphgi共x0兲兴 c i!
−E
冋
1
2 ln p共x; 兲 2
册
.
共17兲
共18兲
0
i
=−
兺N i
phgi共x0兲
−
兺 ln c + 兺 c i
i
i
ln关Nphgi共x0兲兴.
i
共19兲 From Eqs. (3) and (4), we have
x02
册
= − Nph
兺
关gi⬘共x0兲兴2
i
gi共x0兲
+ Nph
兺 g⬙共x 兲. i
0
i
共26兲 From Eq. (20), we have
兺 g⬙共x 兲 = 0, i
共27兲
0
and therefore through expressions (18), (26), and (27) it follows that var共xˆ0兲 ⱖ −E
冋
1
ln p共c;x0兲 2
x02
册
1 = Nph
兺 i
关gi⬘共x0兲兴2
.
gi共x0兲 共28兲
,
兿 p共c ;x 兲 i
冋
2 ln p共c;x0兲
i
where E denotes the expectation operator. It follows from Eq. (17) that ln p共c;x0兲 = ln
E
ci
We can find the accuracy bound of the centroid estimator through the discussion of the CRLB. The CRLB is the low bound of the variance of any unbiased estimator [20]. The variance of any unbiased estimator must satisfy the following expression [20]: var共ˆ 兲 ⱖ
共25兲
The above equation represents the CRLB for an unbiased estimator of x0 under the PSN, which is consistent with Winick’s result [1]. This equation indicates that the low bound of the variance of an unbiased centroid estimator is in inverse proportion to the detected photon numbers Nph under the photon noise assumption. It is also related to the true centroid position x0 and the incident signal power profile. And a result with lower variance can be estimated if the pixel signal gi共x0兲 is more sensitive to the true centroid position x0. As mentioned in Section 2, Eq. (28) can be averaged over the selected pixel regions, thus yielding a low bound on the average mean-squared error. The low bound of the RMS centroid error can be defined as
2042
J. Opt. Soc. Am. A / Vol. 27, No. 9 / September 2010
␦CRLB =
冋冕
0.5
var共xˆ0兲dx0
−0.5
册
Jia et al.
1/2
.
共29兲
4. ACCURACY ANALYSIS OF THE CCoM ALGORITHM In this section, we will numerically analyze the accuracy performance of the CCoM algorithm. Under the PSN assumption, the signal detected by each pixel can be simulated as Ui = R共Nphgi共x0兲兲,
共30兲
where R is random number of Poisson distribution with mean Nphgi共x0兲. The RMS centroid error of the CoM algorithm is calculated from Eqs. (5) and (30) assuming that the centroid x0 has a uniform distribution over 1 pixel. The CRLB of the RMS centroid error is calculated from Eq. (29). In Fig. 4, the results under different detected photon numbers are presented along with the RMS sys-
tematic error of the CoM algorithm calculated in Section 2. The following conclusions can be drawn from Fig. 4: 1. The CRLB of the variance of unbiased centroid estimation increases when the Gaussian width of the spot image increases. The reason is that the signal to noise ratio is reduced when the Gaussian width expands [13]. Furthermore, the CRLB increases dramatically as the Gaussian width decreases to a small value, as discussed in Section 2. 2. The RMS centroid error of CoM is strictly constrained by the systematic error, the error due to noise bounded by the CRLB. The error of the CoM algorithm achieves the CRLB when the Gaussian width is between 0.4 and 0.8 pixel. Through simulation analysis the lowest RMS error is about 1/50 pixel with 1000 incident photons. Additionally the CRLB increases when the detected photon number decreases, as evident through Eq. (28). When 100 photons are detected, the lowest RMS centroid error of CoM is approximately 1/20 pixel. 3. The CRLB decreases when the Gaussian width shrinks and achieves the minimum RMS error as the Gaussian width is around 0.2 pixel. The CoM algorithm fails to attain the CRLB when the Gaussian width is smaller than 0.4 pixel as a result of systematic error. Thus, if systematic error is compensated, the estimation algorithm can achieve better accuracy performance. The RMS centroid error of the CCoM algorithm is calculated similarly from Eqs. (5), (16), and (30). The results under different detected photon numbers are shown in Fig. 5. The figure also includes the CRLB and the RMS centroid error of the CoM algorithm for comparison. Figure 5 indicates that, compared with the CoM algorithm, the CCoM algorithm attains the CRLB for different Gaussian widths and can be asserted as a kind of minimum variance unbiased estimator. The CCoM algorithm attains the highest accuracy of the CRLB close to 1/100 pixel when the Gaussian width is about 0.2–0.3 pixel with 1000 detected photons, which is a significant improvement compared with the CoM algorithm. When the detected photon number decreases to 100, the lowest error is about 1/25 pixel, which still proves to be superior compared to the CoM algorithm. Consequently, the imaging system should be designed to make sure the Gaussian width of target spot is around 0.2–0.3 pixel which is different from previous work when the CCoM algorithm is utilized for better accuracy performance.
5. DISCUSSION
Fig. 4. Simulated RMS error of CoM algorithm and CRLB with PSN under different Gaussian widths. Detected photon numbers (a) Nph = 1000, (b) Nph = 100.
A. Interpixel Cross Talk Some important issues for practical applications of CCoM are discussed in this section, such as interpixel cross talk, additive noises, and multidimension extension. The pixel sensitivity is assumed to be a rectangle function as stated previously in Eq. (6) in an ideal sense, but this is often not the case in practical applications. In fact, the sensitivity function is affected by the interaction between neighboring pixels which is called “cross-talk” phenomenon [21]. A pixel can receive signals from its neighboring pixels due to optical diffraction and the diffusion of photon-generated carriers. The ratio between the active
Jia et al.
Vol. 27, No. 9 / September 2010 / J. Opt. Soc. Am. A
2043
Fig. 6. The pixel sensitivity function with cross talk 共s = 0.048兲. Dashed line: rectangle function; solid line: pixel sensitivity function with cross talk.
equivalent to expanding the Gaussian width of the image spot from to 冑2 + 2s . Therefore, the estimator of CCoM in Eq. (16) should be changed accordingly.
Fig. 5. Simulated RMS error of CCoM algorithm and CRLB with PSN under different Gaussian widths. Detected photon numbers (a) Nph = 1000, (b) Nph = 100.
pixel output and one of its neighbor outputs is reported to be about 1%–2% [22]. A simple model of cross talk is adopted to discuss the possible influence on the performance of the CCoM algorithm. We assume the spread function of detected photons to be a Gaussian function with a Gaussian width s. Thus the sensitivity of pixel can be expressed as the convolution of the rectangle function and the detected photon spread function, 1 P共x兲 = rect共x兲 丢
冑2s
exp共− x2/22s 兲.
共31兲
Given the Gaussian width of s = 0.048 we obtained a cross talk of 2%. Figure 6 models this pixel sensitivity function which is consistent with the spot scan measurements reported in [22]. According to Eq. (31), Eq. (10) can be rewritten as Fe共s兲 = F兵fe共x兲其 = F兵I共x,0兲 丢 P共x兲其 = exp关− 2共s兲2兴sinc共s兲exp关− 2共ss兲2兴 = exp关− 2共s兲2共2 + 2s 兲兴sinc共s兲.
共32兲
Equation (32) indicates that the influence of cross talk is
B. Additive Noise As mentioned previously, besides PSN, there is other noise, including background noise, dark current noise, and readout noise. Taking those noise sources into consideration we have reformulated and simulated results to illustrate the influence of noise upon the centroid estimation. As having no relationship with the photon signal, those noises can be considered as additive noises and the noise in each pixel is represented by a Gaussian random number n with variance 2n and zero mean. Then, from Eq. (30) the signal detected by each pixel can be simulated as Ui = R共Nphgi共x0兲兲 + n共n兲.
共33兲
The RMS centroid error of the CoM algorithm is calculated from Eqs. (5) and (33), and the RMS centroid error of the CCoM algorithm is calculated similarly from Eqs. (5), (16), and (33). In practical situations, with deep cooling of the commercially available CCD with a pixel width of 20 m, the noise n is approximately 3e [17]. Thus we give the simulation with 1000 incident photons under the noise levels n = 10e and n = 5e to evaluate the performance of CCoM. The simulated results are given in Fig. 7. The figure indicates that the additive noise definitely affects the performance of both CCoM and CoM, and the influence is proportional to the sampling window. This is because the total noise increases as the sampling window expands. When a 3 pixel sampling window is utilized, the accuracy performance of CCoM is better than CoM in both noise levels, and the lowest RMS error with n = 10e is 1/50 pixel. Therefore, a small sampling window is recommended when the additive noise is considered. C. Multidimension Extension In practical situations, the proposed CCoM algorithm needs to be applied in two or three dimensions. We ex-
2044
J. Opt. Soc. Am. A / Vol. 27, No. 9 / September 2010
Jia et al.
Fig. 7. Simulated RMS error of CCoM and CoM with PSN and additive noise utilizing different sampling windows. Line with dots: results of CoM; line with open circles: results of CCoM. Horizontal axis: Gaussian width of spot image; vertical axis: RMS error. (a) 3 pixels sampling window, n = 5e; (b) 3 pixels sampling window, n = 10e; (c) 5 pixels sampling window, n = 5e; (d) 5 pixels sampling window, n = 10e.
g共x,y兲 = 共I共x,x0兲I共y,y0兲 丢 rect共x兲rect共y兲兲comb共x,y兲
pand the algorithm to two dimensions in this subsection and three-dimensional situation can be easily developed. From Eqs. (1) and (8), the 2-D systematic error of the CoM algorithm can be given as
␦x =
冕冕 冕冕 ⬁
⬁
−⬁
−⬁
⬁
⬁
−⬁
= fx共x兲fy共y兲comb共x兲comb共y兲 = gx共x兲gy共y兲. Therefore
xg共x,y兲dxdy − x0 = − g共x,y兲dxdy
G共0,0兲/ u 2iG共0,0兲
G共u,v兲 = F兵g共x,y兲其 = Gu共u兲Gv共v兲, then
−⬁
冕冕 冕冕 ⬁
⬁
−⬁
−⬁
⬁
⬁
−⬁
共38兲
− x0 ,
共34兲
␦y =
共37兲
u
− y0 = −
G共0,0兲/ v 2iG共0,0兲
− y0 ,
␦x = −
−⬁
共35兲 where G共u , v兲 is the 2-D Fourier transform of the sampled pixel signal function g共x , y兲, which can be expressed as g共x,y兲 = f共x,y兲comb共x,y兲 = 共I共x,y,x0,y0兲 丢
u
Gu共u兲.
共39兲
Substituting Eqs. (38) and (39) into Eq. (34) yields
yg共x,y兲dxdy
g共x,y兲dxdy
G共u,v兲 = Gv共v兲
rect共x,y兲兲comb共x,y兲.
共36兲
Gu共0兲/ u 2iGu共0兲
共40兲
− x0 .
Therefore, we can get the same CCoM centroid estimator xˆc as Eq. (16), and for the y-axis,
yˆc −
1
兺yU i
兺
⬁ n=1
共− 1兲 exp关− 2共n兲 兴sin共2yˆcn兲/n = n
2
i
兺
i
. Ui
i
Under the Gaussian spot assumption, we can get
共41兲
Jia et al.
Vol. 27, No. 9 / September 2010 / J. Opt. Soc. Am. A
6. CONCLUSION
tion algorithm limited by photon noise for point object,” Opt. Commun. 282, 1526–1530 (2009). S. Thomas, T. Fusco, A. Tokovinin, M. Nicolle, V. Michau, and G. Rousset, “Comparison of centroid computation algorithms in a Shack–Hartmann sensor,” Mon. Not. R. Astron. Soc. 371, 323–336 (2006). A. J. Berglund, M. D. McMahon, J. J. McClelland, and J. A. Liddle, “Fast, bias-free algorithm for tracking single particles with variable size and shape,” Opt. Express 16, 14064–14075 (2008). H. Nobach, N. Damaschke, and C. Tropea, “High-precision sub-pixel interpolation in particle image velocimetry image processing,” Exp. Fluids 39, 299–304 (2005). R. B. Fisher and D. K. Naidu, “A comparison of algorithms for subpixel peak detection,” in Image Technology, Advances in Image Processing, Multimedia and Machine Vision, J. L. C. Sanz, ed. (Springer-Verlag, 1996), pp. 385–404. C. C. Liebe, “Star trackers for attitude determination,” IEEE Aerosp. Electron. Syst. Mag. 10, 10–16 (1995). K. C. N. Brian and F. Alexander, “Elimination of systematic error in subpixel accuracy centroid estimation,” Opt. Eng. (Bellingham) 30, 1320–1330 (1991). B. R. Hancock, R. C. Stirbl, T. J. Cunningham, B. Pain, C. J. Wrigley, and P. G. Ringold, “CMOS active pixel sensor specific performance effects on star tracker/imager position accuracy,” Proc. SPIE 4284, 43–53 (2001). S. B. Grossman and R. B. Emmons, “Performance analysis and size optimization of focal planes for point-source tracking algorithm applications,” Opt. Eng. (Bellingham) 23, 167–176 (1984). H. Kawano, H. Shimoji, S. Yoshikawa, K. Miyatake, K. Hama, and S. Nakamura, “Suppression of sun interference in the star sensor baffling stray light by total internal reflection,” Proc. SPIE 5962, 59621R (2005). B. R. Van, “SIRTF autonomous star tracker,” Proc. SPIE 4850, 108–121 (2002). A. O’Grady, “A comparison of EMCCD, CCD and emerging technologies optimized for low light spectroscopy applications,” Proc. SPIE 6093, 60930S (2006). B. J. Rye and R. M. Hardesty, “Discrete spectral peak estimation in incoherent backscatter heterodyne lidar. I. Spectral accumulation and the Cramer–Rao lower bound,” IEEE Trans. Geosci. Remote Sens. 31, 16–27 (1993). S. Johnson and S. Cain, “Bound on range precision for shotnoise limited ladar systems,” Appl. Opt. 47, 5147–5154 (2008). S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory (Prentice Hall, 1993). K. Hirakawa, “Cross-talk explained,” in Proceedings of IEEE Conference on Image Processing (ICIP 2008) (IEEE, 2008), pp. 677–680. M. Estribeau and P. Magnan, “Pixel crosstalk and correlation with modulation transfer function of CMOS image sensor,” Proc. SPIE 5677, 98–108 (2005).
We propose an unbiased estimation algorithm of subpixel centroid estimation which compensates the systematic error of the center of mass (CoM) algorithm. The CRLB under the condition of photon shot noise (PSN) is derived, and through the numerical simulation it clearly shows that the proposed algorithm is a kind of minimum variance estimator of centroid position estimation. Additional investigation also highlighted that the best Gaussian width of the signal spot is about 0.2–0.3 pixel using this estimation algorithm. The performance in terms of the accuracy had also been improved remarkably compared with the CoM method. Simulation result shows that the algorithm accuracy can achieve as close as 1/100 pixel even only 1000 photons are detected. This is meaningful especially in weak illumination condition such as star tracker where the centroids of dim stars must be estimated. The drawback of this method is that the Gaussian width of the image system should be pre-known in order to compensate the systematic error. Although in usual conditions the Gaussian width can be determined through laboratory calibration, the Gaussian width may change as the physical condition changes, such as the temperature variance. In this way, the robustness of the algorithm proposed in this paper is to be further studied and improved in our future endeavors.
REFERENCES 1. 2. 3. 4.
5.
6.
K. A. Winick, “Cramér–Rao lower bounds on the performance of charge-coupled-device optical position estimators,” J. Opt. Soc. Am. A 3, 1809–1815 (1986). E. Shalom, J. W. Alexander, and R. H. Stanton, “Acquisition and track algorithms for the ASTROS star tracker,” Adv. Astronaut. Sci. 57, 375–398 (1985). C. C. Liebe, “Accuracy performance of star trackers—A tutorial,” IEEE Trans. Aerosp. Electron. Syst. 38, 587–599 (2002). B. M. Quine, V. Tarasyuk, H. Mebrahtu, and R. Hornsey, “Determining star-image location: A new sub-pixel interpolation technique to process image centroids,” Comput. Phys. Commun. 177, 700–706 (2007). H. Li, H. Song, C. Rao, and X. Rao, “Accuracy analysis of centroid calculated by a modified center detection algorithm for Shack-Hartmann wavefront sensor,” Opt. Commun. 281, 750–755 (2008). H. Chen and C. Rao, “Accuracy analysis on centroid estima-
7.
8.
9. 10.
11. 12. 13.
14.
15.
16. 17. 18.
19. 20. 21. 22.
2045
