# Automatic import for live code
import numpy as np
from shrbk.plot import *
from shrbk.data import *
9.1. Moments Basics¶
In this section we will take a look at the basics of using moments to measure shapes.
Note
Try out the interactive version of this notebook on Binder!
# Import dependencies
import numpy as np
from shrbk.plot import *
from shrbk.data import *
from IPython.display import Markdown as md
from shrbk.interact import get_url, make_html_binder_button
# Provide binder badge
url = get_url('moments.ipynb')
md(make_html_binder_button(url))
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
/tmp/ipykernel_4397/2340935711.py in <module>
1 # Import dependencies
2 import numpy as np
----> 3 from shrbk.plot import *
4 from shrbk.data import *
5 from IPython.display import Markdown as md
ModuleNotFoundError: No module named 'shrbk'
9.1.1. Example Data¶
Let’s begin by looking at some dummy images with properties we can easily discern by eye. We will use the data_dict imported from data.py (see appendix) and print the available keys.
for key in data_dict.keys():
print(key)
Now we can examime the content of the matrices that define the images.
data_name = 'Diagonal Line (Down)'
print(data_dict[data_name])
Hint
Try changing data_name to one of the other possible dictionary keys to compare the matrices.
What you should notice is that we have a collection of (very unrealistic) images with a range of orientions and ellipticities. We can use these images to verify in the following sections that we fully understand how to measure the shape of an object using moments. For the purposes of this section we will treat these example images as postage stamps of isolated galaxy images.
9.1.2. Surface Brightness¶
We can define the surface brightness, \(I\), of an object as the amount of light (i.e. flux) per square arcsecond on the sky. \(I(x, y)\) gives the surface brightness at position \(x, y\), which we can relate to a given pixel. To measure the flux density, \(S_\nu\), of an object we need to intragte acoss the whole image
or, for a discrete set of pixels, sum the values of \(I(x, y)\).
Note
Surface brightness is also referred to as specific intensity, spectral intensity, spectral brightness, spectral radiance, or just brightness (see here for more details).
Warning
Here we assume that the object is isolated so that the surface brightness can be measured across the whole image (see e.g. [BartelmannSchneider01]). We will revist this point in the next section.
print('F = {}'.format(np.sum(data_dict['Horizontal Line'])))
9.1.3. Centroid¶
Now that we have established what the surface brightness is we can use it identfy the first moment or centroid \((\bar{x}, \bar{y})\) of a given image.
The following cell provides a simple implementation of equations (9.3) and (9.4).
def get_centroid(data):
# Sum flux over x and y individually
sum_i = np.array([np.sum(data, axis=i) for i in (1, 0)])
# Get range of x and y values
ranges = np.array([np.arange(i) for i in data.shape])
# Calculate centroids
cents = np.sum(sum_i * ranges, axis=1) / np.sum(data)
return cents.astype(int)
Now let’s test that our function works on the example data.
def show_centroid(data_name):
data = data_dict[data_name]
centroid = get_centroid(data)
show_image(data, centroid=centroid)
show_centroid('Circle (Off-Centre)')
We can see that, regardless of the position of the image within the postage stamp, the centroid (highlighted with a red x) is always correctly determined.
9.1.4. Moments¶
Finally, we can define the second moment or quadrupole moments of an image as follows.
where \(k_1,k_2 \in \{x,y\}\), the quadrupoles are \(Q_{xx},Q_{xy}\)and \(Q_{yy}\) (with \(Q_{xy}=Q_{yx}\)). Equation (9.5) is used to determine object shapes. In the following cell we a simple implementation of this equation.
def get_moments(data):
centroid = get_centroid(data)
ranges = np.array([np.arange(i) for i in data.shape])
x = np.outer(ranges[0] - centroid[0], np.ones(data.shape[1]))
y = np.outer(np.ones(data.shape[0]), ranges[1] - centroid[1])
q = np.array([np.sum(data * xi * xj) for xi in (x, y) for xj in (x, y)])
q = (q / np.sum(data)).reshape(2, 2).astype('complex')
return q
We can examine the output of this function for the various example images.
print(get_moments(data_dict['Vertical Line']))
While it is clear that the moments are different for each image, it can be difficult to interpret these numbers.
9.1.5. Measures of Ellipticity¶
It is convenient to convert these moments into a measure of the object’s ellipticity, however this can be done in various ways.
9.1.5.1. Polarisation¶
This is the standard convetion used by methods such as KSB and is often denoted by the greek letter \(\chi\).
9.1.5.2. Ellipticity¶
This is an alternative convention that… and is often denoted by the greek letter \(\epsilon\).
Note
We have adopted the convention of labeling \(Q_{x, y}\) starting from \(x, y\) = 0 to keep the equations consistent with the code.
The following cell provides a simple implementation of both (9.6) and (9.7).
def get_ellipticity(data, method='chi'):
# Calculate moments
q = get_moments(data)
# Calculate the image size.
r2 = q[0, 0] + q[1, 1]
# Calculate the numerator
num = (q[0, 0] - q[1, 1] + 2 * np.complex(0, q[0, 1]))
# Calculate the denominator
den = r2
if method == 'epsilon':
den += 2 * np.sqrt(q[0, 0] * q[1, 1] - q[0, 1] ** 2)
# Calculate the ellipticity/polarisation
ellip = num / den
return np.around([ellip.real, ellip.imag], 3)
The ellipticity makes it a bit easier to interpret the properties of the example images.
def show_ellipticity(data_name):
data = data_dict[data_name]
centroid = get_centroid(data)
chi = get_ellipticity(data)
epsilon = get_ellipticity(data, method='epsilon')
show_image(data, centroid=centroid, chi=chi, epsilon=epsilon)
show_ellipticity('Diagonal Line (Up)')
9.1.5.3. Semi-Major/Minor Axes¶
It is also possible to relate the quadrupole moments to more traditional measures of ellipticity, namely the semi-major axis (\(a\)), semi-minor axis (\(b\)) and angle (\(\theta\)).
The following cell implements equations (9.8), (9.9) and (9.10).
def get_abt(data):
q = get_moments(data)
qq_plus = q[0, 0] + q[1, 1]
qq_minus = q[0, 0] - q[1, 1]
root = np.sqrt(qq_minus ** 2 + 4 * q[0, 1] ** 2)
a = np.around(np.real(np.sqrt(0.5 * (qq_plus + root))), 3)
b = np.around(np.real(np.sqrt(0.5 * (qq_plus - root))), 3)
if qq_minus == 0.0:
theta = -45.0 * np.sign(np.real(q[0, 1]))
else:
theta = np.around(np.real(0.5 * np.arctan(2 * q[0, 1] / qq_minus)) * 180. / np.pi, 3)
if qq_minus > 0.0 and q[0, 1] == 0.0:
theta += 90.0
elif qq_minus > 0.0:
theta -= 90.0 * np.sign(np.real(q[0, 1]))
return a, b, theta
We can test the accuracy of this implementation by generating some images where we define the properties of the ellipse and seeing if we can recover the correct model from the moments.
def show_ellipse(a, b, theta):
data = make_ellipse(a, b, theta)
centroid = get_centroid(data)
chi = get_ellipticity(data)
epsilon = get_ellipticity(data, method='epsilon')
abt = get_abt(data)
show_image(data, centroid=centroid, chi=chi, epsilon=epsilon, abt=abt)
show_ellipse(a=5, b=1, theta=15)
Note
Refer to the appendix to see how the input images are produced.