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