In [ ]:
from astropy.io import fits
from astropy.stats import sigma_clip
import numpy as np
from scipy import ndimage
In [ ]:
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
%matplotlib inline
# Plot image pixels in cartesian ordering (i.e. y-positive == upwards):
plt.rcParams['image.origin'] = 'lower'
# Make plots bigger
plt.rcParams['figure.figsize'] = 10, 10
Load the data. A mask can be applied if necessary - this may be useful e.g. for excluding the region around a bright source, to avoid false detections due to sidelobes.
In [ ]:
# fits_path = '../casapy-simulation-scripts/simulation_output/vla.image.fits'
# hdu0 = fits.open(fits_path)[0]
# img_data = hdu0.data.squeeze()
# imgdata = np.ma.MaskedArray(imgdata, mask=np.zeros_like(imgdata))
# imgdata.mask[900:1100,900:1100] = True
# imgdata.mask.any()
Alternatively, we can simulate a rudimentary image by adding a Gaussian model to background noise, although note that the noise will be uncorrelated, unlike a radio-synthesis image:
In [ ]:
import fastimgproto.fixtures.image as fixture
img_shape = (128,192)
img_data = fixture.uncorrelated_gaussian_noise_background(img_shape,sigma=0.1)
srcs = []
srcs.append(fixture.gaussian_point_source(x_centre=32.5, y_centre=32.66, amplitude=0.5))
srcs.append(fixture.gaussian_point_source(x_centre=64.12, y_centre=48.88, amplitude=1.0))
srcs.append(fixture.gaussian_point_source(x_centre=128.43, y_centre=94.5, amplitude=1.5))
for s in srcs:
fixture.add_gaussian2d_to_image(s, img_data)
In [ ]:
imgmax = np.max(img_data)
plt.imshow(img_data,vmax=imgmax*0.5)
plt.colorbar()
In [ ]:
from fastimgproto.sourcefind.image import SourceFindImage
sfimage = SourceFindImage(data = img_data, detection_n_sigma=5, analysis_n_sigma=3)
Background level and RMS are crudely estimated via median and sigma-clipped std. dev., respectively:
In [ ]:
sfimage.bg_level
In [ ]:
sfimage.rms_est
We use two thresholds when identifying our source 'islands' (connected pixel regions). The high threshold is our detection level, and should be set high enough to avoid false detections due to noise spikes. The lower threshold expands each island, such that it is more likely to contain enough pixels to reasonably fit a Gaussian profile (otherwise the island may consist of only a single pixel over the detection threshold).
Note that this thresholding approach may result in multi-peaked regions (e.g. two distinct but adjacent sources) being assigned to a single island / label. This can be tackled with 'deblending' algorithms if desired, but is not covered in this notebook.
The thresholded data is then run through scipy.ndimage.label which numbers the connected regions:
In [ ]:
plt.imshow(sfimage.label_map,cmap='Paired')
Plotting all data which is merely above the analysis threshold shows the importance of a usefully high detection threshold - there are many noise spikes above the analysis threshold.
In [ ]:
plt.imshow(sfimage.data > sfimage.analysis_n_sigma*sfimage.rms_est, cmap='binary')
Each island is then analysed for peak value, barycentre, etc (and in may be model-fitted in future):
In [ ]:
island = sfimage.islands[1]
island
In [ ]:
plt.imshow(island.data)
plt.xlim(island.extremum.index.x-10,island.extremum.index.x+10,)
plt.ylim(island.extremum.index.y-10,island.extremum.index.y+10,)
moments_fit = island.params.moments_fit
plt.scatter(moments_fit.x_centre,moments_fit.y_centre, marker='*', s=200, c='y',)
In [ ]:
print("Bright source model:")
print(srcs[-1])
print()
print("Island barycentres:")
for i in sfimage.islands:
moments = i.params.moments_fit
print(moments.x_centre, moments.y_centre)