GNU Astronomy Utilities

7.4.4 Measuring elliptical parameters

The shape or morphology of a target is one of the most commonly desired parameters of a target. Here, we will review the derivation of the most basic/simple morphological parameters: the elliptical parameters for a set of labeled pixels. The elliptical parameters are: the (semi-)major axis, the (semi-)minor axis and the position angle along with the central position of the profile. The derivations below follow the SExtractor manual derivations with some added explanations for easier reading.

Let’s begin with one dimension for simplicity: Assume we have a set of \(N\) values \(B_i\) (for example, showing the spatial distribution of a target’s brightness), each at position \(x_i\). The simplest parameter we can define is the geometric center of the object (\(x_g\)) (ignoring the brightness values): \(x_g=(\sum_ix_i)/N\). Moments are defined to incorporate both the value (brightness) and position of the data. The first moment can be written as:

$$\overline{x}={\sum_iB_ix_i \over \sum_iB_i}$$

This is essentially the weighted (by \(B_i\)) mean position. The geometric center (\(x_g\), defined above) is a special case of this with all \(B_i=1\). The second moment is essentially the variance of the distribution:

$$\overline{x^2}\equiv{\sum_iB_i(x_i-\overline{x})^2 \over \sum_iB_i} = {\sum_iB_ix_i^2 \over \sum_iB_i} - 2\overline{x}{\sum_iB_ix_i\over\sum_iB_i} + \overline{x}^2 ={\sum_iB_ix_i^2 \over \sum_iB_i} - \overline{x}^2$$

The last step was done from the definition of \(\overline{x}\). Hence, the square root of \(\overline{x^2}\) is the spatial standard deviation (along the one-dimension) of this particular brightness distribution (\(B_i\)). Crudely (or qualitatively), you can think of its square root as the distance (from \(\overline{x}\)) which contains a specific amount of the flux (depending on the \(B_i\) distribution). Similar to the first moment, the geometric second moment can be found by setting all \(B_i=1\). So while the first moment quantified the position of the brightness distribution, the second moment quantifies how that brightness is dispersed about the first moment. In other words, it quantifies how “sharp” the object’s image is.

Before continuing to two dimensions and the derivation of the elliptical parameters, let’s pause for an important implementation technicality. You can ignore this paragraph and the next two if you do not want to implement these concepts. The basic definition (first definition of \(\overline{x^2}\) above) can be used without any major problem. However, using this fraction requires two runs over the data: one run to find \(\overline{x}\) and another run to find \(\overline{x^2}\) from \(\overline{x}\), this can be slow. The advantage of the last fraction above, is that we can estimate both the first and second moments in one run (since the \(-\overline{x}^2\) term can easily be added later).

The logarithmic nature of floating point number digitization creates a complication however: suppose the object is located between pixels 10000 and 10020. Hence the target’s pixels are only distributed over 20 pixels (with a standard deviation \(<20\)), while the mean has a value of \(\sim10000\). The \(\sum_iB_i^2x_i^2\) will go to very very large values while the individual pixel differences will be orders of magnitude smaller. This will lower the accuracy of our calculation due to the limited accuracy of floating point operations. The variance only depends on the distance of each point from the mean, so we can shift all position by a constant/arbitrary \(K\) which is much closer to the mean: \(\overline{x-K}=\overline{x}-K\). Hence we can calculate the second order moment using:

$$\overline{x^2}={\sum_iB_i(x_i-K)^2 \over \sum_iB_i} - (\overline{x}-K)^2 $$

The closer \(K\) is to \(\overline{x}\), the better (the sums of squares will involve smaller numbers), as long as \(K\) is within the object limits (in the example above: \(10000\leq{K}\leq10020\)), the floating point error induced in our calculation will be negligible. For the most simplest implementation, MakeCatalog takes \(K\) to be the smallest position of the object in each dimension. Since \(K\) is arbitrary and an implementation/technical detail, we will ignore it for the remainder of this discussion.

In two dimensions, the mean and variances can be written as:

