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,
:
),
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
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 = (
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
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)
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)
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
public / stop_motion.git / commitdiff