Siemens Star Analysis

In this notebook I wil show how to use the Siemens Star Target to characterize the Modulation Transfer Function of an optical system.

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.

How to compute the contrast modulation?

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

Now let's draw these circles around the center

shape = data.shape[0]
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)