From 7a977e9ff1177b58fc67799895b2368226bc9cce Mon Sep 17 00:00:00 2001 From: Nick Downing Date: Wed, 9 Feb 2022 21:09:05 +1100 Subject: [PATCH] Improved try-corner algorithm which computes a separate median distance in each corner, improved normalization algorithm which is insensitive to brightness and contrast (previously the high-pass filter only removed the brightness -- causing many positioning errors), improved diagnostics, some tuning of block sizes --- correlate.py | 283 ++++++++++++++++++++++++++++++--------------------- 1 file changed, 168 insertions(+), 115 deletions(-) diff --git a/correlate.py b/correlate.py index 43bbb6b..9162c67 100755 --- a/correlate.py +++ b/correlate.py @@ -20,39 +20,48 @@ XP = 64 YP = 64 # allowable +/- slippage between pairs -XS = 64 -YS = 64 +XS = 128 +YS = 128 -CORNER_CANDIDATES = 16 +CORNER_CANDIDATES = 10 -CUTOFF0 = 64 -CUTOFF1 = 4 +CUTOFF0 = 8 +CUTOFF1 = 2 EPSILON = 1e-6 -def calc_bandpass(image): +def normalize(image): _, _, cs = image.shape - return numpy.stack( + num = image - numpy.stack( [ - scipy.ndimage.gaussian_filter(image[:, :, i], CUTOFF1, mode = 'mirror') - + #scipy.ndimage.gaussian_filter(image[:, :, i], CUTOFF1, mode = 'mirror') - scipy.ndimage.gaussian_filter(image[:, :, i], CUTOFF0, mode = 'mirror') for i in range(cs) ], 2 ) + denom = numpy.sqrt( + scipy.ndimage.gaussian_filter( + numpy.sum(numpy.square(num), 2), + CUTOFF0, + mode = 'mirror' + ) + ) + denom[denom < EPSILON] = EPSILON + return num / denom[:, :, numpy.newaxis] -def calc_match(bandpass0, bandpass1, xc, yc): +def calc_match(normalized0, normalized1, xc, yc): x0 = xc - XM // 2 y0 = yc - YM // 2 x1 = xc - XM // 2 - XS y1 = yc - YM // 2 - YS - block0 = bandpass0[ + block0 = normalized0[ y0:y0 + YM, x0:x0 + XM, : ] - block1 = bandpass1[ + block1 = normalized1[ y1:y1 + YM + YS * 2, x1:x1 + XM + XS * 2, : @@ -74,13 +83,20 @@ def calc_match(bandpass0, bandpass1, xc, yc): ), 0 ) - #temp = corr - numpy.mean(corr) - #temp /= 10. * numpy.sqrt(numpy.mean(numpy.square(temp))) - #gamma.write_image(f'corr_{i:d}_{j:d}.jpg', temp + .5) - # find slippage from correlation yo, xo = numpy.unravel_index(numpy.argmax(corr), corr.shape) max_corr = corr[yo, xo] + if diag: + temp = .5 + .5 * (corr / max_corr) + temp[ + max(yo - 21, 0):max(yo + 22, 0), + max(xo - 1, 0):max(xo + 2, 0) + ] = 0. + temp[ + max(yo - 1, 0):max(yo + 2, 0), + max(xo - 21, 0):max(xo + 22, 0) + ] = 0. + gamma.write_image(f'diag_{xc:d}_{yc:d}_corr.jpg', temp) if ( max_corr < EPSILON or xo < CUTOFF1 or @@ -91,8 +107,8 @@ def calc_match(bandpass0, bandpass1, xc, yc): return None # estimate position within block of feature being matched - block1 = block1[yo:yo + YM, xo:xo + XM, :] - match = numpy.sum(block0 * block1, 2) + block10 = block1[yo:yo + YM, xo:xo + XM, :] + match = numpy.sum(block0 * block10, 2) #print('xxx', numpy.sum(match), max_corr) #assert False xf = ( @@ -101,35 +117,61 @@ def calc_match(bandpass0, bandpass1, xc, yc): yf = ( numpy.sum(match, 1) @ numpy.arange(YM, dtype = numpy.double) ) / max_corr - #if diag: - # x = int(math.floor(xf)) - # y = int(math.floor(yf)) - - # diag0 = block0 + .5 - # diag0[ - # max(y - 21, 0):max(y + 22, 0), - # max(x - 1, 0):max(x + 2, 0), - # : - # ] = 0. - # diag0[ - # max(y - 1, 0):max(y + 2, 0), - # max(x - 21, 0):max(x + 22, 0), - # : - # ] = 0. - # gamma.write_image(f'diag_{xc:d}_{yc:d}_0.jpg', diag0) - - # diag1 = block1 + .5 - # diag1[ - # max(y - 21, 0):max(y + 22, 0), - # max(x - 1, 0):max(x + 2, 0), - # : - # ] = 0. - # diag1[ - # max(y - 1, 0):max(y + 2, 0), - # max(x - 21, 0):max(x + 22, 0), - # : - # ] = 0. - # gamma.write_image(f'diag_{xc:d}_{yc:d}_1.jpg', diag1) + if diag: + x = int(math.floor(xf)) + y = int(math.floor(yf)) + + temp = match + .5 + temp[ + max(y - 21, 0):max(y + 22, 0), + max(x - 1, 0):max(x + 2, 0) + ] = 0. + temp[ + max(y - 1, 0):max(y + 2, 0), + max(x - 21, 0):max(x + 22, 0) + ] = 0. + gamma.write_image(f'diag_{xc:d}_{yc:d}_match.jpg', temp) + + diag0 = block0 + .5 + diag0[ + max(y - 21, 0):max(y + 22, 0), + max(x - 1, 0):max(x + 2, 0), + : + ] = 0. + diag0[ + max(y - 1, 0):max(y + 2, 0), + max(x - 21, 0):max(x + 22, 0), + : + ] = 0. + gamma.write_image(f'diag_{xc:d}_{yc:d}_block0.jpg', diag0) + + diag1 = block10 + .5 + diag1[ + max(y - 21, 0):max(y + 22, 0), + max(x - 1, 0):max(x + 2, 0), + : + ] = 0. + diag1[ + max(y - 1, 0):max(y + 2, 0), + max(x - 21, 0):max(x + 22, 0), + : + ] = 0. + gamma.write_image(f'diag_{xc:d}_{yc:d}_block10.jpg', diag1) + + x += xo + y += yo + diag1 = block1 + .5 + diag1[ + max(y - 21, 0):max(y + 22, 0), + max(x - 1, 0):max(x + 2, 0), + : + ] = 0. + diag1[ + max(y - 1, 0):max(y + 2, 0), + max(x - 21, 0):max(x + 22, 0), + : + ] = 0. + gamma.write_image(f'diag_{xc:d}_{yc:d}_block11.jpg', diag1) # return offset and feature relative to block centre return xo - XS, yo - YS, xf - XM // 2, yf - YM // 2 @@ -170,14 +212,14 @@ sys.stderr.write(f'write {out_jpg:s}\n') gamma.write_image(out_jpg, image0) ys, xs, cs = shape -xb = (xs // 2 - XM - 2 * XS) // XP -yb = (ys // 2 - YM - 2 * YS) // YP +xb = (xs - XM - 2 * XS) // XP // 2 +yb = (ys - YM - 2 * YS) // YP // 2 print('xb', xb, 'yb', yb) -print('bandpass') -bandpass0 = calc_bandpass(image0) +print('normalize') +normalized0 = normalize(image0) if diag: - gamma.write_image('bandpass0.jpg', bandpass0 + .5) + gamma.write_image('normalized0.jpg', normalized0 + .5) # loop through remaining files comparing each to the previous cumulative_A = numpy.identity(3, numpy.double) @@ -188,15 +230,14 @@ for file in range(1, len(files)): image1 = gamma.read_image(in_jpg) assert image1.shape == shape - print('bandpass') - bandpass1 = calc_bandpass(image1) + print('normalize') + normalized1 = normalize(image1) if diag: - gamma.write_image(f'bandpass{file:d}.jpg', bandpass1 + .5) + gamma.write_image(f'normalized{file:d}.jpg', normalized1 + .5) print('find corner candidates') p_all = [] q_all = [] - corner_candidates = [] for i in range(2): for j in range(2): print('i', i, 'j', j) @@ -212,100 +253,112 @@ for file in range(1, len(files)): xc = XS + XM // 2 + l * XP if j: xc = xs - xc - match = calc_match(bandpass0, bandpass1, xc, yc) + match = calc_match(normalized0, normalized1, xc, yc) if match is not None: offsets.append(match [:2]) xo, yo, xf, yf = match + #xf = 0 + #yf = 0 xf0 = xc + xf yf0 = yc + yf xf1 = xf0 + xo yf1 = yf0 + yo - p_all.append(numpy.array([xf0, yf0], numpy.double)) - q_all.append(numpy.array([xf1, yf1], numpy.double)) blocks.append((xf0, yf0, xf1, yf1)) # find the offset trend (median offset per axis) in (i, j)-corner - k = len(blocks) // 2 + k = len(offsets) // 2 xo_median = sorted([xo for xo, _ in offsets])[k] yo_median = sorted([yo for _, yo in offsets])[k] #print('i', i, 'j', j, 'xo_median', xo_median, 'yo_median', yo_median) - # choose CORNER_CANDIDATES blocks closest to trend in (i, j)-corner + # sort all offsets by their closeness to trend in (i,j)-corner u = numpy.array(offsets, numpy.double) v = numpy.array([xo_median, yo_median], numpy.double) dist = numpy.sum(numpy.square(u - v[numpy.newaxis, :]), 1) - corner_candidates.append( - sorted( - [(dist[i],) + blocks[i] for i in range(len(blocks))] - )[:CORNER_CANDIDATES] - ) - p_all = numpy.stack(p_all, 1) - q_all = numpy.stack(q_all, 1) + + p = [] + q = [] + for _, k in sorted([(dist[k], k) for k in range(dist.shape[0])]): + xf0, yf0, xf1, yf1 = blocks[k] + p.append(numpy.array([xf0, yf0], numpy.double)) + q.append(numpy.array([xf1, yf1], numpy.double)) + p = numpy.stack(p, 1) + q = numpy.stack(q, 1) + + p_all.append(p) + q_all.append(q) # try all combinations of the corner candidates print('try corner candidates') - p = numpy.zeros((2, 4), numpy.double) - q = numpy.zeros((2, 4), numpy.double) best_dist = None best_A = None - best_p = None # for diag - for _, xf00, yf00, xf10, yf10 in corner_candidates[0]: - p[0, 0] = xf00 - p[1, 0] = yf00 - q[0, 0] = xf10 - q[1, 0] = yf10 - for _, xf01, yf01, xf11, yf11 in corner_candidates[1]: - p[0, 1] = xf01 - p[1, 1] = yf01 - q[0, 1] = xf11 - q[1, 1] = yf11 - for _, xf02, yf02, xf12, yf12 in corner_candidates[2]: - p[0, 2] = xf02 - p[1, 2] = yf02 - q[0, 2] = xf12 - q[1, 2] = yf12 - for _, xf03, yf03, xf13, yf13 in corner_candidates[3]: - p[0, 3] = xf03 - p[1, 3] = yf03 - q[0, 3] = xf13 - q[1, 3] = yf13 - - A = perspective.calc_transform(p, q) - dist = list( - numpy.sum( - numpy.square( - q_all - perspective.apply_transform_multi(A, p_all) - ), - 0 + for i in range(min(p_all[0].shape[1], CORNER_CANDIDATES)): + for j in range(min(p_all[1].shape[1], CORNER_CANDIDATES)): + for k in range(min(p_all[2].shape[1], CORNER_CANDIDATES)): + for l in range(min(p_all[3].shape[1], CORNER_CANDIDATES)): + A = perspective.calc_transform( + numpy.stack( + [p_all[0][:, i], p_all[1][:, j], p_all[2][:, k], p_all[3][:, l]], + 1 + ), + numpy.stack( + [q_all[0][:, i], q_all[1][:, j], q_all[2][:, k], q_all[3][:, l]], + 1 ) ) - dist = sorted(dist)[len(dist) // 2] # median position error + dist = 0. + for m in range(4): + dists = list( + numpy.sum( + numpy.square( + q_all[m] - perspective.apply_transform_multi(A, p_all[m]) + ), + 0 + ) + ) + dist += sorted(dists)[len(dists) // 2] # median position error if best_dist is None or dist < best_dist: best_dist = dist best_A = A - best_p = numpy.copy(p) # for diag cumulative_A = best_A @ cumulative_A - print('remap') - out_image1 = perspective.remap_image(cumulative_A, image1) if diag: + diag0 = numpy.copy(image0) + diag1 = numpy.copy(image1) for i in range(4): - xf0, yf0 = numpy.floor(best_p[:, i]).astype(numpy.int32) - - out_image1[ - max(yf0 - 21, 0):max(yf0 + 22, 0), - max(xf0 - 1, 0):max(xf0 + 2, 0), - : - ] = 0. - out_image1[ - max(yf0 - 1, 0):max(yf0 + 2, 0), - max(xf0 - 21, 0):max(xf0 + 22, 0), - : - ] = 0. + for j in range(p_all[i].shape[1]): + xf0, yf0 = numpy.floor(p_all[i][:, j]).astype(numpy.int32) + diag0[ + max(yf0 - 21, 0):max(yf0 + 22, 0), + max(xf0 - 1, 0):max(xf0 + 2, 0), + : + ] = 1. + diag0[ + max(yf0 - 1, 0):max(yf0 + 2, 0), + max(xf0 - 21, 0):max(xf0 + 22, 0), + : + ] = 1. + + xf1, yf1 = numpy.floor(q_all[i][:, j]).astype(numpy.int32) + diag1[ + max(yf1 - 21, 0):max(yf1 + 22, 0), + max(xf1 - 1, 0):max(xf1 + 2, 0), + : + ] = 0. + diag1[ + max(yf1 - 1, 0):max(yf1 + 2, 0), + max(xf1 - 21, 0):max(xf1 + 22, 0), + : + ] = 0. + gamma.write_image(f'diag_file{file:d}_0.jpg', diag0) + gamma.write_image(f'diag_file{file:d}_1.jpg', diag1) + + print('remap') + out_image1 = perspective.remap_image(cumulative_A, image1) sys.stderr.write(f'write {out_jpg:s}\n') gamma.write_image(out_jpg, out_image1) image0 = image1 - bandpass0 = bandpass1 + normalized0 = normalized1 -- 2.34.1