An unbiased estimator for the ellipticity from image moments

TL;DR

This paper introduces a statistically unbiased estimator for the ellipticity of objects in noisy images, based on image moments, assuming Gaussian noise and known point-spread function, improving accuracy in shape measurement.

Contribution

The paper presents a novel unbiased estimator for ellipticity derived from image moments, accounting for noise and PSF, with an explicit covariance estimate.

Findings

Estimator is unbiased under specified assumptions

Covariance of the estimator can be explicitly computed

Method improves shape measurement accuracy in noisy images

Abstract

An unbiased estimator for the ellipticity of an object in a noisy image is given in terms of the image moments. Three assumptions are made: i) the pixel noise is normally distributed, although with arbitrary covariance matrix, ii) the image moments are taken about a fixed centre, and iii) the point-spread function is known. The relevant combinations of image moments are then jointly normal and their covariance matrix can be computed. A particular estimator for the ratio of the means of jointly normal variates is constructed and used to provide the unbiased estimator for the ellipticity. Furthermore, an unbiased estimate of the covariance of the new estimator is also given.

Tables

Table 1. Table 1: Monte Carlo results for the unbiased ellipticity estimator

${\chi }_1$	${\chi }_2$	PSF	counts	$⟨{\hat{\chi }}_1⟩$	$⟨{\hat{\chi }}_2⟩$	$\mathrm{Var}{\left[{\hat{\chi }}_1\right]}^{1/2}$	$\mathrm{Var}{\left[{\hat{\chi }}_2\right]}^{1/2}$	${⟨\hat{V}\mathrm{ar}\left[{\hat{\chi }}_1\right]⟩}^{1/2}$	${⟨\hat{V}\mathrm{ar}\left[{\hat{\chi }}_2\right]⟩}^{1/2}$
$0.1107$	$0.0000$	no	500	$0.1113\left(08\right)$	$0.0006\left(11\right)$	$0.762$	$1.112$	$0.762$	$1.111$
			1000	$0.1105\left(02\right)$	$-0.0002\left(02\right)$	$0.216$	$0.214$	$0.216$	$0.214$
			2000	$0.1108\left(01\right)$	$0.0000\left(01\right)$	$0.104$	$0.103$	$0.104$	$0.103$
		yes	500	$0.1115\left(07\right)$	$0.0006\left(06\right)$	$0.717$	$0.598$	$0.717$	$0.599$
			1000	$0.1111\left(02\right)$	$-0.0001\left(02\right)$	$0.216$	$0.214$	$0.215$	$0.214$
			2000	$0.1108\left(01\right)$	$0.0000\left(01\right)$	$0.104$	$0.103$	$0.104$	$0.103$
$0.1820$	$-0.1842$	no	500	$0.1820\left(07\right)$	$-0.1840\left(09\right)$	$0.704$	$0.853$	$0.704$	$0.853$
			1000	$0.1817\left(02\right)$	$-0.1842\left(02\right)$	$0.225$	$0.226$	$0.225$	$0.226$
			2000	$0.1822\left(01\right)$	$-0.1842\left(01\right)$	$0.108$	$0.108$	$0.108$	$0.108$
		yes	500	$0.1802\left(07\right)$	$-0.1840\left(06\right)$	$0.658$	$0.616$	$0.658$	$0.616$
			1000	$0.1819\left(02\right)$	$-0.1845\left(02\right)$	$0.223$	$0.224$	$0.224$	$0.225$
			2000	$0.1819\left(01\right)$	$-0.1841\left(01\right)$	$0.108$	$0.108$	$0.108$	$0.108$
$0.0000$	$0.5492$	no	500	$-0.0005\left(08\right)$	$0.5489\left(15\right)$	$0.793$	$1.451$	$0.794$	$1.450$
			1000	$-0.0001\left(02\right)$	$0.5494\left(03\right)$	$0.248$	$0.323$	$0.248$	$0.322$
			2000	$-0.0002\left(01\right)$	$0.5493\left(01\right)$	$0.117$	$0.149$	$0.117$	$0.149$
		yes	500	$0.0001\left(07\right)$	$0.5473\left(10\right)$	$0.720$	$0.957$	$0.720$	$0.960$
			1000	$0.0002\left(02\right)$	$0.5488\left(03\right)$	$0.244$	$0.315$	$0.244$	$0.316$
			2000	$0.0001\left(01\right)$	$0.5491\left(01\right)$	$0.116$	$0.147$	$0.116$	$0.147$

