Skip to content

Instantly share code, notes, and snippets.

Created December 17, 2020 01:16
Show Gist options
  • Save ajeddeloh/07073bae201f0024e44295d8ecdb6600 to your computer and use it in GitHub Desktop.
Save ajeddeloh/07073bae201f0024e44295d8ecdb6600 to your computer and use it in GitHub Desktop.
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
plt.figure(figsize=(8,6), dpi=100, constrained_layout=True)
parser = ArgumentParser(description = "Generate MTF curves from images")
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]
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_xlabel('Spatial Freq (lp/mm)')
if args.title:
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