Reference: The Ultimate Resolved Source Cheatsheet

2011 December 7

OK, this isn’t actually the *ultimate* resolved source cheatsheet … yet. I’m
just trying to centralize my notes on dealing with extended sources. I’ve been
doing this lately, and I had to rederive a bunch of things that I know I’ve
figured out before. Hopefully, next time I’ll remember that I wrote all this
stuff down here.

# Units🔗

Radio astronomical images typically come in two unit systems: Jy/pixel and Jy/beam. For various questions one usually wants to convert to Jy/arcsec² or Jy/sr. _Jy/pixel is not trivially equivalent to Jy/arcsec², because the pixel size can change as a function of position in many map projections.

There must be something profound about why maps come in Jy/bm and models come in Jy/pixel, but I don’t see it right now.

Source fluxes are best measured as *total fluxes* which come in units of plain
*Jy*. This should be obvious! Flux is conserved! In the general case *peak*
fluxes are always dependent on the image properties.

The important exception, however, is *unresolved* sources which by their
nature have a peak flux equal to their total flux no matter what unit system
you’re using. Total flux of 1 Jy, peak flux of 1 Jy/px spread over 1 pixel, or
peak flux of 1 Jy/bm spread out over 1 beam.

*Resolved* sources will have a peak flux that depends on the particular image
units and projection, because the total flux will be spread out over a number
of pixels that depends on the image particulars. In a Jy/px image, you need to
figure out how many pixels the source subtends, and *mutatis mutandis* for
Jy/bm. Generically, `pkflux = totflux * arcsec2perpix / arcsec2insource`

. If
you have a 2D Gaussian source where the axes are FWHMs in arcsec,
`arcsec2insource = 2 pi major minor / (8 ln 2)`

.

I *think* the best way to work with Jy/bm is to go through Jy/px: the beam
volume in pixels can be derived from the beam parameters and the pixel volume:
`pixperbm = arcsec2perbm / arcsec2perpix`

. For various operation I feel like
this quantity will cancel out, though. Once again, due to the
position-dependent pixel volume in most maps, this parameter is not constant
across the map.

# 2D Gaussian Parametrizations🔗

There are a few ways to parametrize a 2D Gaussian. I’ll omit centering and overall scaling in the following equations for simplicity.

**Tersest.** The simplest expression in
terms of symbols is:

This parametrization is computationally efficient but not intuitive. One
particular danger to note is that the parameters behave a bit oddly, and there
are domain limitations that make it awkward for numerical optimizers (which
are likely to tweak the parameters outside of their domains): *A* and *C* must
be negative, and, if I’m doing my algebra right, you must have `B² < 4AC`

.

**Mathematical Major/Minor/PA.** We often want to express the shape of the 2D Gaussian as a
major axis σ1, minor axis σ2, and position angle θ. Here σ1 ≥ σ2 > 0 and θ is
the rotation between the major axis and +x (with +θ being towards +y). Because
of the symmetry of the shape, there’s a 180° degeneracy in θ. We find

,

, and

. Note that if θ = 0, then

. Inverting, we find

,

, and

. During a roundtrip, θ may flip by 180° depending on the numerics, but as mentioned above this doesn’t change anything.

**Astronomical Major/Minor/PA.** This is the same as the above, except that by
the definition of PA, the *x* axis is *declination* and the *y* axis is *right
ascension*. Be careful to write

.

**Numerically-friendly.** If you’re fitting for
Gaussian parameters, it’s very helpful to numerical stability to have nice
continuous behavior over all of your parameters. The *B* parameters above seem
to be problematic, so here’s a parametrization that is pretty cheap to compute
and behaves nicely:

,

, and

. It shouldn’t be hard to see that

and

. Here, the optimizer may end up with either of the σ values being the larger. If σ2 > σ1, you can exchange them and add (or subtract) π/2 to θ.

**Probabalistic.** Change gears for a second. What if we interpret our 2D
Gaussian as expressing the joint behavior of variables x and y, with standard
deviations σx and σy and covariance Cxy? In this case, we still ignore mean
values, but there is one correct normalization:

where

These unpleasant expressions are difficult to invert. If you’re trying to
solve for a covariance matrix from some data, use
`numpy.cov`

or read
Estimation of covariance matrices.
It is surprising to me how difficult it is to relate σx, σy, and Cxy to other
quantities — everything seems to be done best by going through σ1, σ2, and θ
(see `ellpar`

below).

**Alternate numerically-friendly probabalistic.** The numerically-friendly
version above still requires p, q > 0 and is periodic with θ. If the above
probabalistic approach is workable, use

. For initialization,

**Wikipedia.** The
Wikipedia page on Gaussian Functions
uses:

Here the condition on the parameters is that the matrix `[a b \ b c]`

be
positive definite, or that

(or, equivalently I think, that *a* and *c* be positive and _b_² < *ac*).

Wikipedia gives a conversion formula to this system from major/minor/PA, but
there are two caveats: the Wikipedia system names σ{1,2} as σ{x,y} somewhat
misleadingly, and inverts θ to be a *clockwise* rotation from +x to -y. Avoid
those equations since the ones I give above are clearer.

# Gaussian Convolutions🔗