Equations

$$\chi ={\chi }_1+\mathrm{i}{\chi }_2=\frac{m_{20}-m_{02}+2\mathrm{i}{m}_{11}}{m_{20}+m_{02}}.$$

$$\chi =\frac{{\chi }_{\text{s}}+2g+g^2{\chi }_{\text{s}}^{\ast }}{1+|g{|}^2+2\Re [g{\chi }_{\text{s}}^{\ast }]},$$

$$\chi \approx {\chi }_{\text{s}}+2g.$$

$$⟨\chi ⟩=2⟨g⟩.$$

$${\xi }_{ij}(r)={⟨{\chi }_i(\mathbf{x}){\chi }_j(\mathbf{y})⟩}_{|\mathbf{x}-\mathbf{y}|=r},$$

$$\hat{e}=\frac{{\hat{m}}_{20}-{\hat{m}}_{02}+2\mathrm{i}{\hat{m}}_{11}}{{\hat{m}}_{20}+{\hat{m}}_{02}}.$$

$$\mathbf{d}\sim \mathcal{N}(\mathbf{\mu },\text{mn@boldsymbol}\Sigma ),$$

$$\chi =\frac{m_{20}-m_{02}-m_{00}({\pi }_{20}-{\pi }_{02})+2\mathrm{i}{m}_{11}-2\mathrm{i}{m}_{00}{\pi }_{11}}{m_{20}+m_{02}-m_{00}({\pi }_{20}+{\pi }_{02})},$$

$$u$$

$$v$$

$$s$$

$$(u,v,s)=\text{M}\mathbf{d},$$

$$(u,v,s)\sim \mathcal{N}({\mathbf{\mu }}_{uvs},\text{mn@boldsymbol}{\Sigma }_{uvs}),$$

$${\Sigma }_{uu}$$

$${\Sigma }_{vv}$$

$${\Sigma }_{ss}$$

$$\chi ={\chi }_1+\mathrm{i}{\chi }_2=\frac{{\mu }_u+\mathrm{i}{\mu }_v}{{\mu }_s}.$$

$$p=u-as,$$

$$q=v-bs,$$

$${\hat{\chi }}_1=a+p\hat{g}(s),{\hat{\chi }}_2=b+q\hat{g}(s).$$

$$E[{\hat{\chi }}_1]$$

$$E[{\hat{\chi }}_2]$$

$$\hat{\mathrm{V}}\mathrm{ar}[{\hat{\chi }}_1]$$

$$\hat{\mathrm{V}}\mathrm{ar}[{\hat{\chi }}_2]$$

$$\operatorname{\hat{C}ov}[\hat{\chi}_{1},\hat{\chi}_{2}]=pq\,\big{(}\hat{g}(s)\big{)}^{2}-\big{(}(pq-\Sigma_{uv})/\Sigma_{ss}+ab\big{)}\,\big{(}1-s\,\hat{g}(s)\big{)}\;.$$

$$r={\mu }_x/{\mu }_y.$$

$$w=x-cy,c=\rho {\sigma }_x/{\sigma }_y.$$

$$\hat{g}(y)=\frac{\sqrt{\pi}}{\sqrt{2}\,\sigma_{y}}\,\exp\bigg{(}\frac{y^{2}}{2\,\sigma_{y}^{2}}\bigg{)}\,\operatorname{erfc}\bigg{(}\frac{y}{\sqrt{2}\,\sigma_{y}}\bigg{)}\;.$$

$$\hat{r}=c+w\hat{g}(y).$$

$$E[\hat{r}]=c+E[w]E[\hat{g}(y)]=c+({\mu }_x-c{\mu }_y)/{\mu }_y={\mu }_x/{\mu }_y,$$

$$\operatorname{Var}[\hat{r}]=\big{(}{\operatorname{E}[w]}^{2}+\operatorname{Var}[w]\big{)}\operatorname{Var}[\hat{g}(y)]+\operatorname{Var}[w]\,\operatorname{E}[\hat{g}(y)]^{2}\;.$$

$$\Pr\big{(}|\hat{g}(y)-1/\mu_{y}|<h\big{)}=p\;,$$

$$\operatorname{\hat{V}ar}[\hat{g}(y)]=\big{(}\hat{g}(y)\big{)}^{2}-\big{(}1-y\,\hat{g}(y)\big{)}/\sigma_{y}^{2}$$

