In [ ]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import HTML
HTML('../style/course.css') #apply general CSS
In [ ]:
import matplotlib
from scipy import optimize
import astropy.io.fits
matplotlib.rcParams.update({'font.size': 18})
matplotlib.rcParams.update({'figure.figsize': [12,8]} )
In radio astronomy, source finding is the process through which the attributes of radio sources -- such as flux density and mophorlogy -- are measured from data. In this section we will only cover source finding in the image plane.
Source finding techniques usually involve four steps, i) charecterizing the noise (or background estimation), ii) thresholding the data based on knowledge of the noise, iii) finding regions in the thresholded image with "similar" neighbouring pixels (this is that same as blob detection in image processing), and iv) parameterizing these 'blobs' through a function (usually a 2D Gaussian). The source attributes are then estimated from the parameterization of the blobs.
As mentioned before, the radio data we process with source finders is noisy. To charecterize this noise we need to make a few assumptions about its nature, namely we assume that the niose results from some stochastic process and that it can be described by a normal distribution
$$ G(x \, | \, \mu,\sigma^2) = \frac{1}{\sigma \sqrt{2\pi}}\text{exp}\left( \frac{-(x-\mu)^2}{2\sigma^2}\right) $$
where, $\mu$ is the mean (or expected value) of the variable $x$, and $\sigma^2$ is the variance of the distribution; $\sigma$ is the standard deviation. Hence, the noise can be parameterized through the mean and the standard deviation. Let us illustrate this with an example. Bellow is a noise image from a MeerKAT simulation, along with a histogram of of the pixels (in log space).
In [ ]:
noise_image = "../data/fits/noise_image.fits"
with astropy.io.fits.open(noise_image) as hdu:
data = hdu[0].data[0,0,...]
fig, (image, hist) = plt.subplots(1, 2, figsize=(18,6))
histogram, bins = np.histogram(data.flatten(), bins=401)
dmin = data.min()
dmax = data.max()
x = np.linspace(dmin, dmax, 401)
im = image.imshow(data)
mean = data.mean()
sigma = data.std()
peak = histogram.max()
gauss = lambda x, amp, mean, sigma: amp*np.exp( -(x-mean)**2/(2*sigma**2))
fitdata = gauss(x, peak, mean, sigma)
plt.plot(x, fitdata)
plt.plot(x, histogram, "o")
plt.yscale('log')
plt.ylim(1)
Now, in reality the noise has to measured in the presence of astrophysical emission. Furthermore, radio images are also contaminated by various instrumental effects which can manifest as spurious emission in the image domain. All these factors make it difficult to charercterize the noise in a synthesized image. Since the noise generally dominates the images, the mean and standard deviation of the entire image are still fairly good approximations of the noise. Let us now insert a few sources (image and flux distribution shown below) in the noise image from earlier and then try to estimate noise.
In [ ]:
noise_image = "../data/fits/star_model_image.fits"
with astropy.io.fits.open(noise_image) as hdu:
data = hdu[0].data[0,0,...]
fig, (image, hist) = plt.subplots(1, 2, figsize=(18,6))
histogram, bins = np.histogram(data.flatten(), bins=101)
dmin = data.min()
dmax = data.max()
x = np.linspace(dmin, dmax, 101)
im = image.imshow(data)
mean = data.mean()
sigma_std = data.std()
peak = histogram.max()
gauss = lambda x, amp, mean, sigma: amp*np.exp( -(x-mean)**2/(2*sigma**2))
fitdata_std = gauss(x, peak, mean, sigma_std)
plt.plot(x, fitdata_std, label="STD DEV")
plt.plot(x, histogram, "o", label="Data")
plt.legend(loc=1)
plt.yscale('log')
plt.ylim(1)
The pixel statistics of the image are no longer Gaussian as apparent from the long trail of the flux distribution. Constructing a Gaussian model from the mean and standard deviation results in a poor fit (blue line in the figure on the right). A better method to estimate the variance is to measure the dispersion of the data points about the mean (or median), this is the mean/median absolute deviation (MAD) technique. We will refer to the to median absolute deviation as the MAD Median, and the mean absolute deviation as the MAD Mean. A synthesis imaging specific method to estimate the variance of the noise is to only consider the negative pixels. This works under the assumption that all the astrophysical emission (at least in Stokes I) has a positive flux density. The Figure below shows noise estimates from methods mentioned above.
In [ ]:
mean = data.mean()
sigma_std = data.std()
sigma_neg = data[data<0].std() * 2
mad_mean = lambda a: np.mean( abs(a - np.mean(a) ))
sigma_mad_median = np.median( abs(data - np.median(data) ))
mad_mean = lambda a: np.mean( abs(a - np.mean(a) ))
sigma_mad_mean = mad_mean(data)
peak = histogram.max()
gauss = lambda x, amp, mean, sigma: amp*np.exp( -(x-mean)**2/(2*sigma**2))
fitdata_std = gauss(x, peak, mean, sigma_std)
fitdata_mad_median = gauss(x, peak, mean, sigma_mad_median)
fitdata_mad_mean = gauss(x, peak, mean, sigma_mad_mean)
fitdata_neg = gauss(x, peak, mean, sigma_neg)
plt.plot(x, fitdata_std, label="STD DEV")
plt.plot(x, fitdata_mad_median, label="MAD Median")
plt.plot(x, fitdata_mad_mean, label="MAD Mean")
plt.plot(x, fitdata_neg, label="Negative STD DEV")
plt.plot(x, histogram, "o", label="Data")
plt.legend(loc=1)
plt.yscale('log')
plt.ylim(1)
The MAD and negtive value standard deviation methods produce a better solution to the noise distribution in the presence of sources.
Once the noise has been estimated, the next step is to find and charecterize sources in the image. Generically in image processing this is known as blob detection. In a simple case during synthesis imaging we define a blob as a group contiguous pixels whose spatial intensity profile can be modelled by a 2D Gaussian function. Of course, more advanced functions could be used. Generally, we would like to group together near by pixels, such as spatially 'close' sky model components from deconvolution, into a single complex source. Our interferometric array has finite spatial resolution, so we can further constrain our blobs not to be significantly smaller than the image resolution. We define two further constraints of a blob, the peak and boundary thresholds. The peak threshold, defined as
$$ \sigma_\text{peak} = n * \sigma, $$is the minimum intensity the maximum pixel in a blob must have relative to the image noise. That is, all blobs with peak pixel lower than $\sigma_\text{peak}$ will be excluded from being considered sources. And the boundary threshold
$$ \sigma_\text{boundary} = m * \sigma, $$defines the boundary of a blob, $m$ and $n$ are natural numbers with $m$ < $n$.
We are now in a position to write a simple source finder. To do so we implement the following steps:
In [ ]:
def gauss2D(x, y, amp, mean_x, mean_y, sigma_x, sigma_y):
""" Generate a 2D Gaussian image"""
gx = -(x - mean_x)**2/(2*sigma_x**2)
gy = -(y - mean_y)**2/(2*sigma_y**2)
return amp * np.exp( gx + gy)
def err(p, xx, yy, data):
"""2D Gaussian error function"""
return gauss2D(xx.flatten(), yy.flatten(), *p) - data.flatten()
def fit_gaussian(data, psf_pix):
"""Fit a gaussian to a 2D data set"""
width = data.shape[0]
mean_x, mean_y = width/2, width/2
amp = data.max()
sigma_x, sigma_y = psf_pix, psf_pix
params0 = amp, mean_x, mean_y, sigma_x,sigma_y
npix_x, npix_y = data.shape
x = np.linspace(0, npix_x, npix_x)
y = np.linspace(0, npix_y, npix_y)
xx, yy = np.meshgrid(x, y)
params, pcov, infoDict, errmsg, sucess = optimize.leastsq(err,
params0, args=(xx.flatten(), yy.flatten(),
data.flatten()), full_output=1)
perr = abs(np.diagonal(pcov))**0.5
model = gauss2D(xx, yy, *params)
return params, perr, model
In [ ]:
def source_finder(data, peak, boundary, width, psf_pix):
"""A simple source finding tool"""
# first we make an estimate of the noise. Lets use the MAD mean
sigma_noise = mad_mean(data)
# Use noise estimate to set peak and boundary thresholds
peak_sigma = sigma_noise*peak
boundary_sigma = sigma_noise*boundary
# Pad the image to avoid hitting the edge of the image
pad = width*2
residual = np.pad(data, pad_width=((pad, pad), (pad, pad)), mode="constant")
model = np.zeros(residual.shape)
# Create slice to remove the padding later on
imslice = [slice(pad, -pad), slice(pad,-pad)]
catalog = []
# We will need to convert the fitted sigma values to a width
FWHM = 2*np.sqrt(2*np.log(2))
while True:
# Check if the brightest pixel is at least as bright as the sigma_peak
# Otherwise stop.
max_pix = residual.max()
if max_pix<peak_sigma:
break
xpix, ypix = np.where(residual==max_pix)
xpix = xpix[0] # Get first element
ypix = ypix[0] # Get first element
# Make slice that selects box of size width centred around bright brightest pixel
subim_slice = [ slice(xpix-width/2, xpix+width/2),
slice(ypix-width/2, ypix+width/2) ]
# apply slice to get subimage
subimage = residual[subim_slice]
# blank out pixels below the boundary threshold
mask = subimage > boundary_sigma
# Fit gaussian to submimage
params, perr, _model = fit_gaussian(subimage*mask, psf_pix)
amp, mean_x, mean_y, sigma_x,sigma_y = params
amp_err, mean_x_err, mean_y_err, sigma_x_err, sigma_y_err = perr
# Remember to reposition the source in original image
pos_x = xpix + (width/2 - mean_x) - pad
pos_y = ypix + (width/2 - mean_y) - pad
# Convert sigma values to FWHM lengths
size_x = FWHM*sigma_x
size_y = FWHM*sigma_y
# Add modelled source to model image
model[subim_slice] = _model
# create new source
source = (
amp,
pos_x,
pos_y,
size_x,
size_y
)
# add source to catalogue
catalog.append(source)
# update residual image
residual[subim_slice] -= _model
return catalog, model[imslice], residual[imslice], sigma_noise
Using this source finder we can produce a sky model which contains all 17 sources in our test image from earlier in the section.
In [ ]:
test_image = "../data/fits/star_model_image.fits"
with astropy.io.fits.open(test_image) as hdu:
data = hdu[0].data[0,0,...]
catalog, model, residual, sigma_noise = source_finder(data, 5, 2, 50, 10)
print "Peak_Flux Pix_x Pix_y Size_x Size_y"
for source in catalog:
print " %.4f %.1f %.1f %.2f %.2f"%source
fig, (img, mod, res) = plt.subplots(1, 3, figsize=(24,12))
vmin, vmax = sigma_noise, data.max()
im = img.imshow(data, vmin=vmin, vmax=vmax)
img.set_title("Data")
mod.imshow(model, vmin=vmin, vmax=vmax)
mod.set_title("Model")
res.imshow(residual, vmin=vmin, vmax=vmax)
res.set_title("Residual")
cbar_ax = fig.add_axes([0.92, 0.25, 0.02, 0.5])
fig.colorbar(im, cax=cbar_ax, format="%.2g")
The flux and position on each source varies from the true sky model due to the image noise and distribution. The source finding algorithm we above is heuristic example. It has two major flaws : i) it is capable to handling a situation where two or more sources are close enough to each other that would fall within the same sub-image from which the source parameters are estimated, and ii) the noise in radio images is often non-uniform and 'local' noise estimates are required in order to set thresholds. More advanced source finders are used to work on specific source types such as extended objects and line spectra.
In [ ]: