Skip to main content

How to use a tilted edge in an imaging target to characterize the MTF of the telescope

Slanted Edge method

In this notebook I wil show how to use an image of a slanted edge to characterize the Modulation Transfer Function of an optical system.

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

    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])
    #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]:
<matplotlib.image.AxesImage at 0x7fc15df97990>

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]:
for row in range(n_rows):
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])