Make rectangle calculations more intuitive by always referring to their centres

diff --git a/correlate.py b/correlate.py
index fa9abee..1e93a9f 100755 (executable)
--- a/correlate.py
+++ b/correlate.py
@@ -48,15 +48,20 @@ def bandpass(image):
     2
   )
 
-def correlate(image0_bp, x0, y0, image1_bp, x1, y1):
+def correlate(image0_bp, image1_bp, xc, yc):
+  x0 = xc - BLOCKS * BLOCK_SIZE // 2
+  y0 = yc - BLOCKS * BLOCK_SIZE // 2
+  x1 = xc - (BLOCKS + 1) * BLOCK_SIZE // 2
+  y1 = yc - (BLOCKS + 1) * BLOCK_SIZE // 2
+
   block0 = image0_bp[
     y0:y0 + BLOCK_SIZE * BLOCKS,
     x0:x0 + BLOCK_SIZE * BLOCKS,
     :
   ]
   block1 = image1_bp[
-    y1:y1 + BLOCK_SIZE * (BLOCKS + 1),
-    x1:x1 + BLOCK_SIZE * (BLOCKS + 1),
+    y1:y1 + (BLOCKS + 1) * BLOCK_SIZE,
+    x1:x1 + (BLOCKS + 1) * BLOCK_SIZE,
     :
   ]
 
@@ -82,7 +87,7 @@ def correlate(image0_bp, x0, y0, image1_bp, x1, y1):
 
   y, x = numpy.unravel_index(numpy.argmax(corr), corr.shape)
   return (
-    (x, y)
+    (x - BLOCK_SIZE // 2, y - BLOCK_SIZE // 2)
   if (
     x >= CUTOFF1 and
       x <= BLOCK_SIZE - CUTOFF1 and
@@ -103,18 +108,15 @@ gamma.write_image('image1_bp.jpg', image1_bp + .5)
 print('align')
 buckets = [[[], []], [[], []]]
 for i in range(yb - BLOCKS):
-  y1 = i * BLOCK_SIZE
-  y0 = y1 + BLOCK_SIZE // 2
+  yc = i * BLOCK_SIZE + (BLOCKS + 1) * BLOCK_SIZE // 2
   for j in range(xb - BLOCKS):
-    x1 = j * BLOCK_SIZE
-    x0 = x1 + BLOCK_SIZE // 2
-
-    offset = correlate(image0_bp, x0, y0, image1_bp, x1, y1)
+    xc = j * BLOCK_SIZE + (BLOCKS + 1) * BLOCK_SIZE // 2
+    offset = correlate(image0_bp, image1_bp, xc, yc)
     if offset is not None:
       x, y = offset
       print('i', i, 'j', j, 'x', x, 'y', y)
       buckets[i >= (yb - BLOCKS) // 2][j >= (xb - BLOCKS) // 2].append(
-        (x0, y0, x1, y1, x, y)
+        (xc, yc, xc + x, yc + y)
       )
 
 p = []
@@ -122,35 +124,17 @@ q = []
 for i in range(2):
   for j in range(2):
     k = len(buckets[i][j]) // 2
-    xm = sorted([x for _, _, _, _, x, _ in buckets[i][j]])[k]
-    ym = sorted([y for _, _, _, _, _, y in buckets[i][j]])[k]
+    offset = [(xc1 - xc0, yc1 - yc0) for xc0, yc0, xc1, yc1 in buckets[i][j]]
+    xm = sorted([x for x, _ in offset])[k]
+    ym = sorted([y for _, y in offset])[k]
     #print('i', i, 'j', j, 'xm', xm, 'ym', ym)
-    u = numpy.array(
-      [[x, y] for _, _, _, _, x, y in buckets[i][j]],
-      numpy.double
-    )
+    u = numpy.array(offset, numpy.double)
     v = numpy.array([xm, ym], numpy.double)
     k = numpy.argmin(numpy.sum(numpy.square(u - v[numpy.newaxis, :]), 1))
-    x0, y0, x1, y1, x, y = buckets[i][j][k]
+    xc0, yc0, xc1, yc1 = buckets[i][j][k]
     #print('i', i, 'j', j, 'x', x, 'y', y)
-    p.append(
-      numpy.array(
-        [
-          x0 + (BLOCKS * BLOCK_SIZE // 2),
-          y0 + (BLOCKS * BLOCK_SIZE // 2)
-        ],
-        numpy.double
-      )
-    )
-    q.append(
-      numpy.array(
-        [
-          x1 + x + (BLOCKS * BLOCK_SIZE // 2),
-          y1 + y + (BLOCKS * BLOCK_SIZE // 2)
-        ],
-        numpy.double
-      )
-    )
+    p.append(numpy.array([xc0, yc0], numpy.double))
+    q.append(numpy.array([xc1, yc1], numpy.double))
 p = numpy.stack(p, 1)
 q = numpy.stack(q, 1)
 A = perspective.calc_transform(p, q)