mirror of
https://github.com/storytold/lightcraft.git
synced 2026-10-26 11:09:23 -06:00
LR-EDIT-OPTICS-PROFILE: script to measure Sony distortion against camera JPEGs
Renders each ARW with the embedded correction on and off and compares it with the camera's own JPEG: residual offsets by radius after a similarity fit, checkerboard line straightness, and the ratio of the camera's correction to LightCraft's. Checked on synthetic boards with known radial distortion (corrected: 0.04 px median, ratio 1.000; 80 % corrected: ratio 1.26, expected about 1.25-1.28). Needs numpy, scipy and opencv-python-headless. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015gRtdWF1E9358FfzhzYqCS
This commit is contained in:
Claude
Executable
+361
@@ -0,0 +1,361 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Compare LightCraft's Sony embedded distortion correction with the camera's own JPEG.
|
||||
|
||||
sony_distortion_check.py --cli target/release/lightcraft-cli --work WORKDIR PHOTO_DIR [PHOTO_DIR…]
|
||||
|
||||
Every ARW with a JPG of the same name is rendered by LightCraft twice, embedded correction on
|
||||
and off, and measured against the camera JPEG (shot with Lens Comp. Distortion on Auto):
|
||||
|
||||
* agreement: offsets of matching features after the best similarity (scale, rotation, shift)
|
||||
between the render and the JPEG, by distance from the centre (0 = centre, 1 = corner);
|
||||
* straightness: how far each checkerboard row and column bows away from a straight line.
|
||||
|
||||
Checkerboard shots use the board's corners (sub-pixel), other shots SIFT features (agreement
|
||||
only). Offsets are in pixels of the camera JPEG (6000 px wide on the A7 III); the "@2376"
|
||||
columns rescale them to the 2376 px exports docs/sony-lens-corrections.md reports.
|
||||
Needs numpy, scipy and opencv-python-headless.
|
||||
"""
|
||||
|
||||
import argparse
|
||||
import json
|
||||
import os
|
||||
import subprocess
|
||||
import sys
|
||||
|
||||
import cv2
|
||||
import numpy as np
|
||||
from scipy.spatial import cKDTree
|
||||
|
||||
WORK_LONG_EDGE = 3000 # detection resolution; coordinates are rescaled to the JPEG's pixels
|
||||
BINS = [0.0, 0.25, 0.5, 0.75, 0.9, 1.01]
|
||||
REPORT_WIDTH = 2376
|
||||
|
||||
|
||||
def render(cli, arw, out, on):
|
||||
settings = out + ".settings.json"
|
||||
with open(settings, "w") as f:
|
||||
json.dump({"optics": {"lens_profile": on, "profile_distortion": 100}}, f)
|
||||
if not os.path.exists(out):
|
||||
subprocess.run([cli, "render", arw, "-o", out, "--settings", settings], check=True, capture_output=True)
|
||||
return out
|
||||
|
||||
|
||||
def load_gray(path):
|
||||
"""Grey image at the working resolution and the factor back to its own full-size pixels."""
|
||||
img = cv2.imread(path, cv2.IMREAD_GRAYSCALE)
|
||||
if img is None:
|
||||
raise RuntimeError(f"cannot read {path}")
|
||||
h, w = img.shape
|
||||
k = WORK_LONG_EDGE / max(w, h)
|
||||
small = cv2.resize(img, (round(w * k), round(h * k)), interpolation=cv2.INTER_AREA)
|
||||
return small, 1.0 / k, (w, h)
|
||||
|
||||
|
||||
# ---------------------------------------------------------------- checkerboard corners
|
||||
|
||||
def detect_corners(gray):
|
||||
"""X-corners of a checkerboard: Hessian saddle points, refined with cornerSubPix."""
|
||||
g = cv2.GaussianBlur(gray.astype(np.float32), (0, 0), 3.0)
|
||||
gxx = cv2.Sobel(g, cv2.CV_32F, 2, 0, ksize=5)
|
||||
gyy = cv2.Sobel(g, cv2.CV_32F, 0, 2, ksize=5)
|
||||
gxy = cv2.Sobel(g, cv2.CV_32F, 1, 1, ksize=5)
|
||||
saddle = gxy * gxy - gxx * gyy
|
||||
top = np.percentile(saddle, 99.95)
|
||||
if top <= 0:
|
||||
return np.zeros((0, 2), np.float32)
|
||||
peaks = (saddle == cv2.dilate(saddle, np.ones((15, 15), np.uint8))) & (saddle > 0.08 * top)
|
||||
ys, xs = np.nonzero(peaks)
|
||||
pts = np.stack([xs, ys], 1).astype(np.float32).reshape(-1, 1, 2)
|
||||
if len(pts) == 0:
|
||||
return pts.reshape(-1, 2)
|
||||
crit = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 50, 0.005)
|
||||
pts = cv2.cornerSubPix(gray, pts, (11, 11), (-1, -1), crit)
|
||||
return pts.reshape(-1, 2).astype(np.float64)
|
||||
|
||||
|
||||
def label_lattice(pts, centre):
|
||||
"""Integer (i, j) grid labels for corners connected to the one nearest the centre.
|
||||
i grows along +x, j along +y; local lattice vectors follow the distortion gradually."""
|
||||
if len(pts) < 20:
|
||||
return {}
|
||||
tree = cKDTree(pts)
|
||||
start = int(tree.query(centre)[1])
|
||||
d, idx = tree.query(pts[start], k=5)
|
||||
offsets = pts[idx[1:]] - pts[start]
|
||||
spacing = float(np.median(d[1:]))
|
||||
a = offsets[np.argmax(offsets @ np.array([1.0, 0.0]))]
|
||||
b = offsets[np.argmax(offsets @ np.array([0.0, 1.0]))]
|
||||
if abs(np.dot(a, b)) > 0.5 * np.linalg.norm(a) * np.linalg.norm(b):
|
||||
return {}
|
||||
labels = {start: (0, 0)}
|
||||
taken = {(0, 0): start}
|
||||
vec = {start: (a, b)}
|
||||
queue = [start]
|
||||
while queue:
|
||||
p = queue.pop()
|
||||
i, j = labels[p]
|
||||
va, vb = vec[p]
|
||||
for di, dj, step in ((1, 0, va), (-1, 0, -va), (0, 1, vb), (0, -1, -vb)):
|
||||
key = (i + di, j + dj)
|
||||
if key in taken:
|
||||
continue
|
||||
dist, q = tree.query(pts[p] + step)
|
||||
if dist > 0.3 * np.linalg.norm(step) or q in labels:
|
||||
continue
|
||||
found = pts[q] - pts[p]
|
||||
na = found * di if di else va
|
||||
nb = found * dj if dj else vb
|
||||
labels[q] = key
|
||||
taken[key] = q
|
||||
vec[q] = (na, nb)
|
||||
queue.append(q)
|
||||
return labels if len(labels) > 50 else {}
|
||||
|
||||
|
||||
def board(path):
|
||||
gray, back, size = load_gray(path)
|
||||
pts = detect_corners(gray)
|
||||
labels = label_lattice(pts, np.array(gray.shape[::-1]) / 2.0)
|
||||
full = {lab: pts[k] * back for k, lab in labels.items()}
|
||||
return full, size
|
||||
|
||||
|
||||
def straightness(grid):
|
||||
"""Largest deviation from a fitted straight line over every row/column with 8+ corners."""
|
||||
lines = {}
|
||||
for (i, j), p in grid.items():
|
||||
lines.setdefault(("col", i), []).append(p)
|
||||
lines.setdefault(("row", j), []).append(p)
|
||||
worst = []
|
||||
for pts in lines.values():
|
||||
if len(pts) < 8:
|
||||
continue
|
||||
pts = np.array(pts)
|
||||
keep = np.ones(len(pts), bool)
|
||||
for _ in range(3): # drop stray corners, then refit
|
||||
c = pts[keep].mean(0)
|
||||
_, _, vt = np.linalg.svd(pts[keep] - c)
|
||||
dev = (pts - c) @ vt[1]
|
||||
keep = np.abs(dev) <= max(1.0, 4 * np.median(np.abs(dev[keep])))
|
||||
if keep.sum() < 8:
|
||||
break
|
||||
if keep.sum() >= 8:
|
||||
worst.append(np.abs(dev[keep]).max())
|
||||
if not worst:
|
||||
return None
|
||||
worst = np.array(worst)
|
||||
return {"lines": len(worst), "median": float(np.median(worst)), "max": float(worst.max())}
|
||||
|
||||
|
||||
def match_boards(src, src_size, dst, dst_size):
|
||||
"""Corner pairs between two labellings of the same shot. A board is periodic, so a fit can't
|
||||
tell a one-square shift from the right match; the corners near the centre can, unfitted:
|
||||
the two framings differ by a few percent of zoom, far less than a square there."""
|
||||
sc, dc = np.array(src_size) / 2.0, np.array(dst_size) / 2.0
|
||||
half_diag = np.hypot(*dst_size) / 2.0
|
||||
best = None
|
||||
for oi in range(-3, 4):
|
||||
for oj in range(-3, 4):
|
||||
pairs = [(p, dst[(i + oi, j + oj)]) for (i, j), p in src.items() if (i + oi, j + oj) in dst]
|
||||
central = [np.linalg.norm((p - sc) - (q - dc)) for p, q in pairs if np.linalg.norm(q - dc) < 0.3 * half_diag]
|
||||
if len(pairs) < 50 or len(central) < 10:
|
||||
continue
|
||||
err = float(np.median(central))
|
||||
if best is None or err < best[0]:
|
||||
best = (err, pairs)
|
||||
if best is None:
|
||||
return None, None
|
||||
s, d = map(np.array, zip(*best[1]))
|
||||
return s, d
|
||||
|
||||
|
||||
# ---------------------------------------------------------------- SIFT (other scenes)
|
||||
|
||||
def sift_pairs(src_path, dst_path):
|
||||
"""SIFT matches; both point sets in the destination's full-size pixels."""
|
||||
a, ka, size_a = load_gray(src_path)
|
||||
b, kb, size_b = load_gray(dst_path)
|
||||
sift = cv2.SIFT_create(nfeatures=20000)
|
||||
pa, da = sift.detectAndCompute(a, None)
|
||||
pb, db = sift.detectAndCompute(b, None)
|
||||
if da is None or db is None:
|
||||
return None, None
|
||||
good = [m for m, n in cv2.BFMatcher().knnMatch(da, db, k=2) if m.distance < 0.75 * n.distance]
|
||||
if len(good) < 30:
|
||||
return None, None
|
||||
s = np.array([pa[m.queryIdx].pt for m in good]) * ka * (size_b[0] / size_a[0])
|
||||
d = np.array([pb[m.trainIdx].pt for m in good]) * kb
|
||||
return s, d
|
||||
|
||||
|
||||
# ---------------------------------------------------------------- geometry
|
||||
|
||||
def fit_similarity(s, d):
|
||||
zs = s[:, 0] + 1j * s[:, 1]
|
||||
zd = d[:, 0] + 1j * d[:, 1]
|
||||
ms, md = zs.mean(), zd.mean()
|
||||
c = np.vdot(zs - ms, zd - md) / np.vdot(zs - ms, zs - ms)
|
||||
return c, md - c * ms
|
||||
|
||||
|
||||
def apply(model, p):
|
||||
c, t = model
|
||||
z = c * (p[:, 0] + 1j * p[:, 1]) + t
|
||||
return np.stack([z.real, z.imag], 1)
|
||||
|
||||
|
||||
def agreement(s, d, size):
|
||||
"""Similarity-fit render points s onto JPEG points d; residual statistics by radius."""
|
||||
pairs = {tuple(q): p for q, p in zip(d, s)}
|
||||
keep = np.ones(len(s), bool)
|
||||
# A loose first pass, then reject gross mismatches (SIFT outliers, mislabelled corners).
|
||||
for limit in (0.02 * size[0], None, None):
|
||||
model = fit_similarity(s[keep], d[keep])
|
||||
res = np.linalg.norm(apply(model, s) - d, axis=1)
|
||||
cut = limit if limit else max(5 * np.median(res[keep]), 2.0)
|
||||
keep = res < cut
|
||||
if keep.sum() < 20:
|
||||
return None
|
||||
model = fit_similarity(s[keep], d[keep])
|
||||
s, d = s[keep], d[keep]
|
||||
fitted = apply(model, s)
|
||||
res = fitted - d
|
||||
centre = np.array(size) / 2.0
|
||||
half_diag = np.hypot(*size) / 2.0
|
||||
rel = d - centre
|
||||
r = np.linalg.norm(rel, axis=1) / half_diag
|
||||
radial = np.sum(res * rel, axis=1) / np.maximum(np.linalg.norm(rel, axis=1), 1e-9)
|
||||
mag = np.linalg.norm(res, axis=1)
|
||||
bins = []
|
||||
for lo, hi in zip(BINS, BINS[1:]):
|
||||
m = (r >= lo) & (r < hi)
|
||||
bins.append({"r": f"{lo:.2f}-{min(hi, 1):.2f}", "n": int(m.sum()),
|
||||
"median": float(np.median(mag[m])) if m.any() else None,
|
||||
"radial": float(np.mean(radial[m])) if m.any() else None})
|
||||
return {"n": int(len(s)), "scale": float(abs(model[0])), "median": float(np.median(mag)),
|
||||
"p95": float(np.percentile(mag, 95)), "max": float(mag.max()), "max_r": float(r.max()),
|
||||
"bins": bins, "pairs": pairs}
|
||||
|
||||
|
||||
def correction_fraction(on, off, size):
|
||||
"""How the camera's correction compares with LightCraft's: JPEG ≈ off + α·(on − off) radially,
|
||||
both renders aligned to the JPEG on the same corners. 1 = the same amount as the camera;
|
||||
above 1, the camera corrects more than LightCraft."""
|
||||
keys = [k for k in on["pairs"] if k in off["pairs"]]
|
||||
if len(keys) < 50:
|
||||
return None
|
||||
d = np.array(keys)
|
||||
raw_on = np.array([on["pairs"][k] for k in keys])
|
||||
raw_off = np.array([off["pairs"][k] for k in keys])
|
||||
keep = np.ones(len(d), bool)
|
||||
for _ in range(4):
|
||||
a = apply(fit_similarity(raw_on[keep], d[keep]), raw_on)
|
||||
b = apply(fit_similarity(raw_off[keep], d[keep]), raw_off)
|
||||
rel = d - np.array(size) / 2.0
|
||||
unit = rel / np.maximum(np.linalg.norm(rel, axis=1, keepdims=True), 1e-9)
|
||||
u = np.sum((a - b) * unit, axis=1)
|
||||
v = np.sum((d - b) * unit, axis=1)
|
||||
if np.median(np.abs(u[keep])) < 1.0: # too little distortion at this focal length to tell
|
||||
return None
|
||||
alpha = float(u[keep] @ v[keep] / (u[keep] @ u[keep]))
|
||||
e = np.abs(v - alpha * u)
|
||||
keep = e <= max(0.5, 4 * np.median(e[keep])) # mismatched corners, not the lens
|
||||
return alpha
|
||||
|
||||
|
||||
# ---------------------------------------------------------------- driver
|
||||
|
||||
def jpeg_for(arw):
|
||||
stem = os.path.splitext(arw)[0]
|
||||
for ext in (".JPG", ".jpg", ".JPEG", ".jpeg"):
|
||||
if os.path.exists(stem + ext):
|
||||
return stem + ext
|
||||
return None
|
||||
|
||||
|
||||
def exif_summary(jpg):
|
||||
from PIL import Image
|
||||
tags = Image.open(jpg).getexif().get_ifd(0x8769)
|
||||
lens = tags.get(0xA434, "?")
|
||||
focal = tags.get(0x920A)
|
||||
fnum = tags.get(0x829D)
|
||||
return f"{lens} {float(focal):.0f}mm f/{float(fnum):.1f}" if focal and fnum else str(lens)
|
||||
|
||||
|
||||
def measure(cli, work, arw):
|
||||
jpg = jpeg_for(arw)
|
||||
if not jpg:
|
||||
return None
|
||||
name = os.path.splitext(os.path.basename(arw))[0]
|
||||
on = render(cli, arw, os.path.join(work, name + ".on.png"), True)
|
||||
off = render(cli, arw, os.path.join(work, name + ".off.png"), False)
|
||||
out = {"file": name, "shot": exif_summary(jpg)}
|
||||
jgrid, jsize = board(jpg)
|
||||
out["kind"] = "board" if jgrid else "scene"
|
||||
if jgrid:
|
||||
out["straight_jpeg"] = straightness(jgrid)
|
||||
for tag, path in (("on", on), ("off", off)):
|
||||
if jgrid:
|
||||
grid, size = board(path)
|
||||
if not grid:
|
||||
continue
|
||||
k = jsize[0] / size[0]
|
||||
grid = {key: p * k for key, p in grid.items()}
|
||||
out["straight_" + tag] = straightness(grid)
|
||||
s, d = match_boards(grid, (size[0] * k, size[1] * k), jgrid, jsize)
|
||||
else:
|
||||
s, d = sift_pairs(path, jpg)
|
||||
if s is not None:
|
||||
out["vs_jpeg_" + tag] = agreement(s, d, jsize)
|
||||
if "vs_jpeg_on" in out and "vs_jpeg_off" in out and out["kind"] == "board":
|
||||
out["fraction"] = correction_fraction(out["vs_jpeg_on"], out["vs_jpeg_off"], jsize)
|
||||
for tag in ("on", "off"):
|
||||
out.get("vs_jpeg_" + tag, {}).pop("pairs", None)
|
||||
out["jpeg_width"] = jsize[0]
|
||||
return out
|
||||
|
||||
|
||||
def show(o):
|
||||
k = REPORT_WIDTH / o["jpeg_width"]
|
||||
print(f"\n{o['file']} {o['shot']} [{o['kind']}]")
|
||||
for tag, label in (("on", "LightCraft corrected"), ("off", "LightCraft uncorrected")):
|
||||
a = o.get("vs_jpeg_" + tag)
|
||||
if a:
|
||||
print(f" {label:24s} vs camera JPEG: median {a['median']:.2f} px, p95 {a['p95']:.2f}, max {a['max']:.2f}"
|
||||
f" (@{REPORT_WIDTH}: {a['median'] * k:.2f} / {a['p95'] * k:.2f}) n={a['n']} reach r={a['max_r']:.2f}")
|
||||
print(" by radius: " + " ".join(
|
||||
f"{b['r']}: {b['median']:.2f} ({b['radial']:+.2f})" for b in a["bins"] if b["n"]))
|
||||
for tag in ("jpeg", "on", "off"):
|
||||
st = o.get("straight_" + tag)
|
||||
if st:
|
||||
print(f" straightness {tag:4s}: worst line {st['max']:.2f} px, median line {st['median']:.2f} px ({st['lines']} lines)")
|
||||
if o.get("fraction") is not None:
|
||||
print(f" camera correction / LightCraft correction: {o['fraction']:.3f} (1 = same amount)")
|
||||
|
||||
|
||||
def main():
|
||||
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
|
||||
ap.add_argument("--cli", required=True)
|
||||
ap.add_argument("--work", required=True)
|
||||
ap.add_argument("--json")
|
||||
ap.add_argument("dirs", nargs="+")
|
||||
args = ap.parse_args()
|
||||
os.makedirs(args.work, exist_ok=True)
|
||||
arws = sorted(os.path.join(d, f) for d in args.dirs for f in os.listdir(d) if f.lower().endswith(".arw"))
|
||||
results = []
|
||||
for arw in arws:
|
||||
try:
|
||||
o = measure(args.cli, args.work, arw)
|
||||
except Exception as e: # keep going: one bad file must not stop the batch
|
||||
print(f"\n{arw}: {e}", file=sys.stderr)
|
||||
continue
|
||||
if o:
|
||||
show(o)
|
||||
results.append(o)
|
||||
if args.json:
|
||||
with open(args.json, "w") as f:
|
||||
json.dump(results, f, indent=1)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in New Issue
Block a user