#!/usr/bin/env python3
"""Stitch north + sunrise — v2.

Lessons from v1: full homography needs lots of good inliers.  We had 5, which
over-fit and produced streaking.  Fixes:
  1. Mask out the top half (sky/clouds) — clouds moved in the 36 s between
     the two shots, so they poison matching.
  2. Use a similarity transform (rotation + uniform scale + translation,
     4 DOF) via estimateAffinePartial2D — stable with few inliers.
  3. Narrow the candidate band to the expected overlap sliver.
"""
import sys, cv2, numpy as np

north_path, sunrise_path, out_prefix = sys.argv[1], sys.argv[2], sys.argv[3]

north   = cv2.imread(north_path)
sunrise = cv2.imread(sunrise_path)
H_n, W_n = north.shape[:2]
H_s, W_s = sunrise.shape[:2]
print(f"north   {W_n}x{H_n}")
print(f"sunrise {W_s}x{H_s}")

# 1) Naive side-by-side baseline so the user always has something to look at.
naive = np.hstack([north, sunrise])
cv2.imwrite(f"{out_prefix}-naive.jpg", naive)

# 2) Take only the ground half (skip sky where clouds change) and the
#    band where the two FOVs should overlap:
#       right 40% of north   ×   left 40% of sunrise
sky_cutoff = 0.45           # keep rows from y = 0.45*H downwards
ts_cutoff  = 0.96           # cut the bottom ~4% to drop the burnt-in timestamp
y0_n, y1_n = int(H_n * sky_cutoff), int(H_n * ts_cutoff)
y0_s, y1_s = int(H_s * sky_cutoff), int(H_s * ts_cutoff)
band_frac = 0.40

n_band = north  [y0_n:y1_n, int(W_n * (1 - band_frac)):]
s_band = sunrise[y0_s:y1_s, : int(W_s * band_frac)]
n_ox, n_oy = int(W_n * (1 - band_frac)), y0_n   # band → full-image offset
s_ox, s_oy = 0, y0_s

orb = cv2.ORB_create(nfeatures=8000, scaleFactor=1.2, nlevels=10, edgeThreshold=15)
kp_n, des_n = orb.detectAndCompute(cv2.cvtColor(n_band, cv2.COLOR_BGR2GRAY), None)
kp_s, des_s = orb.detectAndCompute(cv2.cvtColor(s_band, cv2.COLOR_BGR2GRAY), None)
print(f"ORB keypoints: north-band={len(kp_n)}  sunrise-band={len(kp_s)}")

bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=False)
raw = bf.knnMatch(des_n, des_s, k=2)
good = [m for m, n in raw if m.distance < 0.8 * n.distance]
print(f"good Lowe-ratio matches: {len(good)}")

if len(good) < 6:
    print("Not enough matches to align — writing naive hstack as final.")
    cv2.imwrite(f"{out_prefix}-aligned.jpg", naive)
    sys.exit(0)

src = np.float32([[kp_n[m.queryIdx].pt[0] + n_ox,
                   kp_n[m.queryIdx].pt[1] + n_oy] for m in good])
dst = np.float32([[kp_s[m.trainIdx].pt[0] + s_ox,
                   kp_s[m.trainIdx].pt[1] + s_oy] for m in good])

# Similarity transform: rotation + uniform scale + translation.
M, inl = cv2.estimateAffinePartial2D(src, dst, method=cv2.RANSAC,
                                     ransacReprojThreshold=5.0, maxIters=5000,
                                     confidence=0.99)
n_inl = int(inl.sum()) if inl is not None else 0
print(f"similarity inliers: {n_inl}/{len(good)}")

if M is None or n_inl < 4:
    print("Similarity fit too weak — writing naive hstack as final.")
    cv2.imwrite(f"{out_prefix}-aligned.jpg", naive)
    sys.exit(0)

a, b, tx = M[0]; c, d, ty = M[1]
scale = np.sqrt(a * a + c * c)
rot_deg = np.degrees(np.arctan2(c, a))
print(f"est. transform: scale={scale:.4f}  rot={rot_deg:+.2f}°  tx={tx:+.1f}  ty={ty:+.1f}")

# Warp north into sunrise's coordinate frame, then compose the canvas.
M3 = np.vstack([M, [0, 0, 1]])
n_corners = np.float32([[0, 0], [W_n, 0], [W_n, H_n], [0, H_n]]).reshape(-1, 1, 2)
n_warped  = cv2.perspectiveTransform(n_corners, M3)
s_corners = np.float32([[0, 0], [W_s, 0], [W_s, H_s], [0, H_s]]).reshape(-1, 1, 2)
all_c = np.concatenate([n_warped, s_corners])
[xmin, ymin] = np.int32(all_c.min(axis=0).ravel() - 0.5)
[xmax, ymax] = np.int32(all_c.max(axis=0).ravel() + 0.5)
W, H = xmax - xmin, ymax - ymin
shift = np.array([[1, 0, -xmin], [0, 1, -ymin]], dtype=np.float64)
M_shifted = shift @ M3        # 2x3 after slicing
M_shifted = M_shifted[:2]
print(f"canvas {W}x{H}   shift=({-xmin},{-ymin})")

warped_n = cv2.warpAffine(north, M_shifted, (W, H))

canvas = warped_n.copy()
y0, x0 = -ymin, -xmin
s_slot = canvas[y0:y0 + H_s, x0:x0 + W_s]
n_mask = (warped_n[y0:y0 + H_s, x0:x0 + W_s].sum(axis=2) > 0).astype(np.float32)[..., None]
# Where north has content → 60/40 blend favouring sunrise (base).  Else just sunrise.
blended = (sunrise.astype(np.float32) * (0.4 + 0.6 * (1 - n_mask)) +
           s_slot.astype(np.float32) * (0.6 * n_mask)).astype(np.uint8)
# Simpler: just drop sunrise on top wherever it exists.
canvas[y0:y0 + H_s, x0:x0 + W_s] = sunrise

cv2.imwrite(f"{out_prefix}-aligned.jpg", canvas)
np.savetxt(f"{out_prefix}-similarity.txt", M, fmt="%.8f")
print(f"wrote {out_prefix}-aligned.jpg and {out_prefix}-similarity.txt")

# Also draw a debug image showing matched inliers.
dbg = cv2.drawMatches(
    n_band, kp_n, s_band, kp_s,
    [good[i] for i in range(len(good)) if inl[i]],
    None, flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)
cv2.imwrite(f"{out_prefix}-matches.jpg", dbg)
print(f"wrote {out_prefix}-matches.jpg")