$$\operatorname{\hat{V}ar}[\hat{r}]=w^{2}\,\big{(}\hat{g}(y)\big{)}^{2}-\big{(}(w^{2}-\sigma_{x}^{2})/\sigma_{y}^{2}+c^{2}\big{)}\,\big{(}1-y\,\hat{g}(y)\big{)}$$

Full text

††pubyear: 2017††pagerange: An unbiased estimator for the ellipticity from image moments–A

1 Introduction

A number of applications in astronomy require the measurement of shapes of objects from observed images. A commonly used shape descriptor is the ellipticity $\chi$, which is defined in terms of the central image moments $m_{pq}$ as the complex number

[TABLE]

For example, in weak gravitational lensing, the gravitational field distorts the observed shapes of background galaxies, and this shear can be detected in ellipticity measurements. This is possible because the observed ellipticity of a source affected by gravitational lensing is (see e.g. Bartelmann & Schneider, 2001)

[TABLE]

where $\chi_{\text{s}}$ is the intrinsic ellipticity of the source, and $g=g_{1}+\mathop{}\!\mathrm{i}\!\mathop{}g_{2}$ is the so-called reduced shear of the gravitational lens. In the limit of weak lensing, this can be approximated to linear order as

[TABLE]

Averaging over an ensemble of randomly-oriented background sources, i.e. $\langle\chi_{\text{s}}\rangle=0$, the weak lensing equation (3) yields

[TABLE]

The observed ellipticities are thus a direct estimator for the shear $g$ from gravitational lensing. Similar ideas are used in Cosmology (for a recent review, see Kilbinger, 2015). Here, cosmic shear from the large-scale structure of the universe imprints a specific signature onto the ellipticity two-point correlation functions,

[TABLE]

where the average is taken over pairs of sources with the given separation $r$ on the sky. Note that both the one-point function (4) and the two-point function (5) depend on the mean ellipticity over a potentially large sample of sources.

In practice, an estimator $\hat{\chi}$ is used to measure the ellipticities of observed sources. In order to not introduce systematic errors into applications such as the above, the ellipticity estimator $\hat{\chi}$ must be unbiased, i.e. $\operatorname{E}[\hat{\chi}]=\chi$. One of the biggest problems for estimators that directly work on the data is the noise bias (Refregier et al., 2012) arising from pixel noise in the observations. For example, the standard approach to moment-based shape measurement is to obtain estimates $\hat{m}_{pq}$ of the second-order moments from the (noisy) data and use (1) directly as an ellipticity estimate,

[TABLE]

The statistical properties of this estimator have been studied by Melchior & Viola (2012) and Viola et al. (2014), who assumed that the image moments are jointly normal with some given variance and correlation. The estimator (6) then follows the distribution of Marsaglia (1965, 2006) for the ratio of jointly normal variates. None of the moments of this distribution exist, and even for a finite sample, small values in the denominator can quickly cause significant biases and large variances. The estimator (6) is thus generally poorly behaved, unless the signal-to-noise ratio of the data is very high.

Here, a new unbiased estimator for the ellipticity $\chi$ from the second-order image moments is proposed. First, it is shown that for normally-distributed pixel noise with known covariance matrix and a fixed centre, the relevant combinations of image moments are indeed jointly normal and that their covariance matrix can easily be computed. In the appendix, an unbiased estimator for the ratio of the means of jointly normal variates is constructed, which can subsequently be applied to the image moments. This produces the ellipticity estimate, as well as unbiased estimates of its variance and covariance.

2 An unbiased estimator for the ellipticity

It is assumed that the data is a random vector $\bm{d}=(d_{1},d_{2},\dots)$ of pixels following a multivariate normal distribution centred on the unknown true signal $\bm{\mu}$,

[TABLE]

where $\mn@boldsymbol{\Sigma}$ is the covariance matrix for the noise, which is assumed to be known but not restricted to a particular diagonal shape. The observed signal usually involves a point-spread function (PSF), and it is further assumed that this effect can be approximated as a linear convolution of the true signal and the discretised PSF. In this case, definition (1) can be extended to obtain the true ellipticity of the object before convolution,

[TABLE]

from the central moments $\pi_{pq}$ of the (normalised) PSF. Fixing a centre $(\bar{x},\bar{y})$, the relevant combinations of moments to compute the ellipticity from data $\bm{d}$ are thus

[TABLE]

where $w_{i}$ is the window function of the observation. To obtain the true ellipticity estimate of the signal, and for the PSF correction (8) to remain valid, the window function must be unity over the support of the signal.111This is in contrast to the weight functions that are sometimes used in moment-based methods to reduce the influence of noise far from the centre.

