-
Notifications
You must be signed in to change notification settings - Fork 11
/
mtf.py
147 lines (115 loc) · 4.06 KB
/
mtf.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
# From: https://gist.github.com/stefanv/2051954
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage
N = 100
std = 5#3.75
box = 5
# See also: http://www.imatest.com/docs/sharpness/#calc
def findStepEdgeSubpix(x, r=3):
"""Find the position in x that has highest gradient."""
x = abs(ndimage.gaussian_filter(x, 2, order=1))
r = 3
peak = np.argmax(x)
p = np.polyfit(np.r_[peak-r:peak+r+1], np.log(x[peak-r:peak+r+1]), 2)
return -p[1]/(p[0]*2)
def superresEdge(edgeProfiles, n=4, returnBins = False):
"""
Given a bunch of edge profiles, create an average
profile that is n times the size based on the edge positions.
"""
edgePositions = np.array(map(findStepEdgeSubpix, edgeProfiles))
p = np.polyfit(range(len(edgePositions)), edgePositions, 2)
fitPositions = np.polyval(p, range(len(edgePositions)))
meanPos = floor(np.mean(fitPositions))
shifts = (floor(fitPositions) - meanPos)
bins = np.cast[int](np.modf(fitPositions)[0] * n)
edgeProfiles = [ndimage.shift(profile, -round(shft), cval=float('nan'), order=0) for profile, shft in zip(edgeProfiles, shifts)]
toAverage = [[] for i in range(n)]
for bin_i, profile in zip(bins, edgeProfiles):
toAverage[bin_i].append(profile)
result = np.zeros(len(edgeProfiles[0]) * n)
for i in range(n):
if not toAverage[i]: continue
x = np.mean(toAverage[i], 0)
#print 'mean',i, x
result[n-i-1::n] = x # The bins are sorted by increasing order of edge position; we want right-most edge to go in first bin.
# Trim ends -- shifts were marked with nan.
result = result[~isnan(result)]
if returnBins:
return result, edgeProfiles, toAverage
return result, edgeProfiles
def MTF(x, window=True):
"""
Compute the MTF of an edge scan function.
Parameters
----------
x : 1D ndarray
Edge scan function.
window : bool
Whether to apply Hanning windowing to the input.
Notes
-----
The line spread function is the derivative of the edge scan function. The
FFT of the line spread function gives the MTF.
See Also
--------
http://www.cis.rit.edu/research/thesis/bs/2001/perry/thesis.html
"""
y = np.diff(x)
if window:
y = y * np.hanning(len(y))
y = np.append(y, np.zeros(100))
Y = np.fft.fft(y)
return Y[:len(Y) // 2]
if False:
# Generate edge
x = np.zeros(N)
x += 0.1
x[:N // 2] = 1
# Pass through various filters
y0 = x[5:-5]
y1 = ndimage.gaussian_filter1d(x, 3.75)[5:-5]
y2 = np.convolve(x, 1/box * np.ones(box), mode='same')[5:-5]
Y0 = MTF(y0)
Y1 = MTF(y1)
Y2 = MTF(y2)
f, (ax0, ax1) = plt.subplots(1, 2)
ix = np.arange(len(Y1)) / (2 * len(Y1))
ax0.plot(y0, '.-')
ax0.plot(y1, '.-')
ax0.plot(y2, '.-')
ax1.plot(ix, np.abs(Y0), '.-', label='Raw')
ax1.plot(ix, np.abs(Y1), '.-', label='Gaussian %.2f' % std)
ax1.plot(ix, np.abs(Y2), '.-', label='Box %d' % box)
ax1.legend()
plt.show()
if __name__ == '__main__':
x = np.zeros((100,300))
x[:, x.shape[1]//2:] = 1.0
x = ndimage.rotate(x, 5, order=1)
# Mess with the image here:
# x = ndimage.gaussian_filter(x, 1)
x = x[x.shape[0]//2-10:x.shape[0]//2+10, x.shape[1]//2-100:x.shape[1]//2+100]
from pylab import *
figure()
n=4
edges, im, theBins = superresEdge(x, n=n, returnBins=True)
subplot(2,2,1)
imshow(x,interpolation='nearest')
subplot(2,2,3)
imshow(im, interpolation='nearest')
subplot(2,2,2)
plot(edges,'.-')
subplot(2,2,4)
Y = MTF(edges/max(edges))
lpPerPix = n * np.arange(len(Y)) / (2 * len(Y))
plot(lpPerPix, abs(Y),'.-')
if False:
figure()
title('comparing mean of FFT to the ISO 12233 way.')
lpPerPix_ = np.arange(len(MTF(x[0]))) / (2 * len(MTF(x[0])))
plot(lpPerPix_,np.mean([np.abs(MTF(xi)) for xi in x],0) / x.max(),'.-');ylim(0,1.0)
plot(lpPerPix, abs(Y),'g.-')
show()