Resolved sources are often modeled as 2D Gaussians. These are convolved with
the 2D Gaussian of your synthesized beam, so one often wants to compute this
convolution or deconvolution analytically. The file
`src/subs/gaupar.for`

in MIRIAD has routines for going both directions, and they are not trivial.
The deconvolution case needs special handling for sources that are close to
(or smaller than!) point sources. The existing heuristics seem only so-so. See
`pwpy/scilib/astutil.py`

for a Python deconvolution implementation.

To manage flux through such a deconvolution, work in total
flux units. Nothing changes! The *peak* flux is scaled in two ways: once, to
convert from Jy/bm to Jy/px (assuming this is what you’re doing), which
involves the beam volume and the pixel volume; and once to convert the source
area, which involves the original volume and the deconvolved volume.

# Gaussian Coordinate Transformations🔗

When you have a Gaussian model for a source, this usally needs to be rendered into actual pixels or something.

Doing this purely geometrically is difficult. For instance, in a given map
projection, the RA/Dec axes may not be precisely aligned with the X/Y pixel
axes. If you want to draw a Gaussian source as an ellipse, you need to rotate
it a little bit. If the source is very large, these effects change over the
source extent and things get scary. *It’s better to calculate the RA/Dec of
each pixel and brute-force your way through many of these issues.*

To do this brute forcing, you need to find the value of the parameter that goes in the exponential. You can think of this as computing a distance in from the origin in the transformed space where the Gaussian is a circle. You un-translate, then un-rotate, then un-scale:

```
def gaussdist(x0, y0, maj, min, pa, x, y):
# x0, y0: Gaussian center
# maj, min: major/minor axes
# pa: position angle, from +x toward +y
dx, dy = x - x0, y - y0
c, s = np.cos(pa), np.sin(pa)
a = c * dx + s * dy
b = -s * dx + c * dy
d2 = (a / maj)**2 + (b / min)**2
# d2: squared "distance" from Gaussian such
# that z = exp(-0.5 * d2)
return d2
```

The astronomical “East from North” PA
convention means that you can replace *X* with *declination* and *Y* with
*right ascension*. This formulation can be turned into the one in the
Wikipedia article on Gaussian Function
but is a lot simpler to think about.

# A: Converting Jy/pixel to Jy/arcsec²🔗

To do this you need to calculate the size
of a pixel at the position of the source you’re considering. Code to do this
is in [`pwpy/intfbin/msimgen`

](https://github.com/pkgw/pwpy/blob/master/i
ntfbin/msimgen). It looks something like this:

```
def pixelvolume(pixel2world, pixelcoords):
delta = 1e-5
w1 = pixel2world(pixelcoords)
pixelcoords[X] += delta
pixelcoords[Y] += delta
w2 = pixel2world(pixelcoords)
dra = w2[RA] - w1[RA]
ddec = w2[DEC] - w2[DEC]
return (dra**2 + ddec**2) / (2 * delta**2)
```

# B: Error Ellipses From Covariant Variables🔗

This is almost entirely unrelated to everything else here, except that it involves 2D Gaussians. Given two correlated variables, what are the parameters of the ellipse that describes a confidence region?

```
def ellpar(sx, sy, cxy):
# sx: std dev (not variance) of x var
# sy: std dev (not variance) of y var
# cxy: covariance (not corr. coeff.) of x and y
from numpy import arctan2, sqrt
pa = 0.5 * arctan2(2 * cxy, sx**2 - sy**2)
h = sqrt((sx**2 - sy**2)**2 + 4*cxy**2)
maj = sqrt(2 * (sx**2 * sy**2 - cxy**2) / (sx**2 + sy**2 - h))
min = sqrt(2 * (sx**2 * sy**2 - cxy**2) / (sx**2 + sy**2 + h))
# maj/min: major/minor axes of ellipse
# note: sigmas, not FWHMs
# pa: position angle of ellipse, rotating from +x to +y
return maj, min, pa
```

Say we want to draw an actual ellipse representing this bivariate
distribution, so that the area within the ellipse represents a certain
confidence interval on predicted values. Parametrize this with a value f(CL),
where we draw an ellipse that has axes of sizes f*maj and f*min. By
transforming X and Y to have a mean of 0 and stddev of 1, then considering the
probability distribution of r = x²+y², we can find the relationship between
the axis lengths and the confidence interval they subsume. P(r) = r exp
(-r²/2) so

We see that if we just draw the ellipse according to values associated with the underlying distribution, we get the 39% confidence interval. If we want one in terms of the usual σ limits, we must scale by

For 1 sigma this about 1.5 and for 2 sigma it’s about 2.5. **Note:** Earlier I
had this analysis all wrong!

# C: Tracing Ellipses🔗

To trace an
ellipse, all we need to do is express it in parametric form. Let `th`

be an
offset angle from the major axis of the ellipse. In code form:

```
def ellpoint(th, x0, y0, maj, min, pa):
from numpy import cos, sin
x = x0 + maj * cos(th) * cos(pa) - min * sin(th) * sin(pa)
y = y0 + maj * cos(th) * sin(pa) + min * sin(th) * cos(pa)
return x, y
```

That’s all there is to it.