13 # size of block that will be matched
14 # (correlate a block this size against a block with added slippage all around)
18 # pitch between the block centres
22 # allowable +/- slippage between pairs
26 CORNER_CANDIDATES = 16
33 def calc_bandpass(image):
34 _, _, cs = image.shape
37 scipy.ndimage.gaussian_filter(image[:, :, i], CUTOFF1, mode = 'mirror') -
38 scipy.ndimage.gaussian_filter(image[:, :, i], CUTOFF0, mode = 'mirror')
44 def calc_match(bandpass0, bandpass1, xc, yc):
47 x1 = xc - XM // 2 - XS
48 y1 = yc - YM // 2 - YS
61 # note: swapping block1, block0 flips the output (subtracts x and y
62 # from BLOCK_SIZE) and we need this to find matching part of image1
66 scipy.signal.correlate(
77 #temp = corr - numpy.mean(corr)
78 #temp /= 10. * numpy.sqrt(numpy.mean(numpy.square(temp)))
79 #gamma.write_image(f'corr_{i:d}_{j:d}.jpg', temp + .5)
81 # find slippage from correlation
82 yo, xo = numpy.unravel_index(numpy.argmax(corr), corr.shape)
83 max_corr = corr[yo, xo]
87 xo > XS * 2 - CUTOFF1 or
93 # estimate position within block of feature being matched
94 block1 = block1[yo:yo + YM, xo:xo + XM, :]
95 match = numpy.sum(block0 * block1, 2)
96 #print('xxx', numpy.sum(match), max_corr)
99 numpy.sum(match, 0) @ numpy.arange(XM, dtype = numpy.double)
102 numpy.sum(match, 1) @ numpy.arange(YM, dtype = numpy.double)
105 # x = int(math.floor(xf))
106 # y = int(math.floor(yf))
108 # diag0 = block0 + .5
110 # max(y - 21, 0):max(y + 22, 0),
111 # max(x - 1, 0):max(x + 2, 0),
115 # max(y - 1, 0):max(y + 2, 0),
116 # max(x - 21, 0):max(x + 22, 0),
119 # gamma.write_image(f'diag_{xc:d}_{yc:d}_0.jpg', diag0)
121 # diag1 = block1 + .5
123 # max(y - 21, 0):max(y + 22, 0),
124 # max(x - 1, 0):max(x + 2, 0),
128 # max(y - 1, 0):max(y + 2, 0),
129 # max(x - 21, 0):max(x + 22, 0),
132 # gamma.write_image(f'diag_{xc:d}_{yc:d}_1.jpg', diag1)
134 # return offset and feature relative to block centre
135 return xo - XS, yo - YS, xf - XM // 2, yf - YM // 2
138 if len(sys.argv) >= 2 and sys.argv[1] == '--diag':
141 if len(sys.argv) < 3:
142 print(f'usage: {sys.argv[0]} in_prefix out_prefix')
144 in_prefix = sys.argv[1]
145 out_prefix = sys.argv[2]
147 # compile list of input and output files by directory search
148 i = in_prefix.rfind('/') + 1
149 in_prefix0 = in_prefix[:i] # directory part
150 in_prefix1 = in_prefix[i:] # filename part
152 (file, out_prefix + file[len(in_prefix):])
156 for file in os.listdir(in_prefix0 if len(in_prefix0) else '.')
157 if file[:len(in_prefix1)] == in_prefix1
162 # first file is special as no transformation needs to be done
163 in_jpg, out_jpg = files[0]
165 print(f'read {in_jpg:s}')
166 image0 = gamma.read_image(in_jpg)
169 sys.stderr.write(f'write {out_jpg:s}\n')
170 gamma.write_image(out_jpg, image0)
173 xb = (xs // 2 - XM - 2 * XS) // XP
174 yb = (ys // 2 - YM - 2 * YS) // YP
175 print('xb', xb, 'yb', yb)
178 bandpass0 = calc_bandpass(image0)
180 gamma.write_image('bandpass0.jpg', bandpass0 + .5)
182 # loop through remaining files comparing each to the previous
183 cumulative_A = numpy.identity(3, numpy.double)
184 for file in range(1, len(files)):
185 in_jpg, out_jpg = files[file]
187 print(f'read {in_jpg:s}')
188 image1 = gamma.read_image(in_jpg)
189 assert image1.shape == shape
192 bandpass1 = calc_bandpass(image1)
194 gamma.write_image(f'bandpass{file:d}.jpg', bandpass1 + .5)
196 print('find corner candidates')
199 corner_candidates = []
202 print('i', i, 'j', j)
204 # correlate blocks in (i, j)-corner
208 yc = YS + YM // 2 + k * YP
212 xc = XS + XM // 2 + l * XP
215 match = calc_match(bandpass0, bandpass1, xc, yc)
216 if match is not None:
217 offsets.append(match [:2])
218 xo, yo, xf, yf = match
223 p_all.append(numpy.array([xf0, yf0], numpy.double))
224 q_all.append(numpy.array([xf1, yf1], numpy.double))
225 blocks.append((xf0, yf0, xf1, yf1))
227 # find the offset trend (median offset per axis) in (i, j)-corner
229 xo_median = sorted([xo for xo, _ in offsets])[k]
230 yo_median = sorted([yo for _, yo in offsets])[k]
231 #print('i', i, 'j', j, 'xo_median', xo_median, 'yo_median', yo_median)
233 # choose CORNER_CANDIDATES blocks closest to trend in (i, j)-corner
234 u = numpy.array(offsets, numpy.double)
235 v = numpy.array([xo_median, yo_median], numpy.double)
236 dist = numpy.sum(numpy.square(u - v[numpy.newaxis, :]), 1)
237 corner_candidates.append(
239 [(dist[i],) + blocks[i] for i in range(len(blocks))]
240 )[:CORNER_CANDIDATES]
242 p_all = numpy.stack(p_all, 1)
243 q_all = numpy.stack(q_all, 1)
245 # try all combinations of the corner candidates
246 print('try corner candidates')
247 p = numpy.zeros((2, 4), numpy.double)
248 q = numpy.zeros((2, 4), numpy.double)
251 best_p = None # for diag
252 for _, xf00, yf00, xf10, yf10 in corner_candidates[0]:
257 for _, xf01, yf01, xf11, yf11 in corner_candidates[1]:
262 for _, xf02, yf02, xf12, yf12 in corner_candidates[2]:
267 for _, xf03, yf03, xf13, yf13 in corner_candidates[3]:
273 A = perspective.calc_transform(p, q)
277 q_all - perspective.apply_transform_multi(A, p_all)
282 dist = sorted(dist)[len(dist) // 2] # median position error
284 if best_dist is None or dist < best_dist:
287 best_p = numpy.copy(p) # for diag
288 cumulative_A = best_A @ cumulative_A
291 out_image1 = perspective.remap_image(cumulative_A, image1)
294 xf0, yf0 = numpy.floor(best_p[:, i]).astype(numpy.int32)
297 max(yf0 - 21, 0):max(yf0 + 22, 0),
298 max(xf0 - 1, 0):max(xf0 + 2, 0),
302 max(yf0 - 1, 0):max(yf0 + 2, 0),
303 max(xf0 - 21, 0):max(xf0 + 22, 0),
307 sys.stderr.write(f'write {out_jpg:s}\n')
308 gamma.write_image(out_jpg, out_image1)
311 bandpass0 = bandpass1