Generate MTF curves from razorblade images

mkplot.py

from numpy import fft, diff, empty, interp, arange
from imageio import imread
from matplotlib import pyplot as plt

from argparse import ArgumentParser

DPI = 6400
DPMM = DPI / 25.4
SAMPLE_SPACING_MM = 1.0 / DPMM
plt.figure(figsize=(8,6), dpi=100, constrained_layout=True)

parser = ArgumentParser(description = "Generate MTF curves from images")
parser.add_argument('filename')
parser.add_argument('-z', '--horizontal', help='whether the image is horizontal', action='store_true')
parser.add_argument('-t', '--title', help='graph title')
parser.add_argument('-y', '--yvsmtf', help='plot y vs mtf at 10/20/30 lp/mm instead of the mtf function', action='store_true')
args = parser.parse_args()

fname = args.filename
plot_fname = fname.replace('.tiff', '_plot.png')
plot_fname = plot_fname.replace('.tif', '_plot.png')

data = imread(fname)[:,:].astype(float)

# if the *blade edge* runs vertically
if args.horizontal:
    data = data.transpose()

height = data.shape[0]

xs = fft.rfftfreq(data.shape[1]-1, SAMPLE_SPACING_MM)

# indices into the arrays for 10, 20, and 30 lp/mm
x10lpmm_idx = round(interp(10, xs, range(0,len(xs))))
x20lpmm_idx = round(interp(20, xs, range(0,len(xs))))
x30lpmm_idx = round(interp(30, xs, range(0,len(xs))))

if args.yvsmtf:
    xs = arange(0, height) * SAMPLE_SPACING_MM

# Only plot up to ~60 lp/mm, it bottoms out before then anyway, makes the plots
# easier to read
max_idx = len(xs)//2

mtf10 = empty(height)
mtf20 = empty(height)
mtf30 = empty(height)

for line_idx in range(0,data.shape[0]):
    line = diff(data[line_idx,:])
    fftdata = abs(fft.rfft(line))[0:max_idx]
    fftdata /= max(fftdata)
    if args.yvsmtf:
        mtf10[line_idx] = fftdata[x10lpmm_idx]
        mtf20[line_idx] = fftdata[x20lpmm_idx]
        mtf30[line_idx] = fftdata[x30lpmm_idx]
    else:
        plt.plot(xs[:max_idx], fftdata, linewidth = 1, color='black', alpha = .002)

if args.yvsmtf:
    plt.plot(xs, mtf10, label='10 lp/mm')
    plt.plot(xs, mtf20, label='20 lp/mm')
    plt.plot(xs, mtf30, label='30 lp/mm')
    plt.gca().set_xlabel('y position (mm)')
    plt.gca().set_ylabel('MTF')
    plt.legend(ncol=3)
else:
    plt.gca().set_xlabel('Spatial Freq (lp/mm)')
    plt.gca().set_ylabel('MTF')

if args.title:
    plt.title(args.title)

plt.grid()
plt.savefig(plot_fname, pad_inches=0)