diff --git a/images/stitch-preview/cyl-f1550-preview.jpg b/images/stitch-preview/cyl-f1550-preview.jpg new file mode 100644 index 0000000..637be37 Binary files /dev/null and b/images/stitch-preview/cyl-f1550-preview.jpg differ diff --git a/images/stitch-preview/cyl-f1920-preview.jpg b/images/stitch-preview/cyl-f1920-preview.jpg new file mode 100644 index 0000000..9c80fba Binary files /dev/null and b/images/stitch-preview/cyl-f1920-preview.jpg differ diff --git a/stitch-cameras.py b/stitch-cameras.py index 55f6701..b9ae9d1 100644 --- a/stitch-cameras.py +++ b/stitch-cameras.py @@ -1,79 +1,108 @@ #!/usr/bin/env python3 """ -Final stitch: north left, sunrise right. -Translation from phase-correlation (tx=-3324, ty=-85). -Feathered gradient blend across the overlap zone. +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. + +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. """ import sys, cv2, numpy as np -north_p, sunrise_p, out_p = sys.argv[1], sys.argv[2], sys.argv[3] -tx_arg = int(sys.argv[4]) if len(sys.argv) > 4 else -3324 -ty_arg = int(sys.argv[5]) if len(sys.argv) > 5 else -85 +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 north = cv2.imread(north_p) sunrise = cv2.imread(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") -tx, ty = tx_arg, ty_arg # north → sunrise translation -print(f"Using tx={tx} ty={ty}") +# ── 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 + 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) -# ── Canvas ──────────────────────────────────────────────────────────────────── -# Sunrise is at the canvas origin (0,0). -# North is at (tx, ty) relative to sunrise — tx is negative so north is LEFT. -canvas_x0 = min(0, tx) # leftmost pixel -canvas_y0 = min(0, ty) -canvas_x1 = max(W, W + tx) -canvas_y1 = max(H, H + ty) -cW = int(canvas_x1 - canvas_x0) -cH = int(canvas_y1 - canvas_y0) +n_cyl = cylindrical_warp(north, f, cx, cy) +s_cyl = cylindrical_warp(sunrise, f, cx, cy) -# Canvas offsets for each image -s_ox, s_oy = int(-canvas_x0), int(-canvas_y0) # sunrise top-left in canvas -n_ox, n_oy = s_ox + tx, s_oy + ty # north top-left in canvas +cv2.imwrite(f"{prefix}-north-cyl.jpg", n_cyl) +cv2.imwrite(f"{prefix}-sunrise-cyl.jpg", s_cyl) +print("Cylindrical projections saved.") + +# ── 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) + +# 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}") + +# 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}") + +# ── 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}") 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 -# Paint north first (background) -nx0, ny0 = int(n_ox), int(n_oy) -canvas[ny0:ny0+H, nx0:nx0+W] = north - -# Paint sunrise on top (base) -sx0, sy0 = int(s_ox), int(s_oy) -canvas[sy0:sy0+H, sx0:sx0+W] = sunrise - -# ── Feathered blend in overlap zone ────────────────────────────────────────── -# Overlap: columns where both images exist in the canvas -ov_x0 = max(sx0, nx0) -ov_x1 = min(sx0 + W, nx0 + W) -ov_y0 = max(sy0, ny0) -ov_y1 = min(sy0 + H, ny0 + H) - +# 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) if ov_x1 > ov_x0 and ov_y1 > ov_y0: ov_w = ov_x1 - ov_x0 - ov_h = ov_y1 - ov_y0 - print(f"Overlap zone: {ov_w}x{ov_h}px at canvas x=[{ov_x0},{ov_x1}]") + 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") - # Horizontal gradient: sunrise fades OUT on the left (where north takes over) - # north fades OUT on the right (where sunrise takes over) - alpha = np.linspace(0.0, 1.0, ov_w, dtype=np.float32) # 0=north, 1=sunrise - alpha = alpha[np.newaxis, :, np.newaxis] # shape (1, ov_w, 1) +# Crop black margin from ty offset +top = max(sy0, ny0) +bot = min(sy0+H, ny0+H) +canvas = canvas[top:bot, :] - n_patch = north [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) - blended = (s_patch * alpha + n_patch * (1.0 - alpha)).astype(np.uint8) - canvas[ov_y0:ov_y1, ov_x0:ov_x1] = blended - -# ── Crop off the black margin at top/bottom from the ty offset ─────────────── -# After the shift there may be a thin black bar top or bottom — crop it. -top_crop = max(sy0, ny0) # first row where BOTH images exist -bottom_crop = min(sy0+H, ny0+H) # last row where either image exists -canvas = canvas[top_crop:bottom_crop, :] - -cv2.imwrite(out_p, canvas, [cv2.IMWRITE_JPEG_QUALITY, 92]) -print(f"Output: {out_p} size={canvas.shape[1]}x{canvas.shape[0]}") - -# ── Save a 50% scaled preview for quick viewing ─────────────────────────────── +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(out_p.replace(".jpg", "-preview.jpg"), preview, - [cv2.IMWRITE_JPEG_QUALITY, 85]) -print(f"Preview saved.") +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]})")