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 from scipy.signal import savgol_filter
path = '/home/fatima/Desktop/solar_orbiter_project/codes/targets/MTF/'
star = pyfits.getdata(path+'solo_L0_1805230011_20180523T121104.fits') dark = pyfits.getdata(path+'solo_L0_1805230009_20180523T094103.fits') data = star/20. -dark/20.
fig = plt.figure(figsize=(18,12)) plt.imshow(data,origin='lower',cmap='gray')
<matplotlib.image.AxesImage at 0x7f4556699950>
The siemens star target is used in the characterization of the performance of an optical system (from microscopes to telescopes). It consists of alternating black and white stripes who are densely packed near the center of the image and sparsing out towards the edges.¶
The modulation transfer function inspects the contrast variation with spatial frequency. You could get the latter from applying Fourier transform, but in targets like the Siemens star you could measure both the contrast and the spatial frequency. As you could see, the spatial frequency in the siemens star increases from the edges towards the center. Think of the spatial frequency as the inverse of the period which is the distance from white to white or black to black, since this distance increases towards the edges, the spatial frequency decreases. You could also notice how the contrast decreases in the center, since you cannot distinguish black from white at such spatial frequency.¶
The idea is simple, we draw concentric circles around the origin and build the profile along its edges. From this profile (of a frequrency that we know), we compute the contrast modulation. The frequency will be radius dependent and is simply the ratio of the number of line pairs (number of white-black stripes) over the length of the circle. Let's define some functions.¶
## this function will get the values along each circle. We start by defining a range of angles, computing the #coordinates and then finding the values at the nearest pixel position. def get_line(star,theta,radius,x_c,y_c): #theta = np.linspace(0,2*np.pi,N_theta) x = x_c + radius*np.cos(theta) y = y_c + radius*np.sin(theta) x = np.round(x) y = np.round(y) x = x.astype(int) y = y.astype(int) I = star[y,x] return I,x,y
## a function to compute the frequecy for a certain radius def get_radius(freq): N_p = 36 r = N_p/(2*np.pi*freq) return r ## a function to compute the radius for a certain frequency def get_freq(radius): N_p = 36 freq = N_p/(2*np.pi*radius) return freq
shape = data.shape R_max = 1000 N_theta = 300 theta = np.linspace(0,2*np.pi,N_theta) N_radii = 120 N_p = 36 #angular frequency (number of dark (bright) lines) radii = R_max*1/1.05**np.arange(1, N_radii) ## we take a logarithmic range of radii # center of the target (determined manually) x_c = 1089 y_c =996
plt.figure(figsize=(18,12)) plt.imshow(data, origin='lower') for r in radii: freq = get_freq(r) I,x,y = get_line(data,theta,r,x_c,y_c) plt.plot(x,y)