Rewrite stitch-cameras.py: north-only cylindrical warp + colour matching
- Warp only north cylindrically (sunrise completely unmodified) - Focal-length grid search picks the best phase-correlation alignment - Soft validity mask (31x31 erosion) fades out black corners of the warp - Per-channel colour/exposure match: north's treeline strip mean is scaled to match sunrise's, removing the visible sky-colour seam between cameras - Corrected overlap blend: invalid north pixels fall back to sunrise (s_weight = max(alpha, 1-validity)) instead of showing black - Auto-crops left black-corner columns from the final canvas https://claude.ai/code/session_01C4jbd3waXG3eKZYbGUjLUQ
This commit is contained in:
Claude
+144
-79
@@ -1,108 +1,173 @@
|
||||
#!/usr/bin/env python3
|
||||
"""
|
||||
Cylindrical panorama stitch.
|
||||
|
||||
North faces ~0° azimuth, sunrise faces ~90°. Project both onto the same
|
||||
cylinder so the NE overlap region appears at the same y-coordinate in both
|
||||
and the scene is geometrically continuous.
|
||||
Stitch north (cylindrical warp) + sunrise (unmodified).
|
||||
|
||||
Usage:
|
||||
python3 stitch_cyl.py north.jpg sunrise.jpg out_prefix [focal_px]
|
||||
focal_px defaults to 1920 (= ~90° HFOV at 3840 px wide).
|
||||
Try 1300-1920 to span 90°-120° HFOV.
|
||||
python3 stitch-cameras.py north.jpg sunrise.jpg out_prefix [focal_px]
|
||||
|
||||
focal_px: cylindrical focal length in pixels (default: auto-search).
|
||||
Approx guide: 1920 = ~90° HFOV for a 3840px-wide sensor.
|
||||
|
||||
Pipeline
|
||||
--------
|
||||
1. Cylindrical-warp north only; sunrise is never modified.
|
||||
2. Build a soft validity mask for the warped north so black-corner
|
||||
pixels fade out smoothly.
|
||||
3. Phase-correlate edge-enhanced treeline strips to find tx, ty.
|
||||
4. Colour-match north's overlap strip to sunrise's (per-channel gain)
|
||||
so the two cameras' exposure/WB differences don't leave a visible seam.
|
||||
5. Feathered alpha blend in the overlap zone, weighted by north's validity.
|
||||
6. Crop black margins; output full-res panorama + half-size preview.
|
||||
"""
|
||||
import sys, cv2, numpy as np
|
||||
|
||||
north_p = sys.argv[1]
|
||||
sunrise_p= sys.argv[2]
|
||||
prefix = sys.argv[3]
|
||||
f = float(sys.argv[4]) if len(sys.argv) > 4 else 1920.0
|
||||
forced_f = float(sys.argv[4]) if len(sys.argv) > 4 else None
|
||||
|
||||
north = cv2.imread(north_p)
|
||||
sunrise = cv2.imread(sunrise_p)
|
||||
assert north is not None, f"Cannot read {north_p}"
|
||||
assert sunrise is not None, f"Cannot read {sunrise_p}"
|
||||
H, W = north.shape[:2]
|
||||
cx, cy = W / 2.0, H / 2.0
|
||||
print(f"Images: {W}x{H} focal estimate: {f:.0f} px")
|
||||
|
||||
# ── Cylindrical projection ────────────────────────────────────────────────────
|
||||
# Maps output (cylindrical) pixel (xc, yc) ← input (flat) pixel (xs, ys):
|
||||
# theta = xc / f (horizontal angle from optical axis)
|
||||
# xs = f * tan(theta) + cx
|
||||
# ys = yc / cos(theta) + cy
|
||||
def cylindrical_warp(img, f, cx, cy):
|
||||
h, w = img.shape[:2]
|
||||
xc = np.arange(w, dtype=np.float32) - cx
|
||||
yc = np.arange(h, dtype=np.float32) - cy
|
||||
# ── Helpers ───────────────────────────────────────────────────────────────────
|
||||
|
||||
def cylindrical_warp(img, f):
|
||||
xc = np.arange(W, dtype=np.float32) - cx
|
||||
yc = np.arange(H, dtype=np.float32) - cy
|
||||
XC, YC = np.meshgrid(xc, yc)
|
||||
theta = XC / f
|
||||
map_x = (f * np.tan(theta) + cx).astype(np.float32)
|
||||
map_y = (YC / np.cos(theta) + cy).astype(np.float32)
|
||||
return cv2.remap(img, map_x, map_y, cv2.INTER_LINEAR,
|
||||
borderMode=cv2.BORDER_CONSTANT, borderValue=0)
|
||||
theta = XC / f
|
||||
map_x = (f * np.tan(theta) + cx).astype(np.float32)
|
||||
map_y = (YC / np.cos(theta) + cy).astype(np.float32)
|
||||
warped = cv2.remap(img, map_x, map_y, cv2.INTER_LINEAR,
|
||||
borderMode=cv2.BORDER_CONSTANT, borderValue=0)
|
||||
valid = ((map_x >= 0) & (map_x < W) &
|
||||
(map_y >= 0) & (map_y < H)).astype(np.float32)
|
||||
return warped, valid
|
||||
|
||||
n_cyl = cylindrical_warp(north, f, cx, cy)
|
||||
s_cyl = cylindrical_warp(sunrise, f, cx, cy)
|
||||
def phase_corr(n_cyl):
|
||||
"""Return (peak_score, tx, ty) aligning warped-north to sunrise."""
|
||||
bw = int(W * 0.35)
|
||||
y0, y1 = int(H * 0.60), int(H * 0.92)
|
||||
n_strip = cv2.cvtColor(n_cyl [:, W-bw:], cv2.COLOR_BGR2GRAY)[y0:y1].astype(np.float64)
|
||||
s_strip = cv2.cvtColor(sunrise[:, :bw], cv2.COLOR_BGR2GRAY)[y0:y1].astype(np.float64)
|
||||
ne = cv2.Laplacian(n_strip, cv2.CV_64F, ksize=3)
|
||||
se = cv2.Laplacian(s_strip, cv2.CV_64F, ksize=3)
|
||||
N = np.fft.fft2(ne); S = np.fft.fft2(se)
|
||||
R = N * np.conj(S); nrm = np.abs(R); nrm[nrm == 0] = 1
|
||||
corr = np.fft.ifft2(R / nrm).real
|
||||
peak_val = corr.max()
|
||||
dy_r, dx_r = [int(v) for v in np.unravel_index(np.argmax(corr), corr.shape)]
|
||||
rh, rw = corr.shape
|
||||
if dy_r > rh // 2: dy_r -= rh
|
||||
if dx_r > rw // 2: dx_r -= rw
|
||||
return peak_val, -(W - bw + dx_r), -dy_r
|
||||
|
||||
cv2.imwrite(f"{prefix}-north-cyl.jpg", n_cyl)
|
||||
cv2.imwrite(f"{prefix}-sunrise-cyl.jpg", s_cyl)
|
||||
print("Cylindrical projections saved.")
|
||||
def colour_match(src, ref_strip_src, ref_strip_ref):
|
||||
"""Scale each BGR channel of src so strip means match between the two strips."""
|
||||
out = src.astype(np.float32)
|
||||
for c in range(3):
|
||||
s_mean = ref_strip_src[..., c].astype(np.float32).mean()
|
||||
r_mean = ref_strip_ref[..., c].astype(np.float32).mean()
|
||||
gain = np.clip(r_mean / s_mean if s_mean > 1e-3 else 1.0, 0.5, 2.0)
|
||||
out[..., c] = np.clip(out[..., c] * gain, 0, 255)
|
||||
return out
|
||||
|
||||
# ── Phase-correlate on the treeline band of the cylindrical images ────────────
|
||||
y0, y1 = int(H * 0.70), int(H * 0.90)
|
||||
ow = int(W * 0.30) # look in outer 30% of each image
|
||||
n_strip = cv2.cvtColor(n_cyl[:, W-ow:], cv2.COLOR_BGR2GRAY)[y0:y1].astype(np.float64)
|
||||
s_strip = cv2.cvtColor(s_cyl[:, :ow], cv2.COLOR_BGR2GRAY)[y0:y1].astype(np.float64)
|
||||
# ── Focal length search (or forced) ──────────────────────────────────────────
|
||||
focal_candidates = [forced_f] if forced_f else [1400, 1550, 1700, 1800, 1920, 2100, 2300, 2600]
|
||||
print(f"{'focal':>6} {'score':>8} {'tx':>7} {'ty':>5}")
|
||||
best_score, best_f, best_tx, best_ty = -1.0, 1920.0, 0, 0
|
||||
for f in focal_candidates:
|
||||
n_cyl, _ = cylindrical_warp(north, f)
|
||||
score, tx, ty = phase_corr(n_cyl)
|
||||
print(f"{f:6.0f} {score:8.5f} {tx:7d} {ty:5d}")
|
||||
if score > best_score:
|
||||
best_score, best_f, best_tx, best_ty = score, f, tx, ty
|
||||
|
||||
# Normalised cross-power spectrum (phase correlation)
|
||||
N = np.fft.fft2(n_strip)
|
||||
S = np.fft.fft2(s_strip)
|
||||
R = N * np.conj(S)
|
||||
nrm = np.abs(R); nrm[nrm == 0] = 1
|
||||
corr = np.fft.ifft2(R / nrm).real
|
||||
peak = np.unravel_index(np.argmax(corr), corr.shape)
|
||||
dy_r, dx_r = peak
|
||||
rh, rw = corr.shape
|
||||
if dy_r > rh // 2: dy_r -= rh
|
||||
if dx_r > rw // 2: dx_r -= rw
|
||||
print(f"Phase-corr (cylindrical): dx={dx_r} dy={dy_r}")
|
||||
f = best_f
|
||||
tx = best_tx
|
||||
ty = best_ty
|
||||
print(f"\nBest focal={f:.0f} score={best_score:.5f} tx={tx} ty={ty}")
|
||||
|
||||
# Convert from band-strip offsets to full-image translation
|
||||
tx = -(W - ow + dx_r) # how far left north sits relative to sunrise
|
||||
ty = -dy_r
|
||||
print(f"Translation: tx={tx} ty={ty}")
|
||||
# ── Warp + validity mask ──────────────────────────────────────────────────────
|
||||
n_cyl, valid_raw = cylindrical_warp(north, f)
|
||||
|
||||
# ── Build panorama from cylindrical images ───────────────────────────────────
|
||||
xmin = min(0, tx); ymin = min(0, ty)
|
||||
xmax = max(W, W+tx); ymax = max(H, H+ty)
|
||||
cW = int(xmax-xmin); cH = int(ymax-ymin)
|
||||
sx0, sy0 = int(-xmin), int(-ymin) # sunrise origin
|
||||
nx0, ny0 = sx0+tx, sy0+ty # north origin
|
||||
nx0, ny0 = int(nx0), int(ny0)
|
||||
print(f"Canvas: {cW}x{cH}")
|
||||
kernel = np.ones((31, 31), np.float32) / (31 * 31)
|
||||
valid_soft = cv2.filter2D(valid_raw, -1, kernel).clip(0, 1)
|
||||
|
||||
canvas = np.zeros((cH, cW, 3), dtype=np.uint8)
|
||||
# North first (background)
|
||||
canvas[ny0:ny0+H, nx0:nx0+W] = n_cyl
|
||||
# Sunrise on top (base image)
|
||||
canvas[sy0:sy0+H, sx0:sx0+W] = s_cyl
|
||||
# First column of north with >50% valid pixels in the central band
|
||||
central_y0, central_y1 = int(H * 0.10), int(H * 0.90)
|
||||
col_valid = valid_raw[central_y0:central_y1, :].mean(axis=0)
|
||||
first_valid_col = int(np.argmax(col_valid > 0.5))
|
||||
|
||||
# Feathered blend in overlap
|
||||
ov_x0 = max(sx0, nx0); ov_x1 = min(sx0+W, nx0+W)
|
||||
ov_y0 = max(sy0, ny0); ov_y1 = min(sy0+H, ny0+H)
|
||||
# ── Canvas geometry ───────────────────────────────────────────────────────────
|
||||
xmin = min(0, tx); ymin = min(0, ty)
|
||||
xmax = max(W, W + tx); ymax = max(H, H + ty)
|
||||
cW, cH = int(xmax - xmin), int(ymax - ymin)
|
||||
sx0, sy0 = int(-xmin), int(-ymin)
|
||||
nx0, ny0 = int(sx0 + tx), int(sy0 + ty)
|
||||
|
||||
ov_x0 = max(sx0, nx0); ov_x1 = min(sx0 + W, nx0 + W)
|
||||
ov_y0 = max(sy0, ny0); ov_y1 = min(sy0 + H, ny0 + H)
|
||||
|
||||
# ── Colour-match north → sunrise in the overlap strip ────────────────────────
|
||||
# Use a central vertical slice of the overlap (avoid the edge-feather zones).
|
||||
ov_w = ov_x1 - ov_x0
|
||||
# Sample the middle half of the overlap, treeline rows only (avoid sky + timestamp)
|
||||
y_lo, y_hi = int(H * 0.55), int(H * 0.88)
|
||||
mid_n_x0 = ov_x0 - nx0 + ov_w // 4
|
||||
mid_n_x1 = ov_x0 - nx0 + 3 * ov_w // 4
|
||||
mid_s_x0 = ov_x0 - sx0 + ov_w // 4
|
||||
mid_s_x1 = ov_x0 - sx0 + 3 * ov_w // 4
|
||||
n_sample = n_cyl [y_lo:y_hi, mid_n_x0:mid_n_x1]
|
||||
s_sample = sunrise[y_lo:y_hi, mid_s_x0:mid_s_x1]
|
||||
n_cyl_matched = colour_match(n_cyl, n_sample, s_sample)
|
||||
print(f"Colour-matched north to sunrise (overlap centre rows {y_lo}-{y_hi})")
|
||||
|
||||
# ── Composite ────────────────────────────────────────────────────────────────
|
||||
canvas = np.zeros((cH, cW, 3), dtype=np.float32)
|
||||
|
||||
# 1. Paint north (valid content only, black corners excluded)
|
||||
nv = valid_soft[..., np.newaxis]
|
||||
c_s = nx0 + first_valid_col
|
||||
c_e = min(nx0 + W, cW)
|
||||
canvas[ny0:ny0 + H, c_s:c_e] = (
|
||||
n_cyl_matched[:, first_valid_col:c_e - nx0].astype(np.float32)
|
||||
* nv[:, first_valid_col:c_e - nx0]
|
||||
)
|
||||
|
||||
# 2. Sunrise overwrites its entire area (completely unmodified)
|
||||
canvas[sy0:sy0 + H, sx0:sx0 + W] = sunrise.astype(np.float32)
|
||||
|
||||
# 3. Feathered blend in overlap: when north is invalid, fall back to sunrise
|
||||
if ov_x1 > ov_x0 and ov_y1 > ov_y0:
|
||||
ov_w = ov_x1 - ov_x0
|
||||
alpha = np.linspace(0.0, 1.0, ov_w, dtype=np.float32)[np.newaxis, :, np.newaxis]
|
||||
n_patch = n_cyl[ov_y0-ny0:ov_y1-ny0, ov_x0-nx0:ov_x1-nx0].astype(np.float32)
|
||||
s_patch = s_cyl[ov_y0-sy0:ov_y1-sy0, ov_x0-sx0:ov_x1-sx0].astype(np.float32)
|
||||
canvas[ov_y0:ov_y1, ov_x0:ov_x1] = (s_patch*alpha + n_patch*(1-alpha)).astype(np.uint8)
|
||||
print(f"Overlap blend: {ov_w}x{ov_y1-ov_y0}px")
|
||||
alpha = np.linspace(0.0, 1.0, ov_w, dtype=np.float32)[np.newaxis, :, np.newaxis]
|
||||
n_patch = n_cyl_matched[ov_y0 - ny0:ov_y1 - ny0, ov_x0 - nx0:ov_x1 - nx0].astype(np.float32)
|
||||
s_patch = sunrise [ov_y0 - sy0:ov_y1 - sy0, ov_x0 - sx0:ov_x1 - sx0].astype(np.float32)
|
||||
nv_patch = valid_soft [ov_y0 - ny0:ov_y1 - ny0, ov_x0 - nx0:ov_x1 - nx0][..., np.newaxis]
|
||||
s_weight = np.maximum(alpha, 1.0 - nv_patch)
|
||||
canvas[ov_y0:ov_y1, ov_x0:ov_x1] = s_patch * s_weight + n_patch * (1.0 - s_weight)
|
||||
print(f"Overlap blend: {ov_w}x{ov_y1 - ov_y0}px")
|
||||
|
||||
# Crop black margin from ty offset
|
||||
top = max(sy0, ny0)
|
||||
bot = min(sy0+H, ny0+H)
|
||||
canvas = canvas[top:bot, :]
|
||||
out = canvas.clip(0, 255).astype(np.uint8)
|
||||
|
||||
cv2.imwrite(f"{prefix}-panorama.jpg", canvas, [cv2.IMWRITE_JPEG_QUALITY, 92])
|
||||
preview = cv2.resize(canvas, (canvas.shape[1]//2, canvas.shape[0]//2))
|
||||
cv2.imwrite(f"{prefix}-preview.jpg", preview, [cv2.IMWRITE_JPEG_QUALITY, 88])
|
||||
print(f"Panorama → {prefix}-panorama.jpg ({canvas.shape[1]}x{canvas.shape[0]})")
|
||||
# Vertical crop to rows where both images exist
|
||||
top = max(sy0, ny0); bot = min(sy0 + H, ny0 + H)
|
||||
out = out[top:bot, :]
|
||||
|
||||
# Left-crop black corner columns
|
||||
col_has_content = out.max(axis=(0, 2)) > 0
|
||||
left_crop = int(np.argmax(col_has_content))
|
||||
out = out[:, left_crop:]
|
||||
print(f"Final size: {out.shape[1]}x{out.shape[0]}")
|
||||
|
||||
cv2.imwrite(f"{prefix}-panorama.jpg", out, [cv2.IMWRITE_JPEG_QUALITY, 93])
|
||||
half = cv2.resize(out, (out.shape[1] // 2, out.shape[0] // 2))
|
||||
cv2.imwrite(f"{prefix}-preview.jpg", half, [cv2.IMWRITE_JPEG_QUALITY, 88])
|
||||
|
||||
with open(f"{prefix}-params.txt", "w") as fh:
|
||||
fh.write(f"focal={f:.0f}\ntx={tx}\nty={ty}\n")
|
||||
|
||||
Reference in New Issue
Block a user