Global page assembly with loop closures and lens undistortion

The frame-by-frame stitcher accumulated drift, so later passes landed rotated
and offset, doubling text. myscan/pagegraph.py instead solves all frame poses
together: ECC odometry between consecutive frames (bridging single corrupt
frames), SIFT+RANSAC loop closures between overlapping passes, and a robust
Gauss-Newton pose graph. Loop matches must agree over the whole overlap, which
rejects matches on repeated text. Off-page frames and frames that disagree
with the mosaic are left out of the render.

calib.Undistort applies the EEPROM lens model (OpenCV k1..k3, p1, p2 plus
the tilt homography, as decoded in EEPROM.md); up to 8.8 px at the corners.

scan.py registers frames on a worker thread so capture no longer drops
frames, then assembles the page after the stop press (~11 s for a full
page). stitch.py uses the new assembly by default; --live keeps the old one.

EEPROM.md documents the full calibration blob layout, checksum and signature.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
erik
2026-10-02 13:06:45 +00:00
co-authored by Claude Opus 5.5
parent e70d56d957
commit 026aa6c894
8 changed files with 781 additions and 36 deletions
+27
View File
@@ -13,6 +13,7 @@ white -> colorDestBrightImage = 250).
import os
import struct
import cv2
import numpy as np
EEPROM_BASE = 0x3000
@@ -67,6 +68,8 @@ def parse(data):
"laser_left_m": (d(0x5B), d(0x63)),
"laser_right_m": (d(0x73), d(0x7B)),
"focal_px": d(0x85), "principal_px": (d(0x8D), d(0x95)),
"dist_coeffs": struct.unpack_from("<5d", blob, 0xA5), # OpenCV k1, k2, p1, p2, k3
"tilt_h": struct.unpack_from("<8d", blob, 0xD1) + (1.0,), # 3x3 row-major, pixel space
"pixel_size_m": d(0x111),
"black": samples[0], "white": samples[1], # columns x, y, b, g, r
"gains_rgb": struct.unpack_from("<3i", blob, end + 4), # sensor regs 1, 2, 3
@@ -125,6 +128,30 @@ class FlatField:
return (f - b) * (DEST_WHITE / np.maximum(w - b, 1))
class Undistort:
"""Lens + tilt correction, as the app's generateMapFromParams (EEPROM.md): for each
output pixel p, the source pixel is K * distort(inv(K) * inv(H) * p)."""
def __init__(self, cal, shape=(296, 640)):
h, w = shape
f, (cx, cy) = cal["focal_px"], cal["principal_px"]
K = np.array([[f, 0, cx], [0, f, cy], [0, 0, 1.0]])
M = np.linalg.inv(K) @ np.linalg.inv(np.reshape(cal["tilt_h"], (3, 3)))
k1, k2, p1, p2, k3 = cal["dist_coeffs"]
ys, xs = np.mgrid[0:h, 0:w].astype(np.float64)
X, Y, W = (M @ np.stack([xs.ravel(), ys.ravel(), np.ones(xs.size)])).reshape(3, h, w)
x, y = X / W, Y / W
r2 = x * x + y * y
radial = 1 + k1 * r2 + k2 * r2 ** 2 + k3 * r2 ** 3
xd = x * radial + 2 * p1 * x * y + p2 * (r2 + 2 * x * x)
yd = y * radial + p1 * (r2 + 2 * y * y) + 2 * p2 * x * y
self.map = cv2.convertMaps((f * xd + cx).astype(np.float32), (f * yd + cy).astype(np.float32),
cv2.CV_16SC2)
def __call__(self, frame):
return cv2.remap(frame, *self.map, cv2.INTER_LINEAR, borderMode=cv2.BORDER_REPLICATE)
def load(x=None, cache=CACHE):
if x is None and not (cache and os.path.exists(cache)):
from .xu import XU
+383
View File
@@ -0,0 +1,383 @@
"""Global page assembly: place all frames of a scan at once (pose-graph optimisation).
The live Stitcher places each frame relative to the canvas built so far, so its error
accumulates and a later pass that comes back over earlier content can lock on a wrong
copy. Here every frame keeps its own pose and the poses are solved together:
1. Odometry edges between consecutive frames: laser prediction (with the phase-correlation
override, as in Stitcher), refined by ECC between the two frames when that converges.
After a frame that ECC could not register, the next frame is also registered against
the one before it, which bridges single corrupt frames.
2. Loop-closure edges between non-consecutive keyframes whose current poses overlap:
SIFT matches + RANSAC (rigid), so large drift is recovered, then ECC for sub-pixel.
3. Robust Gauss-Newton over (x, y, theta) per frame. Loop edges that still disagree with
the solution are dropped (repeated content gives false matches) and the solve repeats;
rounds of 2-3 continue with the improved poses so passes far apart find each other.
4. Frames linked to the main page by image-verified edges are blended at their final
pose; the rest (off the page, lifted, corrupt) only had the lasers and are left out,
as are frames that disagree with the mosaic of all the others.
A pose M (3x3) maps frame pixels to canvas pixels; an edge (i, j, T) says M_j ~ M_i @ T.
"""
import os
from concurrent.futures import ThreadPoolExecutor
import numpy as np
import cv2
from scipy.sparse import coo_matrix, diags
from scipy.sparse.linalg import spsolve
from .stitch import (MIN_PHASE_RESP, ROT_GAIN, _aff, _feather, _highpass, compose,
laser_deltas, laser_shift, relative_motion)
ECC_CRIT = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 50, 1e-4)
SIGMA_ECC = 1.0 # px, edge registered by ECC
SIGMA_LASER = 6.0 # px, odometry edge from lasers / phase correlation only
SIGMA_LOOP = 1.5 # px, loop closure
HUBER = 3.0 # robust kernel width, in sigmas
PRUNE = 5.0 # drop loop edges whose residual exceeds this many sigmas
MAX_STEP_ROT = np.radians(2.0) # larger laser rotation between two frames is a miscount
def _rigid(W):
"""2x3 (possibly with a little scale) -> 3x3 rigid transform."""
th = np.arctan2(W[1, 0] - W[0, 1], W[0, 0] + W[1, 1])
return _aff(th, W[0, 2], W[1, 2])
def _ecc(template, image, T):
"""Refine T (image pixels -> template pixels). Returns (T, correlation) or (None, 0)."""
W = np.linalg.inv(T)[:2].astype(np.float32) # ECC warps template coords into image coords
try:
cc, W = cv2.findTransformECC(template, image, W, cv2.MOTION_EUCLIDEAN, ECC_CRIT, None, 5)
except cv2.error:
return None, 0.0
return np.linalg.inv(np.vstack([W, [0, 0, 1]]).astype(np.float64)), cc
def _warp_add(dst, img, M):
"""dst += img warped by M (3x3, img pixels -> dst pixels), touching only the covered box."""
h, w = img.shape[:2]
q = (M @ np.array([[0, 0, 1], [w, 0, 1], [0, h, 1], [w, h, 1.0]]).T)[:2]
x0, y0 = np.maximum(np.floor(q.min(1)).astype(int) - 1, 0)
x1, y1 = np.minimum(np.ceil(q.max(1)).astype(int) + 2, (dst.shape[1], dst.shape[0]))
if x1 > x0 and y1 > y0:
dst[y0:y1, x0:x1] += cv2.warpAffine(img, (_aff(0, -x0, -y0) @ M)[:2], (x1 - x0, y1 - y0),
flags=cv2.INTER_LINEAR)
def _close(A, B, max_deg, max_px):
D = np.linalg.solve(B, A)
return abs(np.degrees(np.arctan2(D[1, 0], D[0, 0]))) <= max_deg and np.hypot(D[0, 2], D[1, 2]) <= max_px
class PageGraph:
def __init__(self, frame_shape, flat=None, undistort=None, rot_gain=ROT_GAIN, min_cc=0.6):
"""flat, undistort: calib.FlatField / calib.Undistort (raw frame -> corrected), or None."""
self.h, self.w = frame_shape[:2]
self.flat, self.undistort, self.rot_gain, self.min_cc = flat, undistort, rot_gain, min_cc
self.center = (self.w / 2, self.h / 2)
self.window = cv2.createHanningWindow((self.w, self.h), cv2.CV_32F)
self.frames, self.hp, self.anchored = [], [], []
self.odo, self.bridges, self.loops = [], [], [] # edges: (i, j, T, sigma)
self.prev_v = None
self.stats = {"frames": 0, "odo_ecc": 0, "odo_laser": 0, "bridges": 0}
@property
def n(self):
return len(self.frames)
def _prep(self, frame):
f = frame.astype(np.float32) if self.flat is None else self.flat(frame)
return f if self.undistort is None else self.undistort(f)
def add(self, frame, v):
"""frame: HxW or HxWx3 uint8 raw frame; v: its record's six raw values."""
f = self._prep(frame)
gray = np.ascontiguousarray(f.mean(axis=2) if f.ndim == 3 else f, np.float32)
hp = _highpass(gray)
v = np.asarray(v)
i = self.n
self.frames.append(frame)
self.hp.append(hp)
self.anchored.append(False)
self.stats["frames"] += 1
if i:
d = laser_deltas(np.stack([self.prev_v, v]))[1]
shift = None
(dx, dy), resp = cv2.phaseCorrelate(self.hp[i - 1], hp, self.window)
if resp > MIN_PHASE_RESP and np.hypot(*(laser_shift(d) - (dx, dy))) > 2:
shift = (dx, dy)
gain = self.rot_gain if abs(self.rot_gain * (d[1] - d[3])) < MAX_STEP_ROT else 0.0
T0 = relative_motion(d, self.center, gain, shift)
T, cc = _ecc(self.hp[i - 1], hp, T0)
if T is not None and cc > self.min_cc and _close(T, T0, 3.0, 15):
self.odo.append((i - 1, i, T, SIGMA_ECC))
self.anchored[i - 1] = self.anchored[i] = True
self.stats["odo_ecc"] += 1
else:
self.odo.append((i - 1, i, T0, SIGMA_LASER))
self.stats["odo_laser"] += 1
if i >= 2 and self.odo[-2][3] != SIGMA_ECC:
Tb0 = self.odo[-2][2] @ self.odo[-1][2]
Tb, cc = _ecc(self.hp[i - 2], hp, Tb0)
if Tb is not None and cc > self.min_cc and _close(Tb, Tb0, 6.0, 60):
self.bridges.append((i - 2, i, Tb, SIGMA_ECC))
self.anchored[i - 2] = self.anchored[i] = True
self.stats["bridges"] += 1
self.prev_v = v
# --- poses -------------------------------------------------------------------
def _chain(self):
P = [np.eye(3)]
for _, _, T, _ in self.odo:
P.append(P[-1] @ T)
return np.array([[M[0, 2], M[1, 2], np.arctan2(M[1, 0], M[0, 0])] for M in P])
def _optimise(self, x, edges, iters=20):
"""Robust Gauss-Newton with frame 0 held fixed. x: (N, 3) -> (x, residual per edge in sigmas)."""
n, E = len(x), len(edges)
a = 150.0 # lever arm of the two residual points
P2 = np.array([[self.center[0] - a, self.center[1]], [self.center[0] + a, self.center[1]]])
I = np.array([e[0] for e in edges])
J = np.array([e[1] for e in edges])
sig = np.array([e[3] for e in edges], np.float64)
Tp = np.einsum("eab,pb->epa", np.array([e[2] for e in edges])[:, :2],
np.hstack([P2, np.ones((2, 1))])) # (E, 2 points, 2): in frame i
Pj = np.broadcast_to(P2, (E, 2, 2))
def rot(c, s, q):
return np.stack([c[:, None] * q[..., 0] - s[:, None] * q[..., 1],
s[:, None] * q[..., 0] + c[:, None] * q[..., 1]], -1)
def drot(c, s, q):
return np.stack([-s[:, None] * q[..., 0] - c[:, None] * q[..., 1],
c[:, None] * q[..., 0] - s[:, None] * q[..., 1]], -1)
pin = np.zeros(n * 3)
pin[:3] = 1e9
rows = np.arange(E * 4).reshape(E, 2, 2)
for _ in range(iters):
ci, si, cj, sj = np.cos(x[I, 2]), np.sin(x[I, 2]), np.cos(x[J, 2]), np.sin(x[J, 2])
r = (rot(cj, sj, Pj) + x[J, None, :2] - rot(ci, si, Tp) - x[I, None, :2]) / sig[:, None, None]
nrm = np.sqrt((r ** 2).sum((1, 2)) / 2)
w = np.where(nrm <= HUBER, 1.0, HUBER / np.maximum(nrm, 1e-9))
dj, di = drot(cj, sj, Pj) / sig[:, None, None], -drot(ci, si, Tp) / sig[:, None, None]
R, C, V = [], [], []
for c in range(2):
for frame, sign, dth in ((J, 1.0, dj), (I, -1.0, di)):
for p in range(2):
R += [rows[:, p, c], rows[:, p, c]]
C += [frame * 3 + c, frame * 3 + 2]
V += [sign / sig, dth[:, p, c]]
Jm = coo_matrix((np.concatenate(V), (np.concatenate(R), np.concatenate(C))),
shape=(E * 4, n * 3)).tocsr()
Wd = diags(np.repeat(w, 4))
H = (Jm.T @ Wd @ Jm + diags(pin + 1e-6)).tocsc()
dx = spsolve(H, -(Jm.T @ (Wd @ r.reshape(-1))))
x = x + dx.reshape(n, 3)
if np.abs(dx).max() < 1e-3:
break
return x, nrm
# --- loop closures ---------------------------------------------------------------
def _keyframes(self, x, step_px=60.0):
keys, last = [], None
for i in range(self.n):
if not self.anchored[i]:
continue
if last is None or np.hypot(*(x[i, :2] - x[last, :2])) > step_px or i - last > 15:
keys.append(i)
last = i
return keys
def _detect(self, i):
g = cv2.normalize(self.hp[i], None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
return cv2.SIFT_create(nfeatures=600).detectAndCompute(g, None)
def _features(self, i):
if i not in self._feat:
self._feat[i] = self._detect(i)
return self._feat[i]
def _match(self, a, b):
"""Rigid transform b -> a from features + ECC, or None."""
ka, da = self._features(a)
kb, db = self._features(b)
if da is None or db is None or len(ka) < 15 or len(kb) < 15:
return None
good = [m for m, n2 in (p for p in cv2.BFMatcher(cv2.NORM_L2).knnMatch(db, da, k=2) if len(p) == 2)
if m.distance < 0.75 * n2.distance]
if len(good) < 12:
return None
pb = np.float32([kb[m.queryIdx].pt for m in good])
pa = np.float32([ka[m.trainIdx].pt for m in good])
A, inl = cv2.estimateAffinePartial2D(pb, pa, method=cv2.RANSAC, ransacReprojThreshold=3.0)
if A is None or inl.sum() < 12 or abs(np.hypot(A[0, 0], A[1, 0]) - 1) > 0.05:
return None
T = _rigid(A)
T2, cc = _ecc(self.hp[a], self.hp[b], T)
if T2 is not None and cc > 0.5 and _close(T2, T, 1.0, 4) and self._agree(a, b, T2):
return T2
return None
def _agree(self, a, b, T, min_ncc=0.85, min_tiles=0.8, tile=32):
"""Whole-overlap check. Repeated text (two lines ending in the same words) gives
matches that ECC confirms on the shared words while the rest of the overlap
disagrees: require high NCC over the overlap and agreement in most textured tiles."""
hb = cv2.warpAffine(self.hp[b], T[:2], (self.w, self.h))
m = cv2.warpAffine(np.ones((self.h, self.w), np.float32), T[:2], (self.w, self.h)) > 0.99
m = cv2.erode(m.astype(np.uint8), np.ones((15, 15), np.uint8)) > 0
if m.sum() < 0.05 * m.size or np.corrcoef(self.hp[a][m], hb[m])[0, 1] < min_ncc:
return False
th, tw = self.h // tile, self.w // tile
cut = lambda im: im[:th * tile, :tw * tile].reshape(th, tile, tw, tile).swapaxes(1, 2).reshape(th, tw, -1) # noqa: E731
A, B, M = cut(self.hp[a]), cut(hb), cut(m.astype(np.float32))
A = A - A.mean(-1, keepdims=True)
B = B - B.mean(-1, keepdims=True)
sa, sb = A.std(-1), B.std(-1)
use = (M.mean(-1) > 0.9) & (np.maximum(sa, sb) > 8)
if not use.any():
return True
ncc = (A * B).mean(-1) / np.maximum(sa * sb, 1e-6)
return (ncc[use] > 0.5).mean() >= min_tiles
def _close_loops(self, x, keys, radius, min_gap=12, tries_per_run=3):
"""Match each keyframe against overlapping later keyframes. Candidates are grouped
into runs of consecutive keyframes (one pass over the spot); each run is tried
nearest first and stops at its first match. Runs are matched in parallel."""
P = x[keys, :2]
jobs = []
for ia, a in enumerate(keys):
d = np.hypot(*(P[ia + 1:] - P[ia]).T)
near = [ib for ib in np.nonzero(d < radius)[0] + ia + 1 if keys[ib] - a >= min_gap]
runs, cur = [], []
for ib in near:
if cur and ib - cur[-1] > 2:
runs.append(cur)
cur = []
cur.append(ib)
runs += [cur] if cur else []
for run in runs:
bs = [keys[ib] for ib in sorted(run, key=lambda k: d[k - ia - 1])[:tries_per_run]]
bs = [b for b in bs if (a, b) not in self._tried]
if bs:
jobs.append((a, bs))
def run_job(job):
a, bs = job
for k, b in enumerate(bs):
T = self._match(a, b)
if T is not None:
return (a, b, T, SIGMA_LOOP), k + 1
return None, len(bs)
new = 0
for (a, bs), (edge, tries) in zip(jobs, self._pool.map(run_job, jobs)):
self._tried.update((a, b) for b in bs[:tries])
self.stats["loop_tries"] += tries
if edge is not None:
self.loops.append(edge)
new += 1
return new
def solve(self, rounds=4):
"""Poses (N, 3) after loop closing and optimisation."""
self._feat, self._tried, self.loops = {}, set(), []
self.stats.update(loop_tries=0, loops=0, pruned=0)
x = self._chain()
keys = self._keyframes(x)
self._pool = ThreadPoolExecutor(os.cpu_count())
self._feat.update(zip(keys, self._pool.map(self._detect, keys)))
base = self.odo + self.bridges
for r in range(rounds):
# Wide search first (drift can be hundreds of px), tighter once poses are good.
new = self._close_loops(x, keys, radius=700 if r == 0 else 500)
for _ in range(5):
x, nrm = self._optimise(x, base + self.loops)
bad = nrm[len(base):] > PRUNE
if not bad.any():
break
self.loops = [e for e, b in zip(self.loops, bad) if not b]
self.stats["pruned"] += int(bad.sum())
if r and not new:
break
self.anchored = self._main_component()
self._verify(x)
self._pool.shutdown()
self._pool, self._feat = None, {}
self.stats["loops"] = len(self.loops)
self.stats["rendered"] = int(sum(self.anchored))
self.residuals = nrm
return x
def _main_component(self):
"""Frames linked to the largest group by image-verified edges. Off-page stretches
register frame to frame but connect to the page only through the lasers."""
parent = list(range(self.n))
def find(k):
while parent[k] != k:
parent[k] = parent[parent[k]]
k = parent[k]
return k
for i, j, _, s in self.odo + self.bridges + self.loops:
if s != SIGMA_LASER:
parent[find(i)] = find(j)
roots = np.array([find(k) for k in range(self.n)])
main = np.bincount(roots).argmax()
return list(roots == main)
def _verify(self, x, min_ncc=0.5, scale=0.5):
"""Drop frames that disagree with the mosaic of all frames (a stretch held only by a
wrong registration). Compared on high-passed images at half resolution."""
use = [i for i in range(self.n) if self.anchored[i]]
S = np.diag([scale, scale, 1.0])
Ms, off, W, H = self._layout(x, use, S)
fe = cv2.resize(_feather(self.h, self.w), None, fx=scale, fy=scale)
small = list(self._pool.map(
lambda i: cv2.resize(self.hp[i], None, fx=scale, fy=scale, interpolation=cv2.INTER_AREA), use))
A = [off @ M @ np.linalg.inv(S) for M in Ms]
acc = np.zeros((H, W), np.float32)
ws = np.zeros((H, W), np.float32)
for im, M in zip(small, A):
_warp_add(acc, im * fe, M)
_warp_add(ws, fe, M)
mos = acc / np.maximum(ws, 1e-6)
h, w = fe.shape
def score(k):
back = cv2.warpAffine(mos, np.linalg.inv(A[k])[:2], (w, h))
return np.corrcoef(small[k][10:-10, 10:-10].ravel(), back[10:-10, 10:-10].ravel())[0, 1]
self.ncc = np.full(self.n, np.nan)
self.ncc[use] = list(self._pool.map(score, range(len(use))))
for i in use:
if not self.ncc[i] >= min_ncc:
self.anchored[i] = False
self.stats["inconsistent"] = int(np.sum(self.ncc < min_ncc))
def _layout(self, x, use, S=np.eye(3)):
"""Poses of frames `use` scaled by S, the offset that puts them all on a canvas, and its size."""
Ms = [S @ _aff(x[i, 2], x[i, 0], x[i, 1]) for i in use]
c = np.array([[0, 0, 1], [self.w, 0, 1], [0, self.h, 1], [self.w, self.h, 1.0]]).T
q = np.concatenate([(M @ c)[:2] for M in Ms], 1)
off = _aff(0, 2 - q[0].min(), 2 - q[1].min())
W, H = int(np.ceil(np.ptp(q[0]))) + 4, int(np.ceil(np.ptp(q[1]))) + 4
return Ms, off, W, H
# --- rendering ---------------------------------------------------------------
def image(self, x=None, white_balance=True, background=255):
x = self.solve() if x is None else x
use = [i for i in range(self.n) if self.anchored[i]]
Ms, off, W, H = self._layout(x, use)
ch = self.frames[0].shape[2] if self.frames[0].ndim == 3 else 0
acc = np.zeros((H, W, ch) if ch else (H, W), np.float32)
wsum = np.zeros((H, W), np.float32)
feather = _feather(self.h, self.w)
for i, M in zip(use, Ms):
f = self._prep(self.frames[i])
_warp_add(acc, f * (feather[..., None] if f.ndim == 3 else feather), off @ M)
_warp_add(wsum, feather, off @ M)
return compose(acc, wsum, white_balance, background)
+17 -12
View File
@@ -177,18 +177,23 @@ class Stitcher:
self.n += 1
def image(self, white_balance=True, background=255):
covered = self.wsum > 0.05
ws = np.maximum(self.wsum, 1e-6)
out = self.acc / (ws[..., None] if self.ch > 1 else ws)
if white_balance:
out = white_patch(out, covered)
m = covered[..., None] if self.ch > 1 else covered
out = np.where(m, out, background).clip(0, 255).astype(np.uint8)
ys, xs = np.nonzero(covered)
if len(ys):
out = out[ys.min():ys.max() + 1, xs.min():xs.max() + 1]
# Frames arrive transposed and mirrored relative to the page.
return cv2.flip(cv2.transpose(out), -1)
return compose(self.acc, self.wsum, white_balance, background)
def compose(acc, wsum, white_balance=True, background=255):
"""Weighted canvas (acc, wsum) -> cropped page image in page orientation."""
covered = wsum > 0.05
ws = np.maximum(wsum, 1e-6)
out = acc / (ws[..., None] if acc.ndim == 3 else ws)
if white_balance:
out = white_patch(out, covered)
m = covered[..., None] if acc.ndim == 3 else covered
out = np.where(m, out, background).clip(0, 255).astype(np.uint8)
ys, xs = np.nonzero(covered)
if len(ys):
out = out[ys.min():ys.max() + 1, xs.min():xs.max() + 1]
# Frames arrive transposed and mirrored relative to the page.
return cv2.flip(cv2.transpose(out), -1)
def white_patch(img, mask, pct=97, white=245):