$$\overline{x}={\sum_iB_ix_i\over B_i}, \quad \overline{x^2}={\sum_iB_ix_i^2 \over \sum_iB_i} - \overline{x}^2$$ $$\overline{y}={\sum_iB_iy_i\over B_i}, \quad \overline{y^2}={\sum_iB_iy_i^2 \over \sum_iB_i} - \overline{y}^2$$ $$\quad\quad\quad\quad\quad\quad\quad\quad\quad \overline{xy}={\sum_iB_ix_iy_i \over \sum_iB_i} - \overline{x}\times\overline{y}$$

If an elliptical profile’s major axis exactly lies along the \(x\) axis, then \(\overline{x^2}\) will be directly proportional with the profile’s major axis, \(\overline{y^2}\) with its minor axis and \(\overline{xy}=0\). However, in reality we are not that lucky and (assuming galaxies can be parameterized as an ellipse) the major axis of galaxies can be in any direction on the image (in fact this is one of the core principles behind weak-lensing by shear estimation). So the purpose of the remainder of this section is to define a strategy to measure the position angle and axis ratio of some randomly positioned ellipses in an image, using the raw second moments that we have calculated above in our image coordinates.

Let’s assume we have rotated the galaxy by \(\theta\), the new second order moments are:

$$\overline{x_\theta^2} = \overline{x^2}\cos^2\theta + \overline{y^2}\sin^2\theta - 2\overline{xy}\cos\theta\sin\theta $$ $$\overline{y_\theta^2} = \overline{x^2}\sin^2\theta + \overline{y^2}\cos^2\theta + 2\overline{xy}\cos\theta\sin\theta$$ $$\overline{xy_\theta} = \overline{x^2}\cos\theta\sin\theta - \overline{y^2}\cos\theta\sin\theta + \overline{xy}(\cos^2\theta-\sin^2\theta)$$

The best \(\theta\) (\(\theta_0\), where major axis lies along the \(x_\theta\) axis) can be found by:

$$\left.{\partial \overline{x_\theta^2} \over \partial \theta}\right|_{\theta_0}=0$$ Taking the derivative, we get: $$2\cos\theta_0\sin\theta_0(\overline{y^2}-\overline{x^2}) + 2(\cos^2\theta_0-\sin^2\theta_0)\overline{xy}=0$$ When \(\overline{x^2}\neq\overline{y^2}\), we can write: $$\tan2\theta_0 = 2{\overline{xy} \over \overline{x^2}-\overline{y^2}}.$$

MakeCatalog uses the standard C math library’s atan2 function to estimate \(\theta_0\), which we define as the position angle of the ellipse. To recall, this is the angle of the major axis of the ellipse with the \(x\) axis. By definition, when the elliptical profile is rotated by \(\theta_0\), then \(\overline{xy_{\theta_0}}=0\), \(\overline{x_{\theta_0}^2}\) will be the extent of the maximum variance and \(\overline{y_{\theta_0}^2}\) the extent of the minimum variance (which are perpendicular for an ellipse). Replacing \(\theta_0\) in the equations above for \(\overline{x_\theta}\) and \(\overline{y_\theta}\), we can get the semi-major (\(A\)) and semi-minor (\(B\)) lengths:

$$A^2\equiv\overline{x_{\theta_0}^2}= {\overline{x^2} + \overline{y^2} \over 2} + \sqrt{\left({\overline{x^2}-\overline{y^2} \over 2}\right)^2 + \overline{xy}^2}$$

$$B^2\equiv\overline{y_{\theta_0}^2}= {\overline{x^2} + \overline{y^2} \over 2} - \sqrt{\left({\overline{x^2}-\overline{y^2} \over 2}\right)^2 + \overline{xy}^2}$$

As a summary, it is important to remember that the units of \(A\) and \(B\) are in pixels (the standard deviation of a positional distribution) and that they represent the spatial light distribution of the object in both image dimensions (rotated by \(\theta_0\)). When the object cannot be represented as an ellipse, this interpretation breaks down: \(\overline{xy_{\theta_0}}\neq0\) and \(\overline{y_{\theta_0}^2}\) will not be the direction of minimum variance.