GNU Astronomy Utilities Defining an ellipse and ellipsoid

The PSF, see Point spread function, and galaxy radial profiles are generally defined on an ellipse. Therefore, in this section we will start defining an ellipse on a pixelated 2D surface. Labeling the major axis of an ellipse \(a\), and its minor axis with \(b\), the axis ratio is defined as: \(q\equiv b/a\). The major axis of an ellipse can be aligned in any direction, therefore the angle of the major axis with respect to the horizontal axis of the image is defined to be the position angle of the ellipse and in this book, we show it with \(\theta\).

Our aim is to put a radial profile of any functional form \(f(r)\) over an ellipse. Hence we need to associate a radius/distance to every point in space. Let’s define the radial distance \(r_{el}\) as the distance on the major axis to the center of an ellipse which is located at \(i_c\) and \(j_c\) (in other words \(r_{el}\equiv{a}\)). We want to find \(r_{el}\) of a point located at \((i,j)\) (in the image coordinate system) from the center of the ellipse with axis ratio \(q\) and position angle \(\theta\). First the coordinate system is rotated226 by \(\theta\) to get the new rotated coordinates of that point \((i_r,j_r)\):

$$i_r(i,j)=+(i_c-i)\cos\theta+(j_c-j)\sin\theta$$ $$j_r(i,j)=-(i_c-i)\sin\theta+(j_c-j)\cos\theta$$

Recall that an ellipse is defined by \((i_r/a)^2+(j_r/b)^2=1\) and that we defined \(r_{el}\equiv{a}\). Hence, multiplying all elements of the ellipse definition with \(r_{el}^2\) we get the elliptical distance at this point point located: \(r_{el}=\sqrt{i_r^2+(j_r/q)^2}\). To place the radial profiles explained below over an ellipse, \(f(r_{el})\) is calculated based on the functional radial profile desired.

An ellipse in 3D, or an ellipsoid, can be defined following similar principles as before. Labeling the major (largest) axis length as \(a\), the second and third (in a right-handed coordinate system) axis lengths can be labeled as \(b\) and \(c\). Hence we have two axis ratios: \(q_1\equiv{b/a}\) and \(q_2\equiv{c/a}\). The orientation of the ellipsoid can be defined from the orientation of its major axis. There are many ways to define 3D orientation and order matters. So to be clear, here we use the ZXZ (or \(Z_1X_2Z_3\)) proper Euler angles to define the 3D orientation. In short, when a point is rotated in this order, we first rotate it around the Z axis (third axis) by \(\alpha\), then about the (rotated) X axis by \(\beta\) and finally about the (rotated) Z axis by \(\gamma\).

Following the discussion in Merging multiple warpings, we can define the full rotation with the following matrix multiplication. However, here we are rotating the coordinates, not the point. Therefore, both the rotation angles and rotation order are reversed. We are also not using homogeneous coordinates (see Linear warping basics) since we are not concerned with translation in this context:

$$\left[\matrix{i_r\cr j_r\cr k_r}\right] = \left[\matrix{cos\gamma&sin\gamma&0\cr -sin\gamma&cos\gamma&0\cr 0&0&1}\right] \left[\matrix{1&0&0\cr 0&cos\beta&sin\beta\cr 0&-sin\beta&cos\beta }\right] \left[\matrix{cos\alpha&sin\alpha&0\cr -sin\alpha&cos\alpha&0\cr 0&0&1}\right] \left[\matrix{i_c-i\cr j_c-j\cr k_c-k}\right] $$

Recall that an ellipsoid can be characterized with \((i_r/a)^2+(j_r/b)^2+(k_r/c)^2=1\), so similar to before (\(r_{el}\equiv{a}\)), we can find the ellipsoidal radius at pixel \((i,j,k)\) as: \(r_{el}=\sqrt{i_r^2+(j_r/q_1)^2+(k_r/q_2)^2}\).

MakeProfiles builds the profile starting from the nearest element (pixel in an image) in the dataset to the profile center. The profile value is calculated for that central pixel using Monte Carlo integration, see Sampling from a function. The next pixel is the next nearest neighbor to the central pixel as defined by \(r_{el}\). This process goes on until the profile is fully built upto the truncation radius. This is done fairly efficiently using a breadth first parsing strategy227 which is implemented through an ordered linked list.

Using this approach, we build the profile by expanding the circumference. Not one more extra pixel has to be checked (the calculation of \(r_{el}\) from above is not cheap in CPU terms). Another consequence of this strategy is that extending MakeProfiles to three dimensions becomes very simple: only the neighbors of each pixel have to be changed. Everything else after that (when the pixel index and its radial profile have entered the linked list) is the same, no matter the number of dimensions we are dealing with.



Do not confuse the signs of \(sin\) with the rotation matrix defined in Linear warping basics. In that equation, the point is rotated, here the coordinates are rotated and the point is fixed.