Skip to content

Instantly share code, notes, and snippets.

@ajeddeloh

ajeddeloh/mkplot.py

Created Dec 17, 2020
Embed
What would you like to do?
Generate MTF curves from razorblade images
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment