mtf.py

# based on https://gist.github.com/stefanv/2051954
from __future__ import division

import numpy as np
from pylab import *
from scipy import ndimage
import time

ion()

def MTF(x, window=True):
    """
    Compute the MTF of an edge scan function.

    Parameters
    ----------
    x : 1D ndarray
        Edge scan function.
    window : bool
        Whether to apply Hanning windowing to the input.

    Notes
    -----
    The line spread function is the derivative of the edge scan function.  The
    FFT of the line spread function gives the MTF.

    See Also
    --------
    http://www.cis.rit.edu/research/thesis/bs/2001/perry/thesis.html

    """
    y = np.diff(x)

    if window:
        y = y * np.hanning(len(y))

    y = np.append(y, np.zeros(100))
    Y = np.fft.fft(y)

    return Y[:len(Y) // 2]
    
N = 100
# Generate edge
x = np.zeros(N)
x[:N // 2] = 1

plt.figure()
for (std,box) in zip(np.arange(1,10,0.1),np.arange(1,10,0.1)):
    plt.clf()

    # Pass through various filters
    y1 = ndimage.gaussian_filter1d(x, std)[box:-box]
    y2 = np.convolve(x, 1/box * np.ones(box), mode='same')[box:-box]

    Y1 = MTF(y1)
    Y2 = MTF(y2)

    plt.subplot(121)
    ix = np.arange(len(Y1)) / (2 * len(Y1))

    plt.plot(y1,label='Gaussian')
    plt.plot(y2,label='Box')
    plt.legend(loc='best')
    plt.xlim([0,100])
    plt.ylim([0,1.1])  

    plt.subplot(122)
    plt.plot(ix, np.abs(Y1), label='Gaussian %0.2f' % std)
    plt.plot(ix, np.abs(Y2), label='Box %0.2f' % box)
    plt.legend()
    plt.ylim([0,1])
    plt.ylim([0,1.1])
    
    
    plt.draw()
    time.sleep(0.0125)
    
ioff()
plt.show()