Skip to main content

Fatima Kahil

How to use the circular grid calibration target to estimate the geometrical distortion in the image plane of a telescope

Fatima Kahil

Using the circles calibration image

In [2]:
import numpy as np
import pyfits
import matplotlib.pyplot as plt
import skimage
from skimage.feature import blob_dog, blob_doh, blob_log, canny
from skimage.color import rgb2gray
from skimage.feature import corner_harris, corner_subpix, corner_peaks
from skimage.segmentation import slic
from skimage.filters import sobel
from scipy.signal import convolve2d
from scipy.ndimage import gaussian_filter
from skimage import measure
from scipy.optimize import curve_fit
import matplotlib.ticker as mtick

Determining the circles centers coordinates

In [3]:
circles = pyfits.getdata('calibration_circles.fits')
In [4]:
plt.figure(figsize=(18,12))

plt.imshow(circles, origin='lower')
plt.savefig('circles',dpi=300)
In [5]:
circles = circles/circles.mean()
In [6]:
hdu = pyfits.PrimaryHDU(circles)
#hdu.writeto('circles_norm.fits')
In [7]:
mask = np.zeros((2048,2048))
for i in range(2048):
    for j in range(2048):
        if circles[i][j] < 0.6:
            mask[i][j] = 1
        else:
            mask[i][j] = 0
In [8]:
hdu = pyfits.PrimaryHDU(mask)
#hdu.writeto('circles_mask.fits')
In [9]:
plt.figure(figsize=(18,12))

plt.imshow(mask,cmap='gray', origin='lower')
Out[9]:
<matplotlib.image.AxesImage at 0x7f7251ede7d0>
In [10]:
from scipy.ndimage.filters import gaussian_filter as gf
mask_smooth = gf(mask,sigma=5,mode='nearest')
mask_smooth.max()
fig = plt.figure(figsize=(18,12))
plt.imshow(mask_smooth,cmap='gray',origin='lower')
hdu = pyfits.PrimaryHDU(mask_smooth)
#hdu.writeto('mask_smooth_circles.fits')
In [11]:
mask_2 = np.zeros((2048,2048))
for i in range(2048):
         for j in range(2048):
             if mask_smooth[i][j] > 0.6:

                   mask_2[i][j] = 1
             else:
                   mask_2[i][j] = 0
In [12]:
fig = plt.figure(figsize=(18,12))
plt.imshow(mask_2,cmap='gray',origin='lower')
hdu = pyfits.PrimaryHDU(mask_2)
plt.savefig('binary_circles.png',dpi=300)
#hdu.writeto('extended_circles.fits')

Labeling features using the sequential region labeling in skimage

In [16]:
from skimage import measure
from skimage.measure import label
from skimage.measure import regionprops

labels, numb = measure.label(mask_2,background=0, return_num=True)
fig = plt.figure(figsize=(18,12))

plt.imshow(labels, origin='lower')
plt.savefig('labels_circles.png',dpi=300)
In [17]:
props = skimage.measure.regionprops(labels)
centers = np.array([prop.centroid for prop in props])
eccen = np.array([prop.eccentricity for prop in props])
R_eq =  np.array([prop.equivalent_diameter for prop in props])
moments_central = np.array([prop.moments for prop in props])
diameter = np.array([prop.equivalent_diameter for prop in props])
area = np.array([prop.area for prop in props])
major =  np.array([prop.major_axis_length for prop in props])
minor = np.array([prop.minor_axis_length for prop in props])
print major/minor
#print eccen
#print moments_central[0].shape
#print diameter
#print 2*np.sqrt(area/np.pi)

[ 1.81742676  1.83985983  1.02236881  1.01181635  1.00727974  1.01200482
  1.02189062         inf  1.01274483  1.00441373  1.00111527  1.00625171
  1.01528526  1.0114089   1.00418396  1.00320184  1.0050895   1.01204674
  1.01186028  1.00440516  1.00105013  1.00096328  1.00898163  1.01900053
  1.00767101  1.00753621  1.0071054   1.01505011  1.63299316]

/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in divide
  # Remove the CWD from sys.path while we load stuff.
In [14]:
for i in np.arange(numb):
    y, x = props[i].coords[:,0], props[i].coords[:,1]
    yc, xc = props[i].centroid
    plt.plot(x,y,'bo')
    plt.plot(xc,yc,'rx')
    #plt.title('Feature# '+str(label),fontsize=22)
    plt.xlabel('Pixels')
    plt.ylabel('Pixels')
    plt.title('xc = '+str(xc)+', yc = '+str(yc))
    plt.show()
#plt.show()