ajeddeloh · December 17, 2020 01:16 · Dec 17, 2020
diff --git a/mkplot.py b/mkplot.py
@@ -0,0 +1,75 @@
+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)