Due to the linearity in the pixel values $d_{i}$, the vector $(u,v,s)$ can be written in matrix form,

[TABLE]

where the three rows $\alpha_{i}$, $\beta_{i}$, $\gamma_{i}$ of matrix M are defined by (9). The random vector $(u,v,s)$ is hence normally distributed,

[TABLE]

with unknown mean $\bm{\mu}_{uvs}=(\mu_{u},\mu_{v},\mu_{s})=\textbf{{M}}\,\bm{\mu}$ and known $3\times 3$ covariance matrix $\mn@boldsymbol{\Sigma}_{uvs}=\textbf{{M}}\,\mn@boldsymbol{\Sigma}\,\textbf{{M}}^{\mathrm{T}}$ with entries

[TABLE]

where $\Sigma_{ij}$ are the entries of the covariance matrix $\mn@boldsymbol{\Sigma}$ of the pixel noise. Hence the covariance matrix $\mn@boldsymbol{\Sigma}_{uvs}$ of the moments can be computed if the pixel noise statistics are known.

The true ellipticity (8) of the signal can be written in terms of the mean values of the variates $u$, $v$ and $s$ defined in (9) as

[TABLE]

The problem is to find an unbiased estimate of $\chi_{1}$ and $\chi_{2}$ from the observed values of $u$, $v$ and $s$. In appendix A, an unbiased estimator is given for the ratio of the means of two jointly normal random variables. It can be applied directly to the ellipticity (13). First, two new variates $p$ and $q$ are introduced,

[TABLE]

where $a$, $b$ are constants. This definition corresponds to (20) in the univariate case, and therefore both $(p,s)$ and $(q,s)$ are pairs of independent normal variates. The desired estimator for the ellipticity is then

[TABLE]

Because the mean $\mu_{s}$ is always positive for a realistic signal, the function $\hat{g}(s)$ is given by (21). From the expectation

[TABLE]

it follows that $\hat{\chi}=\hat{\chi}_{1}+\mathop{}\!\mathrm{i}\!\mathop{}\hat{\chi}_{2}$ is indeed an unbiased estimator for the true ellipticity $\chi$ of the signal.

In addition, an unbiased estimate of the the variance of $\hat{\chi}_{1}$ and $\hat{\chi}_{2}$ is provided by (27),

[TABLE]

Similarly, there is an unbiased estimator for the covariance,

[TABLE]

It follows that the individual estimates of the ellipticity components are in general not independent. However, for realistic pixel noise, window functions and PSFs, the correlations between $u$, $v$ and $s$, and hence $\hat{\chi}_{1}$ and $\hat{\chi}_{2}$, can become very small.

To demonstrate that the proposed estimator is in fact unbiased under the given assumptions of i) normal pixel noise with known covariance, ii) a fixed centre for the image moments, and iii) the discrete convolution with a known PSF, a Monte Carlo simulation was performed using mock observations of an astronomical object. The images are postage stamps of $49\times 49$ pixels, containing a centred source that is truncated near the image borders. A circular aperture is used as window function. The source is elliptical and follows the light profile of de Vaucouleurs (1948). The ellipse containing half the total light has a 10 pixel semi-major axis and ellipticity as specified. Where indicated, the signal is convolved with a Gaussian PSF with 5 pixel FWHM. The pixel noise is uncorrelated with unit variance. The normalisation $N$ of the object (i.e. the total number of counts) varies to show the effect of the signal-to-noise ratio on the variance of the estimator.222To compare the results to a given data set, it is then possible to scale the data so that the noise has unit variance, and compare the number of counts. For example, the control-ground-constant data set of the GREAT3 challenge (Mandelbaum et al., 2014) has mostly $N=500$–$1000$.

The signal ellipticity $\chi$ is computed from the image before the PSF is applied and noise is added. The mean of the estimator $\hat{\chi}$ is computed from $10^{6}$ realisations of noise for the same signal. Also computed are the square root of the sample variance $\operatorname{Var}[\hat{\chi}]$ and the mean of the estimated variance $\operatorname{\hat{V}ar}[\hat{\chi}]$, respectively. The results shown in Table 1 indicate that the estimator performs as expected.

3 Conclusion & discussion

The unbiased estimator (15) provides a new method for ellipticity measurement from noisy images. Its simple and analytic form allows quick implementation and fast evaluation, and statistics for the results can be obtained directly with unbiased estimates of the variance (17) and covariance (18).

