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:
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]:
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)
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()
Fitting detected circles and deriving their centers¶
In [67]:
from scipy import optimize
def calc_R(x,y, xc, yc):
""" calculate the distance of each 2D points from the center (xc, yc) """
return np.sqrt((x-xc)**2 + (y-yc)**2)
def f(c, x, y):
""" calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc) """
Ri = calc_R(x, y, *c)
return Ri - Ri.mean()
def leastsq_circle(x_m,y_m,x,y):
# coordinates of the barycenter
#x_m = np.mean(x)
#y_m = np.mean(y)
center_estimate = x_m, y_m
#print x_m,y_m
center, pcov = optimize.leastsq(f, center_estimate, args=(x,y))
xc, yc = center
s_sq= np.sum(f(center, x,y)**2)/(len(x)-len(center))
pcov = pcov * s_sq
#err = np.sqrt(np.diag(pcov))
Ri = calc_R(x, y, *center)
R = Ri.mean()
residu = np.sum((Ri - R)**2)
return xc, yc, R
def plot_data_circle(x,y, xc, yc, R):
f = plt.figure( facecolor='white') #figsize=(7, 5.4), dpi=72,
plt.axis('equal')
theta_fit = np.linspace(-np.pi, np.pi, 180)
x_fit = xc + R*np.cos(theta_fit)
y_fit = yc + R*np.sin(theta_fit)
plt.xlabel('x')
plt.ylabel('y')
# plot data
plt.plot(x, y, 'rx', label='data', mew=1)
plt.plot(x_fit, y_fit, 'b-' , label="fitted circle", lw=2)
plt.plot([xc], [yc], 'bD', mec='y', mew=1)
plt.legend(loc='best',labelspacing=0.1 )
plt.grid()
#plt.title('Least Squares Circle')
## loop on features:
for i in np.arange(numb):
y,x = props[i].coords[:,0], props[i].coords[:,1]
y_m, x_m = props[i].centroid
#yc,xc = props[i].centroid
#R_eq = props[i].major_axis_length
#plt.plot(x,y,'rx')
#plt.Circle((xc,yc),R_eq/2,color='blue')
#plt.show()
xc, yc, R = leastsq_circle(x_m,y_m,x,y)
print R
# print R
#print xc, yc
plot_data_circle(x,y,xc,yc,R)
In [71]:
for i in np.arange(numb):
y,x = props[i].coords[:,0], props[i].coords[:,1]
R = props[i].equivalent_diameter/2.
yc, xc = props[i].centroid
plot_data_circle(x,y,xc,yc,R)
In [17]:
center_estimate = x_m, y_m
pfit, pcov = optimize.leastsq(f, center_estimate, args=(x,y))
xc, yc = center
s_sq= np.sum(f(center, x,y)**2)/(len(x)-len(center))
s_sq
#pcov = pcov * s_sq
#pcov
#err = np.sqrt(np.diag(pcov))
In [72]:
'''
F = open('Centers_python_circles','w')
for i in range(numb):
F.write(str(centers[i][1]) + ' '+ str(centers[i][0])+'\n')
F.close()
'''
Out[72]:
Distortion work¶
In [4]:
import os
#p = os.getcwd()
fig = plt.figure(figsize=(15,12))
lines = np.loadtxt('/home/fatima/Desktop/solar_orbiter_project/codes/dont_touch/centers_python_lines_circles')
columns = np.loadtxt('/home/fatima/Desktop/solar_orbiter_project/codes/dont_touch/centers_python_columns_circles')
plt.plot(lines[:,0], lines[:,1],'rx',markersize=8,mew=3)
plt.xlim(0,2048)
plt.ylim(0,2048)
#plt.title('Calibration image',fontsize=22)
#plt.savefig('coordinates_circles.png', dpi=300)
plt.show()
In [6]:
plt.clf()
fig = plt.figure(figsize=(15,12))
image = pyfits.getdata('/home/fatima/Desktop/solar_orbiter_project/codes/targets/calibration_circles.fits')
Ny, Nx = image.shape
i,j = np.indices([Ny,Nx])
center = np.array([(i.max()-i.min())/2.0, (j.max()-j.min())/2.0])
xc = center[1]
yc = center[0]
plt.scatter(xc,yc,marker='x',color='b',s=50,alpha=0.8)
plt.plot(lines[:,0], lines[:,1],'rx',markersize=8,mew=3)
#plt.title('Position of image center',fontsize=22)
plt.xlim(0,2048)
plt.ylim(0,2048)
#plt.plot(columns[:,0], columns[:,1],'bo')
Ny, Nx,xc, yc
#plt.savefig('center_circles.png',dpi=300)
Out[6]:
In [7]:
LINES = []; LINES.append(lines[0:5]); LINES.append(lines[5:10]);LINES.append(lines[10:15]); LINES.append(lines[15:20]);LINES.append(lines[20:25])
In [8]:
LINES
Out[8]:
In [9]:
COLUMNS = []; COLUMNS.append(columns[0:5]); COLUMNS.append(columns[5:10]);COLUMNS.append(columns[10:15]);COLUMNS.append(columns[15:20]); COLUMNS.append(columns[20:25])
In [10]:
len(COLUMNS)
Out[10]:
In [11]:
def parabola_lines(x,a,b,c):
y = a*x**2 + b*x + c
return y
def parabola_columns(y,m,n,p):
x = m*y**2 + n*y + p
return x
In [12]:
plt.clf()
fig=plt.figure(figsize=(20,15))
ax = fig.add_subplot(111)
n=0
A = np.zeros(len(LINES))
B = np.zeros(len(LINES))
C = np.zeros(len(LINES))
for line in LINES:
x = line[:,0]; y = line[:,1]
temp1 = np.array([[x[0]**2,x[1]**2,x[2]**2], [x[0],x[1],x[2]],[1,1,1]])
temp2 = np.array([y[0], y[1],y[2]])
p0 = np.linalg.solve(temp1,temp2)
popt, pcov = curve_fit(parabola_lines, x,y,p0)
## GOODNESS OF FIT
chi_sq_w_lines = np.sum((parabola_lines(x,*popt)-y)**2)
red_chi_sq_lines = chi_sq_w_lines/(len(x)-len(popt))
print chi_sq_w_lines
A[n] = popt[0]; B[n] = popt[1]; C[n] = popt[2]
xnew = np.linspace(100, 1800, 100)
ynew = parabola_lines(xnew, *popt)
ax.plot(x,y,'x',markersize=8,mew=3)
ax.plot(xnew, ynew,label=str(n))
ax.set_xlim(0,2048)
ax.set_ylim(0,2048,2400)
n=n+1
#plt.legend()
plt.savefig('fits_lines.png',dpi=300)
plt.show()
In [13]:
plt.clf()
fig=plt.figure(figsize=(20,15))
ax = fig.add_subplot(111)
n=0
M = np.zeros(len(COLUMNS))
N = np.zeros(len(COLUMNS))
P = np.zeros(len(COLUMNS))
for col in COLUMNS:
x = col[:,0]; y = col[:,1]
temp1 = np.array([[x[0]**2,x[1]**2,x[2]**2], [x[0],x[1],x[2]],[1,1,1]])
temp2 = np.array([y[0], y[1],y[2]])
p0 = np.linalg.solve(temp1,temp2)
popt, pcov = curve_fit(parabola_columns, y,x,p0)
chi_sq_w_cols = np.sum((parabola_columns(y,*popt)-x)**2)
red_chi_sq_cols = chi_sq_w_cols/(len(x)-len(popt))
print chi_sq_w_cols ,red_chi_sq_cols
M[n]=popt[0]; N[n] = popt[1]; P[n] = popt[2]
ynew = np.linspace(100, 1800, 100)
xnew = parabola_columns(ynew, *popt)
ax.plot(x,y,'x',markersize=8,mew=3)
ax.plot(xnew, ynew,label=str(n))
ax.set_xlim(0,2048)
ax.set_ylim(0,2048)
n=n+1
#plt.legend()
plt.savefig('fits_cols.png',dpi=300)
plt.show()
In [14]:
fig=plt.figure(figsize=(20,15))
ax = fig.add_subplot(111)
for col, line in zip(COLUMNS,LINES):
x = col[:,0]; y = col[:,1]
temp1 = np.array([[x[0]**2,x[1]**2,x[2]**2], [x[0],x[1],x[2]],[1,1,1]])
temp2 = np.array([y[0], y[1],y[2]])
p0 = np.linalg.solve(temp1,temp2)
popt, pcov = curve_fit(parabola_columns, y,x,p0)
ynew = np.linspace(0, 2048, 100)
xnew = parabola_columns(ynew, *popt)
#ax.plot(x,y,'x',markersize=8,mew=3)
ax.plot(xnew, ynew,label=str(n))
x = line[:,0]; y = line[:,1]
temp1 = np.array([[x[0]**2,x[1]**2,x[2]**2], [x[0],x[1],x[2]],[1,1,1]])
temp2 = np.array([y[0], y[1],y[2]])
p0 = np.linalg.solve(temp1,temp2)
popt, pcov = curve_fit(parabola_lines, x,y,p0)
xnew = np.linspace(0, 2048, 100)
ynew = parabola_lines(xnew, *popt)
#ax.plot(x,y,'x',markersize=8,mew=3)
ax.plot(xnew, ynew,label=str(n))
ax.set_xlim(0,2048)
ax.set_ylim(0,2048,2400)
n+1
#plt.savefig('grid_circles_dis.png',dpi=300)
plt.show()
Determining the center of curvature¶
In [15]:
fig = plt.figure(figsize=(15,12))
ax = fig.add_subplot(111)
ax.plot(C,A,'rx',label='lines',markersize=8,mew=3)
ax.plot(P,M,'bx',label='columns',markersize=8,mew=3)
ax.set_xlabel('offset from origin',fontsize=22)
ax.set_ylabel('curvature',fontsize=22)
#ax.set_xlim(P.min(), P.max())
#ax.set_ylim(A.min(), A.max())
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
plt.legend(loc='best')
#plt.savefig('cvso_origin.png',dpi=300)
Out[15]:
In [16]:
## Shifting to center of image to find the distortion center
C_2 = A*xc**2 + B*xc + C - yc #lines
P_2 = M*yc**2 + N*yc + P - xc #columns
In [17]:
fig = plt.figure(figsize=(15,12))
ax = fig.add_subplot(111)
ax.plot(C_2,A,'rx',label='lines',markersize=8,mew=3)
ax.plot(P_2,M,'bx',label='columns',markersize=8,mew=3)
ax.set_xlabel('offset from image center',fontsize=22)
ax.set_ylabel('curvature',fontsize=22)
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
plt.legend(loc='best')
#plt.savefig('cvso.png',dpi=300)
Out[17]:
In [18]:
## lines
t = np.polyfit(C_2, A, 1)
tnew = np.linspace(C_2.min(),C_2.max(),100)
model = np.poly1d(t)
Tnew = model(tnew)
## GOODNESS OF FIT
chi_sq_w_lines = np.sum((model(A)-A)**2)
red_chi_sq_lines = chi_sq_w_lines/(len(A)-len(t))
print t
print chi_sq_w_lines,red_chi_sq_lines
## columns
t2 = np.polyfit(P_2, M, 1)
tnew2 = np.linspace(P_2.min(),P_2.max(),100)
model2 = np.poly1d(t2)
Tnew2 = model2(tnew2)
print t2
## goodness of fit
chi_sq_w_cols = np.sum((model2(M)-M)**2)
red_chi_sq_cols = chi_sq_w_cols/(len(M)-len(t2))
print chi_sq_w_cols ,red_chi_sq_cols
fig = plt.figure(figsize=(18,12))
ax = fig.add_subplot(111)
ax.plot(C_2,A,'rx',label='lines',markersize=8,mew=3)
ax.plot(tnew,Tnew,'r')
ax.plot(P_2,M,'bx',label='columns',markersize=8,mew=3)
ax.plot(tnew2, Tnew2, 'b')
ax.set_xlabel('offset from image center',fontsize=22)
ax.set_ylabel('curvature',fontsize=22)
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
plt.legend(loc='best')
#plt.savefig('cvso_fits.png',dpi=300)
Out[18]:
In [19]:
## Finding distortion center:
del_lines = -t[1]/t[0]
del_columns = -t2[1]/t2[0]
y_d = del_lines+yc
x_d = del_columns + xc
print 'the distortion center is:', (x_d,y_d)
del_lines, del_columns
Out[19]:
In [20]:
plt.clf()
fig = plt.figure(figsize=(15,12))
plt.scatter(xc,yc,marker='x', color='m', s=60,label='image center')
plt.scatter(x_d,y_d,marker='x',color='k', s=60,label='distortion center')
#plt.plot(lines[:,0], lines[:,1],'rx')
plt.xlim(0,2048)
plt.ylim(0,2048)
plt.plot(columns[:,0], columns[:,1],'xr', markersize=10)
plt.legend(loc='best', prop={'size': 15})
print xc
Refining the position of the distortion center¶
In [21]:
Ny, Nx = image.shape
i,j = np.indices([Ny,Nx])
center = np.array([(i.max()-i.min())/2.0, (j.max()-j.min())/2.0])
xd = center[1]
yd = center[0]
C_2 = A*xd**2 + B*xd + C - yd #lines
P_2 = M*yd**2 + N*yd + P - xd #columns
trials = 10
i=0
temp_x = temp_y = np.zeros((trials))
In [22]:
for n in range(trials):
t = np.polyfit(C_2, A, 1)
t2 = np.polyfit(P_2, M, 1)
del_lines = -t[1]/t[0]
del_columns = -t2[1]/t2[0]
yd = del_lines+yd
xd = del_columns + xd
print 'the distortion center is:', (xd,yd)
print del_lines, del_columns
C_2 = A*xd**2 + B*xd + C - yd#lines
P_2 = M*yd**2 + N*yd + P - xd #columns
temp_x[i] = xd
temp_y[i] = yd
i=i+1
In [23]:
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('trials')
ax1.set_ylabel('X_c', color=color)
ax1.plot(np.arange(trials), temp_x, 'rx',color=color)
#ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.set_ylabel('Y_c', color=color) # we already handled the x-label with ax1
ax2.plot(np.arange(trials), temp_y,'bx' ,color=color)
#ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
plt.clf()
In [24]:
## The final center of distortion:
x_d = xd
y_d = yd
In [25]:
## parabola coefficient in the new reference of distortion center
C_3 = A*x_d**2 + B*x_d + C - y_d
P_3 = M*y_d**2 + N*y_d + P - x_d
B_3 = 2*A*x_d + B
N_3 = 2*M*y_d + N
In [26]:
t = np.polyfit(C_3, A, 1)
tnew = np.linspace(C_3.min(),C_3.max(),100)
model = np.poly1d(t)
Tnew = model(tnew)
## GOODNESS OF FIT
chi_sq_w_lines = np.sum((model(A)-A)**2)
red_chi_sq_lines = chi_sq_w_lines/(len(A)-len(t))
print t
print chi_sq_w_lines,red_chi_sq_lines
## columns
t2 = np.polyfit(P_3, M, 1)
tnew2 = np.linspace(P_3.min(),P_3.max(),100)
model2 = np.poly1d(t2)
Tnew2 = model2(tnew2)
print t2
## goodness of fit
chi_sq_w_cols = np.sum((model2(M)-M)**2)
red_chi_sq_cols = chi_sq_w_cols/(len(M)-len(t2))
print chi_sq_w_cols ,red_chi_sq_cols
fig = plt.figure(figsize=(18,12))
ax = fig.add_subplot(111)
ax.plot(C_3,A,'rx',label='lines',markersize=8,mew=3)
ax.plot(tnew,Tnew,'r')
ax.plot(P_3,M,'bx',label='columns',markersize=8,mew=3)
ax.plot(tnew2, Tnew2, 'b')
ax.set_xlabel('offset from distortion center',fontsize=22)
ax.set_ylabel('curvature',fontsize=22)
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
plt.legend(loc='best')
plt.savefig('cvso_final.png',dpi=300)
In [27]:
plt.clf()
fig = plt.figure(figsize=(15,12))
plt.scatter(1023,1023,marker='x', color='m', s=60,label='image center')
plt.scatter(x_d,y_d,marker='x',color='k', s=60,label='distortion center')
#plt.plot(lines[:,0], lines[:,1],'rx')
plt.xlim(0,2048)
plt.ylim(0,2048)
plt.plot(columns[:,0], columns[:,1],'xr', markersize=10)
plt.legend(loc='best', prop={'size': 15})
print x_d, y_d
#plt.savefig('center.png',dpi=300)
Correcting for aspect ratio¶
In [28]:
## slope of curvature versus offset in the lines
s_x = t[0]
## slope of curvature versus offset in columns
s_y = t2[0]
## aspect ratio
AR = np.sqrt(np.abs(s_x/s_y))
print AR
In [29]:
## correcting the coefficients:
M_4 = M/AR
N_4 = N_3/AR
P_4 = P_3/AR
A_4 = A*AR**2
B_4 = B_3*AR
C_4 = C_3
In [30]:
plt.figure(figsize=(18,12))
plt.plot(P_4,M_4,'x',label='cols new')
plt.plot(P_3,M,label='cols old')
plt.plot(C_4,A_4,'x',label='lines new')
plt.plot(C_3,A,label='lines old')
plt.legend()
plt.show()
computing the distortion coefficient¶
In [31]:
k_lines = (-A)/(C_3*(3*A*C_3 + 3*B_3**2 + 1))
plt.figure(figsize=(10,5))
plt.plot(C_3, k_lines,'x')
plt.ylabel('distortion coefficient',fontsize=20)
plt.xlabel('offset',fontsize=20)
plt.title('Lines',fontsize=20)
print k_lines.mean()
print np.average(k_lines,weights=np.abs(C_3))
#plt.savefig('k_lines.png')
In [32]:
k_lines_4 = (-A_4)/(C_4*(3*A_4*C_4 + 3*B_4**2 + 1))
plt.figure(figsize=(10,5))
plt.plot(C_4, k_lines_4,'x')
plt.ylabel('distortion coefficient',fontsize=20)
plt.xlabel('offset',fontsize=20)
plt.title('Lines',fontsize=20)
print k_lines_4.mean()
print np.average(k_lines_4,weights=np.abs(C_4))
In [117]:
'''
## removing the outlier:
np.where(k_lines==k_lines.min())
k_lines_2 = np.delete(k_lines,[2])
print k_lines_2, k_lines_2.shape
'''
Out[117]:
In [33]:
plt.figure(figsize=(10,5))
plt.ylabel('distortion coefficient',fontsize=20)
plt.xlabel('offset',fontsize=20)
plt.title('Columns',fontsize=20)
k_cols = (-M)/(P_3*(3*M*P_3+ 3*N_3**2 + 1))
plt.plot(P_3,k_cols,'x')
print k_cols.mean()
print np.average(k_cols,weights=np.abs(P_3))
#plt.savefig('k_cols.png')
In [35]:
plt.figure(figsize=(10,5))
plt.ylabel('distortion coefficient',fontsize=20)
plt.xlabel('offset',fontsize=20)
plt.title('Columns',fontsize=20)
k_cols_4 = (-M_4)/(P_4*(3*M*P_4+ 3*N_4**2 + 1))
plt.plot(P_4,k_cols,'x')
print k_cols_4.mean()
print np.average(k_cols_4,weights=np.abs(P_4))
In [35]:
plt.plot(P_4,k_cols,'x')
plt.plot(C_4, k_lines,'x')
Out[35]:
In [39]:
K = np.array([k_lines.mean() ,k_cols.mean()])
print K.mean()
K = np.array([k_lines_4.mean() ,k_cols_4.mean()])
print K.mean()
K = np.array([np.average(k_lines, weights=np.abs(C_3)), np.average(k_cols, weights=np.abs(P_3))])
print K.mean()
K = np.array([np.average(k_lines_4, weights=np.abs(C_4)), np.average(k_cols_4, weights=np.abs(P_4))])
print K.mean()
In [40]:
K = np.array([np.average(k_lines_4, weights=np.abs(C_4)), np.average(k_cols_4, weights=np.abs(P_4))])
print "the distortion coefficient is k =", K.mean()
In [42]:
x_d = xd
y_d = yd
print x_d,y_d
k = K.mean()
# correcting for distortion
## loading all data points
f = np.loadtxt('/home/fatima/Desktop/solar_orbiter_project/codes/dont_touch/centers_python_lines_circles')
X_d = f[:,0]
Y_d = f[:,1]
In [43]:
X_u = X_d + (X_d - x_d )*(k*((X_d-x_d)**2 + (Y_d-y_d)**2))
Y_u = Y_d + (Y_d - y_d)*(k*((X_d-x_d)**2 + (Y_d-y_d)**2))
In [44]:
fig = plt.figure(figsize=(15,12))
plt.style.use('ggplot')
ax = fig.add_subplot(111)
#plt.scatter(xc,yc,marker='x',color='m',label='image center')
plt.scatter(x_d,y_d,marker='x',color='k',s=150,alpha=1,label='distortion center')
#plt.plot(lines[:,0], lines[:,1],'rx')
ax.plot(X_d, Y_d,'bo', markersize=8 ,label='distorted')
ax.plot(X_u,Y_u, 'rx', markersize=8,mew=2, label='undistorted')
plt.xlim(0,2048)
plt.ylim(0,2048)
plt.xlabel('PIXELS',fontsize=22)
plt.ylabel('PIXELS',fontsize=22)
#rec = plt.Rectangle((0,0), width=Nx,height=Ny, fill='False')
#ax.add_patch(rec)
rms = 4.5
percent = -0.4
k=8e-09
textstr = '\n'.join((
r'RMS error=%.2f' % (rms, ) +' pixels',
r'percentage of Distortion=%.2f' % (percent, )+'%',
r'Distortion coefficient k=%.2E' % (k, )))
text = r'$r_u = r_d \times (1+k\times r_d^2)$'
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
#ax.text(0.05, 0.95, textstr, transform=ax.transAxes, fontsize=14,verticalalignment='top', bbox=props)
#ax.text(0.4, 0.95, text, transform=ax.transAxes, fontsize=20,verticalalignment='top', bbox=props)
plt.legend(loc='upper right',prop={'size': 15} )
#plt.savefig('distortion.png',dpi=300)
Out[44]:
Distortion measure¶
In [45]:
plate_scale = 0.5
d_m_pixels = np.sqrt((1./len(X_u))*np.sum((X_u-X_d)**2+(Y_u-Y_d)**2))
d_m_arc = d_m_pixels*plate_scale
d_m_km = d_m_arc*725
print "the distortion is around:",'\n', d_m_pixels, 'pixels', '\n', d_m_arc, 'arcseconds', '\n',d_m_km, 'km'
In [46]:
np.sqrt((X_u-X_d)**2+(Y_u-Y_d)**2)
Out[46]:
In [50]:
# disortionpercent:
D=((np.sqrt((X_d-x_d)**2+(Y_d-y_d)**2) - np.sqrt((X_u-x_d)**2+(Y_u-y_d)**2))/np.sqrt((X_u-x_d)**2+(Y_u-y_d)**2))*100
print D.mean()
print np.abs(D).max()
In [51]:
np.sqrt(((X_u-X_d)**2+(Y_u-Y_d)**2)).max()
Out[51]:
Stright lines fit according to the computed K¶
In [54]:
plt.clf()
fig=plt.figure(figsize=(18,12))
ax = fig.add_subplot(111)
n=0
for line in LINES:
X_d = line[:,0] -x_d
Y_d = line[:,1] -y_d
X_u = X_d + (X_d)*(k*((X_d)**2 + (Y_d)**2))
Y_u = Y_d + (Y_d )*(k*((X_d)**2 + (Y_d)**2))
params = np.array([B_3[n]*(3*k*C_3[n]**2 +1), C_3[n]*(k*C_3[n]**2+1)])
t_fit = np.polyfit(Y_u,X_u, 1)
print params
print t_fit
xnew = np.linspace(-2000, 2000, 100)
ynew = params[0]*xnew + params[1]
ax.plot(X_d,Y_d,'o',markersize=8,label='distorted')
ax.plot(X_u,Y_u,'x',markersize=10,label='undistorted')
ax.plot(xnew, ynew,label='straight')
ax.set_xlim(-2000,2400)
ax.set_ylim(80-y_d,1800-y_d)
n=n+1
#plt.legend()
plt.savefig('polyfits_lines.png',dpi=300)
plt.show()
In [55]:
plt.clf()
fig=plt.figure(figsize=(18,12))
ax = fig.add_subplot(111)
n=0
for col in COLUMNS:
X_d = col[:,0] - x_d
Y_d = col[:,1] - y_d
X_u = X_d + (X_d)*(k*((X_d)**2 + (Y_d)**2))
Y_u = Y_d + (Y_d )*(k*((X_d)**2 + (Y_d)**2))
params = np.array([N_3[n]*(3*k*P_3[n]**2 +1), P_3[n]*(k*P_3[n]**2+1)])
ynew = np.linspace(-2000,2000, 100)
xnew = params[0]*ynew + params[1]
ax.plot(X_d,Y_d,'o',markersize=8,label='distorted')
ax.plot(X_u,Y_u,'x', markersize=10,label='undistorted')
ax.plot(xnew, ynew,label='straight')
ax.set_xlim(-1000,1000)
ax.set_ylim(-1800,2400)
n=n+1
#plt.legend()
plt.savefig('polyfits_cols.png',dpi=300)
plt.show()
In [77]:
plt.clf()
fig=plt.figure(figsize=(18,12))
ax = fig.add_subplot(111)
n=0
for col,line in zip(COLUMNS,LINES):
X_d = col[:,0] - x_d
Y_d = col[:,1] - y_d
X_u = X_d + (X_d)*(k*((X_d)**2 + (Y_d)**2))
Y_u = Y_d + (Y_d )*(k*((X_d)**2 + (Y_d)**2))
params = np.array([N_3[n]*(3*k*P_3[n]**2 +1), P_3[n]*(k*P_3[n]**2+1)])
ynew = np.linspace(-2000,2000, 100)
xnew = params[0]*ynew + params[1]
#ax.plot(X_d,Y_d,'o',markersize=8,label='distorted')
#ax.plot(X_u,Y_u,'x', markersize=10,label='undistorted')
ax.plot(xnew, ynew,label='straight')
X_d = line[:,0] -x_d
Y_d = line[:,1] -y_d
X_u = X_d + (X_d)*(k*((X_d)**2 + (Y_d)**2))
Y_u = Y_d + (Y_d )*(k*((X_d)**2 + (Y_d)**2))
params = np.array([B_3[n]*(3*k*C_3[n]**2 +1), C_3[n]*(k*C_3[n]**2+1)])
t_fit = np.polyfit(Y_u,X_u, 1)
xnew = np.linspace(-2000, 2000, 100)
ynew = params[0]*xnew + params[1]
#ax.plot(X_d,Y_d,'o',markersize=8,label='distorted')
#ax.plot(X_u,Y_u,'x',markersize=10,label='undistorted')
ax.plot(xnew, ynew,label='straight')
ax.set_xlim(-1023,1023)
ax.set_ylim(-1023,1023)
n=n+1
#plt.legend()
plt.savefig('grid_circles_und.png',dpi=300)
plt.show()
center of distortion is image center¶
In [185]:
Ny, Nx = image.shape
i,j = np.indices([Ny,Nx])
center = np.array([(i.max()-i.min())/2.0, (j.max()-j.min())/2.0])
xc = center[1]
yc = center[0]
In [186]:
C_2 = A*xc**2 + B*xc + C - yc #lines
P_2 = M*yc**2 + N*yc + P - xc #columns
B_2 = 2*A*xc + B
N_2 = 2*M*yc + N
In [187]:
k_lines = (-A)/(C_2*(3*A*C_2 + 3*B_2**2 + 1))
plt.figure(figsize=(10,5))
plt.plot(C_3, k_lines,'x')
plt.ylabel('distortion coefficient',fontsize=20)
plt.xlabel('offset',fontsize=20)
plt.title('Lines',fontsize=20)
print k_lines
print k_lines.mean()
print np.average(k_lines,weights=np.abs(C_2))
#plt.savefig('k_lines.png')
In [189]:
k_cols = (-M)/(P_2*(3*M*P_2+ 3*N_2**2 + 1))
plt.plot(P_2,k_cols,'x')
print k_cols
print k_cols.mean()
print np.average(k_cols,weights=np.abs(P_2))
In [190]:
K = np.array([np.average(k_lines, weights=np.abs(C_2)), np.average(k_cols, weights=np.abs(P_2))])
print "the distortion coefficient is k =", K.mean()
In [192]:
f = np.loadtxt(p+'/dont_touch/centers_python_lines_circles')
X_d = f[:,0]
Y_d = f[:,1]
x_d = xc
y_d =yc
In [193]:
X_u = X_d + (X_d - x_d )*(k*((X_d-x_d)**2 + (Y_d-y_d)**2))
Y_u = Y_d + (Y_d - y_d)*(k*((X_d-x_d)**2 + (Y_d-y_d)**2))
In [194]:
plate_scale = 0.5
d_m_pixels = np.sqrt((1./len(X_u))*np.sum((X_u-X_d)**2+(Y_u-Y_d)**2))
d_m_arc = d_m_pixels*plate_scale
d_m_km = d_m_arc*725
print "the distortion is around:",'\n', d_m_pixels, 'pixels', '\n', d_m_arc, 'arcseconds', '\n',d_m_km, 'km'
In [197]:
R = np.sqrt(((X_u-X_d)**2+(Y_u-Y_d)**2))
print np.where(R==R.max())
print R.max()
print X_u[0], X_d[0]
print Y_u[0], Y_d[0]
D[0]
Out[197]:
equivalence¶
In [253]:
f = np.loadtxt(p+'/dont_touch/centers_python_lines_circles')
X_d = f[:,0]
Y_d = f[:,1]
x_d = 1017.3841738107818
y_d = 906.05325667308955
k = 7.99207263147e-09
P1 = -k*x_d
P2 = -k*y_d
def model_x(params,X_d,Y_d,X_u,Y_u):
xc = params[0]
yc = params[1]
r_d = np.sqrt((X_d-xc)**2 + (Y_d-yc)**2)
x_u = X_d+(X_d-xc)*k*(r_d)**2 +P1*(r_d**2+2*(X_d-xc)**2)+2*P2*(X_d-xc)*(Y_d-yc)
y_u = Y_d + (Y_d - yc)*k*(r_d)**2 +2*P1*(X_d-xc)*(Y_d-yc) + P2*(r_d**2+2*(Y_d-yc)**2)
return np.sum((x_u-X_u)**2)
In [267]:
import scipy
p0 = [1023,1023]
def Con(params):
xc = params[0]
yc = params[1]
r_d = np.sqrt((X_d-xc)**2 + (Y_d-yc)**2)
x_u = X_d+(X_d-xc)*k*(r_d)**2 +P1*(r_d**2+2*(X_d-xc)**2)+2*P2*(X_d-xc)*(Y_d-yc)
y_u = Y_d + (Y_d - yc)*k*(r_d)**2 +2*P1*(X_d-xc)*(Y_d-yc) + P2*(r_d**2+2*(Y_d-yc)**2)
d_m_pixels = np.sqrt((1./len(x_u))*np.sum((x_u-X_d)**2+(y_u-Y_d)**2))
return 5-d_m_pixels
cons = ({'type':'ineq','fun':Con})
mini_x = scipy.optimize.minimize(model_x, p0, args=(X_d,Y_d,X_u,Y_u),constraints=cons,method='SLSQP')
In [268]:
mini_x.x
Out[268]:
In [269]:
xc = mini_x.x[0]
yc = mini_x.x[1]
In [270]:
r_d = np.sqrt((X_d-xc)**2 + (Y_d-yc)**2)
In [271]:
x_u = X_d+(X_d-xc)*k*(r_d)**2 +P1*(r_d**2+2*(X_d-xc)**2)+2*P2*(X_d-xc)*(Y_d-yc)
In [272]:
y_u = Y_d + (Y_d - yc)*k*(r_d)**2 +2*P1*(X_d-xc)*(Y_d-yc) + P2*(r_d**2+2*(Y_d-yc)**2)
In [273]:
d_m_pixels = np.sqrt((1./len(x_u))*np.sum((x_u-X_d)**2+(y_u-Y_d)**2))
In [274]:
d_m_pixels
Out[274]:
In [ ]:
In [ ]: