Posts tagged photometry

Pixelating a 2-D Gaussian with Python

4 Sep 2009

Posted by Anthony in Research

No comments

They're coming thick and fast now.

Here's a Python function to accompany the previous post. It's not maximally efficient, but should make sense...

from scipy import stats
def gaussian_pixel(minxy, maxxy, sigma, meanxy=(0.,0.), norm=None):
    """Return the value of a pixel sampling a 2D Gaussian,
    normalized such that the area under the Gaussian is 1
    (default) or such that the peak is given by norm."""
    x1, y1 = minxy
    x2, y2 = maxxy
    x0, y0 = meanxy
    if norm is None:
        norm = 1. / 2 / math.pi / sigma ** 2
    return norm * 2 * math.pi * sigma ** 2 / (x2 - x1) / (y2 - y1) * (
        (1 - stats.erfc((x2 - x0) / math.sqrt(2) / sigma)) / 2.
        - (1 - stats.erfc((x1 - x0) / math.sqrt(2) / sigma)) / 2.) * (
        (1 - stats.erfc((y2 - y0) / math.sqrt(2) / sigma)) / 2.
        - (1 - stats.erfc((y1 - y0) / math.sqrt(2) / sigma)) / 2.)
, ,

On the normalization of PRFs

4 Sep 2009

Posted by Anthony in Research

1 comment

Yesterday I said that the PRF for a map in Jy/beam (or similar) should be normalized so that that peak is 1. But this is true only for an idealised (not pixelated) PRF, or if the map has infinitesimally small pixels.

If the pixels are larger than infinitesimal, as is generally the case, then the maximum value of the pixelated PRF will be the average value over the pixel, which will be less than 1.

For example, if the PRF is a two-dimensional Gaussian, centred on (x_0, y_0), with standard deviation \sigma, then the value in a pixel with x_1 < x < x_2 and y_1 < y < y_2 will be

A = \frac{1}{(x_2 - x_1)(y_2 - y_1)} \int_{x_1}^{x_2} \int_{y_1}^{y_2} e^{- \left( \frac{(x-x_0)^2 + (y-y_0)^2}{2\sigma^2} \right)} \mathrm{d}x \mathrm{d}y.

which is

A = \frac{2 \pi \sigma^2}{(x_2 - x_1)(y_2 - y_1)} \left( \frac{1}{2} \mathrm{erf} \left( \frac{x_2 - x_0}{\sqrt{2}\sigma} \right) - \frac{1}{2} \mathrm{erf} \left( \frac{x_1 - x_0}{\sqrt{2}\sigma} \right) \right)

\times \left( \frac{1}{2} \mathrm{erf} \left( \frac{y_2 - y_0}{\sqrt{2}\sigma} \right) - \frac{1}{2} \mathrm{erf} \left( \frac{y_1 - y_0}{\sqrt{2}\sigma} \right) \right).

Ugh. Let's make that simpler. For a PRF centred on (0,0), and a pixel (\pm r, \pm r), this is

A = \frac{\pi \sigma^2}{2 r^2} \left( \mathrm{erf} \left( \frac{r}{\sqrt{2}\sigma} \right) \right)^2.

As an example, the fairly-Gaussian beam for the Herschel Space Observatory SPIRE instrument has an FWHM of around 18", which corresponds to a standard deviation of around 18. If we make a Jy/beam map with pixel size 6", then the peak value for a 1 Jy point source in the centre of a pixel will be

A = \frac{\pi (7.64


No big deal really...

,

Estimating the flux of a point source

3 Sep 2009

Posted by Anthony in Research

2 comments

You have a map and you know what a point source looks like. How do you filter the map so that the value of each pixel is now the most likely flux of a point source centred on that pixel? (An isolated point source, to be more precise.)

Easy.

First, find P_i, which is the point response function (PRF), telling you what a point source of flux 1 will look like in the map. This may be normalized so that the peak is 1 (if your map is in Jy/beam or similar), or so that \sum P_i = 1 (if your map is in Jy/pixel or similar). If your map is in MJy/sr ... well, figure it out and add a comment below. Basically, if you normalize your PRF correctly, you won't need to worry about the map units in what follows. Phew.

Now the measured value of each pixel around the point source, d_i, will be

d_i = f P_i + n_i,

where f is the flux of the source and n_i is the noise, drawn from a normal distribution with mean zero and standard deviation \sigma_i.

Now the badness of the fit is measured by the \chi^2, which is given by

\chi^2 = \sum_i \left( \frac{d_i - fP_i}{\sigma_i} \right)^2.

At the maximum likelihood value of the flux, f, \chi^2 will be at a minimum, so

\frac{\mathrm{d}\chi^2}{\mathrm{d}f} = 0.

Hence

\sum_i (-2d_i + 2 f P_i^2) / \sigma_i^2 = 0.

Solving this for f, we find the maximum likelihood solution

f = \frac{\sum_i d_i P_i / \sigma_i^2}{\sum_i P_i^2 / \sigma_i^2}.

Now just do this for each pixel in the map (corresponding to a point source centred on each pixel) and you're done.

Worked example. P_i is 0.5, 1.0 and 0.5, for three adjacent pixels (you'll have realised that the map is in Jy/beam or similar), and d_i is 1, 2 and 1 Jy/beam, for three adjacent pixels, with the same (tiny!) value of \sigma_i for each pixel (in this case, we can ignore the value of \sigma_i in what follows). So the flux at the central pixel is estimated to be

f = (0.5 \times 1 + 1.0 \times 2 + 0.5 \times 1) / (0.5^2 + 1.0^2 + 0.5^2) = 2 \,\mathrm{Jy},

which is no surprise, since the maximum value of the map in Jy/beam is 2 at that position.

This is an example of a matched filter (I haven't read the page, but hopefully including the link will make me look clever). And, given that point sources are under no particular obligation to align themselves with the centres of the pixels of your map, P_i can easily be re-estimated for a source with a certain offset from the pixel centre.

,
flicker_image

Visualizing noisy images

13 May 2009

Posted by Anthony in Computing

No comments

You have an image. Each pixel has a value with some uncertainty. How do you visualize the uncertainty in each pixel? Like this:

flicker_image

Here's the Python code

import numpy as np
from matplotlib import pyplot as plt
 
class FlickerImage(object):
    def __init__(self, im, err):
        self.im = im.copy()
        self.err = err.copy()
        finite = np.isfinite(self.im + self.err)
        self.vmin = (self.im - 2 * self.err)[finite].min()
        self.vmax = (self.im + 2 * self.err)[finite].max()
        self.im[np.invert(finite)] = self.vmax
        self.err[np.invert(finite)] = 0
    def flicker(self):
        fg = plt.imshow(np.zeros(self.im.shape),
                        interpolation='nearest',
                        vmin=self.vmin,
                        vmax=self.vmax)
        while True:
            ran = np.random.normal(size=im.shape)
            fg.set_data(im + err * ran)
            plt.draw()

And here's an example script:

import pyfits
f = pyfits.open('file.fits')
im = f["IMAGE"].data
err = f["ERROR"].data
flicker_image = FlickerImage(im, err)
flicker_image.flicker()
, ,