Using an unbiased estimator for the ellipticity of the signal $\bm{\mu}$ eliminates the influence of noise from a subsequent analysis of the results (the so-called “noise bias” in weak lensing, Refregier et al., 2012). However, depending on the application, other kinds of biases might exist even for a noise-free image. For example, due to the discretisation of the image, the signal ellipticity can differ from the intrinsic ellipticity of the observed object. This “pixellation bias” (Simon & Schneider, 2016) remains an issue in applications such as weak lensing, where the relevant effects must be measured from the intrinsic ellipticity of the objects.

Furthermore, a fixed centre $(\bar{x},\bar{y})$ for the moments has been assumed throughout. For a correct ellipticity estimate, this must be the centroid of the signal, which is usually estimated from the data itself. Centroid errors (Melchior & Viola, 2012) might therefore ultimately bias the ellipticity estimator or increase its variance, although this currently does not seem to be a significant effect.

In practice, additional biases might arise when the assumed requirements for the window function and PSF are not fulfilled by the data. The estimator should therefore always be carefully tested for the application at hand.

Lastly, the ellipticity estimate might be improved by suitable filtering of the observed image. A linear filter with matrix A can be applied to the image before estimating the ellipticity, since the transformed pixels remain multivariate normal with mean $\bm{\mu}^{\prime}=\textbf{{A}}\,\bm{\mu}$ and covariance matrix $\mn@boldsymbol{\Sigma}^{\prime}=\textbf{{A}}\,\mn@boldsymbol{\Sigma}\,\textbf{{A}}^{\mathrm{T}}$. Examples of viable filters are nearest-neighbour or bilinear interpolation, as well as convolution. A combination of these filters could be used to perform PSF deconvolution on the observed image, as an alternative to the algebraic PSF correction (8).

Acknowledgements

NT would like to thank S. Bridle for encouragement and many conversations about shape measurement. The author acknowledges support from the European Research Council in the form of a Consolidator Grant with number 681431.

Appendix A A ratio estimator for normal variates

Let $x$ and $y$ be jointly normal variates with unknown means $\mu_{x}$ and $\mu_{y}$, and known variances $\sigma_{x}^{2}$ and $\sigma_{y}^{2}$ and correlation $\rho$. The goal here is to find an unbiased estimate of the ratio $r$ of their means,

[TABLE]

Under the additional assumption that the sign of $\mu_{y}$ is known, an unbiased estimator for $r$ can be found in two short steps.

First, the transformation of Marsaglia (1965, 2006) is used to construct a new variate,

[TABLE]

The constant $c$ is chosen so that $\operatorname{E}[wy]=\operatorname{E}[w]\operatorname{E}[y]$. It is clear that $w$ is normal, and that variates $w$ and $y$ are jointly normal, uncorrelated, and thus independent. Note that $c=0$ and $w=x$ for independent $x$ and $y$.

Secondly, Voinov (1985) derived an unbiased estimator for the inverse mean of the normal variate $y$, i.e. a function $\hat{g}(y)$ with $\operatorname{E}[\hat{g}(y)]=1/\mu_{y}$. For the relevant case of an unknown but positive mean $\mu_{y}>0$, this function is given by

[TABLE]

It is then straightforward to construct an estimator for the ratio (19),

[TABLE]

Since $w$ and $y$ are independent, the expectation is

[TABLE]

which shows that $\hat{r}$ is in fact an unbiased estimator for $r$.

The variance of the ratio estimator $\hat{r}$ is formally given by

[TABLE]

As pointed out by Voinov (1985), the variance $\operatorname{Var}[\hat{g}(y)]$ does not exist for function (21) due to a divergence at infinity. The confidence interval $h$ with probability $p$,

[TABLE]

however, is well-defined, and the variance of $\hat{g}(y)$ remains finite in applications where infinite values of $y$ are not observed. In this case, an unbiased estimator

[TABLE]

exists and, together with Voinov’s estimator for $1/\mu_{y}^{2}$, yields

[TABLE]

as an unbiased estimate of the variance of the estimator $\hat{r}$.

When $y$ is significantly larger than its standard deviation, e.g. $y/\sigma_{y}>10$, the function $\hat{g}(y)$ given by (21) is susceptible to numerical overflow. However, in this regime, it is also very well approximated by its series expansion,

[TABLE]

For even larger values $y\gg\sigma_{y}$, this approaches $1/y$, as expected.