Add cylindrical panorama previews (f=1920 and f=1550)

North is projected onto a cylinder so its right edge curves naturally into
sunrise's left edge, giving a continuous panoramic sweep from North to East.

cyl-f1920-preview.jpg  — f=1920 px (~90° HFOV): most natural perspective
cyl-f1550-preview.jpg  — f=1550 px (~105° HFOV): more pronounced curve

stitch-cameras.py updated to use cylindrical projection.  Focal length is
a tunable argv[4]; refine once the user confirms which curve looks right.
The 36-second inter-shot gap causes cloud drift at the seam — live video
frames captured simultaneously will not have this artefact.

https://claude.ai/code/session_01C4jbd3waXG3eKZYbGUjLUQ

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]})")