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

## 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:
else:
```
In [8]:
```hdu = pyfits.PrimaryHDU(mask)
```
In [9]:
```plt.figure(figsize=(18,12))

```
Out[9]:
`<matplotlib.image.AxesImage at 0x7f7251ede7d0>`
In [10]:
```from scipy.ndimage.filters import gaussian_filter as gf
fig = plt.figure(figsize=(18,12))
```
In [11]:
```mask_2 = np.zeros((2048,2048))
for i in range(2048):
for j in range(2048):

else:
```
In [12]:
```fig = plt.figure(figsize=(18,12))
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

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()
```