"""Fit a pinhole camera to the user's concept picture (reference/user_concept.jpg, 1408 x 768) from landmarks, so the model can be laid over it and compared (a measuring tool only; the concept's scene is not built). python build/ae_concept_fit.py -> build/concept_cam.json (Blender camera: location, rotation XYZ, lens for a 36 mm sensor) Landmarks: model points (from the parts' own geometry) and where the concept shows them (read from the 4x ruler zooms reference/z4_c*.png). Levenberg-Marquardt over position, rotation and focal length; principal point at the centre. """ import os import sys import json import math import numpy as np HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, HERE) import ae_lib as L import ae_parts as A W_IMG, H_IMG = 1408, 768 def mat_centroids(mesh, mat): V = np.array(mesh.v) return np.array([V[list(f)].mean(0) for f, m in zip(mesh.f, mesh.m) if m == mat]) def cluster(P, k_expected): """Split points into k groups along their largest spread (simple 1D k-means on a projection).""" P = np.asarray(P) d = P - P.mean(0) _, _, vt = np.linalg.svd(d, full_matrices=False) t = d @ vt[0] cuts = np.quantile(t, np.linspace(0, 1, k_expected + 1)) out = [] order = np.argsort(t) for g in np.array_split(order, k_expected): out.append(P[g].mean(0)) return out def model_marks(): m = {} belt = A.nose_belt() lamps = mat_centroids(belt, 'Lamp') # four lamp discs (two per block): group by side, then by height for side, sg in (('port', 1), ('stbd', -1)): s = lamps[np.sign(lamps[:, 0]) == sg] hi = s[s[:, 2] > np.median(s[:, 2])].mean(0); lo = s[s[:, 2] <= np.median(s[:, 2])].mean(0) m[f'belt_lamp_{side}_up'] = hi; m[f'belt_lamp_{side}_lo'] = lo ch = mat_centroids(A.chin(), 'Lamp') m['chin_lamp_stbd'] = ch[ch[:, 0] < 0].mean(0); m['chin_lamp_port'] = ch[ch[:, 0] > 0].mean(0) sl = mat_centroids(A.side_lamp(), 'Lamp').mean(0) m['side_lamp_port'] = sl m['fan_low_port'] = L.to_world(L.NACELLE['X'], 312.0, L.NACELLE['y']) m['fan_low_stbd'] = L.to_world(-L.NACELLE['X'], 312.0, L.NACELLE['y']) m['fan_up_port'] = L.to_world(L.POD['X'], 548.0, 274.0) m['fin_top_fore'] = L.to_world(0, 577.0, 96.0) m['fin_top_aft'] = L.to_world(0, 600.0, 96.0) m['fin_root_fore'] = L.to_world(0, 516.0, 136.0) m['hatch_top'] = L.to_world(0, L.HATCH[0], 116.0) def skin(v, y): vv, ph = L.side_to_param(np.array([v]), np.array([y])) return L.surf_world(vv, ph)[0] bx = L.TEXT_NAME['box'] m['name_A_bottom'] = skin(bx[0], bx[3]) m['name_9_bottom'] = skin(bx[2], bx[3]) m['name_A_top'] = skin(bx[0], bx[1]) return m # concept pixels (full image), from the 4x zooms IMG = { 'belt_lamp_stbd_up': (671.25, 260.0), 'belt_lamp_stbd_lo': (673.75, 271.25), 'belt_lamp_port_up': (755.0, 276.25), 'belt_lamp_port_lo': (758.0, 287.5), 'chin_lamp_stbd': (731.25, 365.0), 'chin_lamp_port': (753.75, 371.25), 'side_lamp_port': (861.5, 281.5), 'fan_low_port': (857.5, 357.5), 'fan_low_stbd': (687.5, 321.0), 'fan_up_port': (1025.0, 265.0), 'fin_top_fore': (946.25, 83.75), 'fin_top_aft': (953.0, 83.75), 'fin_root_fore': (905.0, 128.75), 'hatch_top': (872.5, 116.25), 'name_A_bottom': (915.0, 172.5), 'name_9_bottom': (986.25, 161.25), 'name_A_top': (915.5, 160.0), } def euler_xyz(rx, ry, rz): cx, sx, cy, sy, cz, sz = math.cos(rx), math.sin(rx), math.cos(ry), math.sin(ry), math.cos(rz), math.sin(rz) Rx = np.array([[1, 0, 0], [0, cx, -sx], [0, sx, cx]]) Ry = np.array([[cy, 0, sy], [0, 1, 0], [-sy, 0, cy]]) Rz = np.array([[cz, -sz, 0], [sz, cz, 0], [0, 0, 1]]) return Rz @ Ry @ Rx # Blender 'XYZ': X first def project(p, X3): C = p[:3]; R = euler_xyz(*p[3:6]); f = math.exp(p[6]) pc = (X3 - C) @ R # world -> camera (camera looks down -Z, up +Y) x = W_IMG / 2 + f * pc[:, 0] / -pc[:, 2] y = H_IMG / 2 - f * pc[:, 1] / -pc[:, 2] return np.stack([x, y], 1) def look_euler(cam, tgt): d = np.asarray(tgt, float) - np.asarray(cam, float); d /= np.linalg.norm(d) # camera -Z along d, +Y up-ish: yaw about Z, pitch about X rx = math.acos(np.clip(-d[2], -1, 1)) # tilt from looking straight down rz = math.atan2(-d[0], d[1]) return np.array([rx, 0.0, rz]) def fit(): mm = model_marks() names = [n for n in IMG if n in mm] X3 = np.array([mm[n] for n in names]); UV = np.array([IMG[n] for n in names], float) best = None for az in range(-150, -10, 10): # camera round the nose-port quarter for el in (5, 15, 25): for dist in (8.0, 11.0, 15.0): tgt = np.array([0, -0.5, 1.4]) a, e = math.radians(az), math.radians(el) cam = tgt + dist * np.array([math.cos(e) * math.cos(a), math.cos(e) * math.sin(a), math.sin(e)]) cam[0] = abs(cam[0]) # the port side faces the camera p = np.concatenate([cam, look_euler(cam, tgt), [math.log(1400.0)]]) p = lm(p, X3, UV) r = np.sqrt(((project(p, X3) - UV) ** 2).sum(1)) if best is None or r.mean() < best[1]: best = (p, r.mean(), r) p, mean, r = best return p, dict(zip(names, np.round(r, 2).tolist())), mean, names, X3, UV def lm(p, X3, UV, iters=200): lam = 1e-2 def res(q): return (project(q, X3) - UV).ravel() r = res(p); cost = r @ r for _ in range(iters): J = np.zeros((len(r), len(p))) for k in range(len(p)): dq = np.zeros(len(p)); dq[k] = 1e-6 if k < 3 else 1e-7 J[:, k] = (res(p + dq) - r) / dq[k] H = J.T @ J; g = J.T @ r step = np.linalg.solve(H + lam * np.diag(np.diag(H) + 1e-9), -g) pn = p + step; rn = res(pn); cn = rn @ rn if cn < cost: p, r, cost = pn, rn, cn; lam *= 0.4 if np.linalg.norm(step) < 1e-9: break else: lam *= 4 return p if __name__ == '__main__': p, per, mean, names, X3, UV = fit() f_px = math.exp(p[6]) lens = f_px * 36.0 / W_IMG out = dict(location=p[:3].tolist(), rotation_euler_xyz=p[3:6].tolist(), focal_px=f_px, lens_mm_36=lens, resolution=[W_IMG, H_IMG], mean_px=round(float(mean), 2), per_mark_px=per) json.dump(out, open(os.path.join(HERE, 'concept_cam.json'), 'w'), indent=1) print(json.dumps(out, indent=1))