How to use a tilted edge in an imaging target to characterize the MTF of the telescope
The basics of this method is to find a slightly tilted edge in a target which you image with your optical system. The Edge Profile Function (EPF) will inform you on the change of contrast across this edge (by computing the gradient along the edge, which we will call the Line Transfer Function, LTF). The Fourier Transform of the LTF will be the 1D Modulation Transfer Function (MTF). The reason we take a tilted edge is to have more sample data points after projecting them to the direction normal to the edge for a precise computation of the MTF.¶
In [14]:
import numpy as np
from scipy import special
import pyfits
import matplotlib.pyplot as plt
import scipy
import math
from numpy import fft
from numpy.fft import fft, ifft, fftshift, ifftshift
from scipy.optimize import curve_fit
import statsmodels.api as sm
lowess = sm.nonparametric.lowess
In [2]:
## Define Erfc function
def Erfc(x,sigma,a):
y = a*special.erfc(-x/(sigma*np.sqrt(2)))
return y
In [3]:
# Compute the first derivative
def process_limb(x,y):
frac = 1/30.
d1 = np.diff(y)/np.diff(x)
xd1 = 0.5*(x[:-1]+x[1:])
xd2 = 0.5*(xd1[:-1]+xd1[1:])
d2 = np.diff(d1)/np.diff(xd1)
w = lowess(d1, xd1, frac)
s_x, s_y = w[:,0], w[:,1]
#ind = np.where(d1==d1.max())
#print x[ind]
ind = np.where(s_y==s_y.max())
#temp = x[ind]
temp = s_x[ind]
x = (x- temp)#x[ind]
print temp
plt.plot(xd1,d1)
plt.plot(s_x,s_y)
plt.axvline(x=temp)
plt.show()
return temp,x#, y
In [4]:
def detect_edge(x,y):
frac = 1/30.
d1 = np.diff(y)/np.diff(x)
xd1 = 0.5*(x[:-1]+x[1:])
w = lowess(d1, xd1, frac)
s_x, s_y = w[:,0], w[:,1]
ind = np.where(s_y==s_y.max())
plt.plot(x,y)
plt.plot(s_x,s_y)
plt.axvline(x[ind])
plt.show()
#plt.clf()
#print ind
return s_x[ind],y[ind]
In [5]:
def merge(x,y):
x = np.concatenate(x)
y = np.concatenate(y)
return x,y
In [6]:
def sort(x,y):
points = zip(x,y)
points.sort()
xx = np.array([p[0] for p in points])
yy = np.array([p[1] for p in points])
return xx,yy
In [7]:
def compute_lsf(x,y):
d1 = np.abs(np.diff(y)/np.diff(x))
xd1 = 0.5*(x[:-1]+x[1:])
return xd1, d1
In [8]:
def compute_mtf(lsf_x,lsf_y):
lsf = fftshift(lsf_y)
mtf = fft(lsf)
f = (lsf_x.size) * np.arange(0, lsf_x.size, step = 1)
#freq = np.linspace(0,N*0.5,size/2)
#freq = np.arange(0,N*0.5,1./(size))
mtf= np.abs(mtf)
mtf = mtf/mtf.max()
plt.plot(f[f <= 1 / 2], mtf[f <= 1 / 2])
plt.show()
#return freq, mtf[0:size/2]
In [37]:
def Lowess(y,x,frac):
w = sm.nonparametric.lowess(y, x, frac)
s_x, s_y = w[:,0], w[:,1]
return s_x,s_y
In [9]:
path = '/home/fatima/Desktop/solar_orbiter_project/codes/targets/MTF/'
edge = pyfits.getdata(path+'Edge_0_FS_3.fits')
theta = 3
In [10]:
plt.figure(figsize=(8,10))
plt.imshow(edge,origin='lower',cmap='gray')
Out[10]:
I will first compute the edge angle by plotting the edge profiles and computing the inflection point.¶
In [11]:
rows = np.arange(edge.shape[0])
n_rows = len(rows)
i = 0
profile = np.zeros((n_rows,edge.shape[1]))
for row in rows:
profile[i] = edge[row,:]
i = i+1
In [12]:
plt.figure(figsize=(15,12))
for row in range(n_rows):
plt.plot(profile[row],'x')
In [15]:
coords_x = []
coords_y = []
x = np.arange(edge.shape[1])
for row in range(n_rows):
x_e,y_e = detect_edge(x,profile[row])
coords_x.append(x_e)
coords_y.append(row)
In [17]:
x,y = sort(coords_x,coords_y)
In [19]:
from scipy.optimize import curve_fit
x = np.asarray(x)
y = np.asarray(y)
def f(x, A, B): # this is your 'straight line' y=f(x)
return A*x + B
p0 = [1,1]
popt, pcov = scipy.optimize.curve_fit(f,x[:,0],y,p0)#bin_edges[:-1],bin_means) #
In [25]:
angle =np.rad2deg(np.arctan(popt[0])) #np.rad2deg(np.arctan2(ynew[-1] - ynew[0], xnew[-1] - xnew[0]))
print 'the angle is', 90-angle, 'deg'
In [21]:
plt.figure(figsize=(8,10))
plt.imshow(edge,origin='lower',cmap='gray')
plt.plot(x[:,0],y,'rx')
Out[21]:
It is recommended to take the number of rows equal to¶
In [26]:
1/np.tan(math.radians(theta))
Out[26]:
Now we should shift the profiles so that the position of the inflection point is zero¶
In [27]:
rows = np.arange(30,30+19,1)#
n_rows = len(rows) #ed.shape[0]
i = 0
profile = np.zeros((n_rows,edge.shape[1]))
for row in rows:
profile[i] = edge[row,:]
i = i+1
In [28]:
plt.figure(figsize=(18,15))
for row in range(n_rows):
plt.plot(profile[row],'x')
In [29]:
x = np.arange(edge.shape[1])
shifted_x =[]
Temp = []
temp = process_limb(x,profile[0])[0]
shifted_x.append(process_limb(x,profile[0])[1])
Temp.append(temp)
In [30]:
#plt.figure(figsize=(10,10))
#shifted_x = np.zeros((n_rows,edge.shape[1]))
for row in range(1,n_rows):
#shifted_x[row] = process_limb(x,profile[row])
x = np.arange(edge.shape[1])
temp =temp + np.tan(math.radians(theta)) #temp + np.tan(math.radians(theta))
x = x - temp
print temp
shifted_x.append(x)
Temp.append(temp)
#plt.axvline(x=temp)
print temp
#plt.savefig('edge_coords.png',dpi=250)
transforming the coordinate ROF to that normal to the inclined edge¶
In [31]:
shifted_x = np.asarray(shifted_x)
In [32]:
plt.figure(figsize=(10,10))
for row in range(n_rows):
plt.plot(shifted_x[row],profile[row],'x')
In [33]:
shifted_x = shifted_x*np.cos(math.radians(theta))
In [34]:
all_x, all_y = merge(shifted_x,profile)
In [35]:
all_x, all_y = sort(all_x,all_y)
Smoothing the Edge funciton¶
In [38]:
s_x,s_y = Lowess(all_y,all_x,1/40.)
In [39]:
esfx =s_x#all_x
esfy =s_y
In [40]:
plt.figure(figsize=(10,10))
plt.plot(all_x,all_y,'x')
plt.plot(esfx,esfy)
Out[40]:
computing the LSF¶
In [41]:
lsfx,lsfy = compute_lsf(esfx,esfy)
bringing the boundaries of the LSF to zero¶
In [42]:
lsfy2 = lsfy*np.hanning(len(lsfy))
In [43]:
plt.figure(figsize=(10,10))
plt.plot(lsfx,lsfy)
plt.plot(lsfx,lsfy2)
Out[43]:
Computing the MTF¶
In [47]:
plt.figure(figsize=(10,10))
size = lsfx.shape[0]#edge.shape[1]
size2 = edge.shape[1]
lsf = fftshift(lsfy)
mtf = fft(lsf)
f_s = 1./19
f = (1./(f_s*size) * np.arange(0, size, step = 1))
f2 = (1./(size2) * np.arange(0, size2, step = 1))
mtf= np.abs(mtf)
mtf = mtf/mtf.max()
mtf2 = mtf[:edge.shape[1]]
plt.plot(f[f <= n_rows / 2], mtf[f <= n_rows / 2])
plt.scatter(f2, mtf2,label='Slanted edge')
plt.xlim(0,0.5)
#plt.plot(X,Y,label='Siemens star')
plt.legend()
#plt.savefig('MTF.png',dpi=200)
Out[47]:
In [ ]: