5 portraits placed along a cubic Bézier curve

// 5 portraits placed along a cubic Bézier curve
// demo_bezier_placer.js
//!OUTPUT: OUTPUT

// demo_bezier_placer.js
// Places multiple copies of an image tangent-aligned along a cubic Bézier curve.
//
// Algorithm:
//  1. Build an arc-length LUT for uniform spacing along the curve.
//  2. At each placement, compute tangent + normal from the LUT.
//  3. Derive the 4-corner quad for the image in canvas space.
//  4. Solve the dst→src homography (canvas quad → portrait corners) directly,
//     using coordinate normalization to ensure numerical stability even for
//     very large source images (e.g. 2048×2048 portrait).
//  5. Composite each warped tile onto the canvas using Porter-Duff "Over" (mode 10).

const CANVAS_W = 1200;
const CANVAS_H = 600;
const NUM_IMAGES = 5;
const TILE_W = 200; // rendered width on canvas (pixels along tangent)
const TILE_H = 133; // rendered height on canvas (pixels along normal)
const BLEND_OVER = BlendMode.Over; // Porter-Duff Over mode

// Cubic Bézier control points (S-curve across the canvas)
const P0 = { x: 80,   y: 300 };
const P1 = { x: 380,  y: 30  };
const P2 = { x: 820,  y: 570 };
const P3 = { x: 1120, y: 300 };

// ── Bézier math ──────────────────────────────────────────────────────────────

function bezierPoint(t) {
  const mt = 1 - t, mt2 = mt * mt, mt3 = mt2 * mt;
  const t2 = t * t,  t3  = t2 * t;
  return {
    x: mt3*P0.x + 3*mt2*t*P1.x + 3*mt*t2*P2.x + t3*P3.x,
    y: mt3*P0.y + 3*mt2*t*P1.y + 3*mt*t2*P2.y + t3*P3.y,
  };
}

function bezierTangent(t) {
  const mt = 1 - t;
  const dx = 3*mt*mt*(P1.x-P0.x) + 6*mt*t*(P2.x-P1.x) + 3*t*t*(P3.x-P2.x);
  const dy = 3*mt*mt*(P1.y-P0.y) + 6*mt*t*(P2.y-P1.y) + 3*t*t*(P3.y-P2.y);
  const len = Math.sqrt(dx*dx + dy*dy);
  return len < 1e-6 ? { x: 1, y: 0 } : { x: dx/len, y: dy/len };
}

// ── Arc-length LUT ───────────────────────────────────────────────────────────

const LUT_STEPS = 400;
const arcLUT = [{ t: 0, dist: 0 }];
let totalLen = 0;
{
  let prev = bezierPoint(0);
  for (let i = 1; i <= LUT_STEPS; i++) {
    const t  = i / LUT_STEPS;
    const pt = bezierPoint(t);
    const dx = pt.x - prev.x, dy = pt.y - prev.y;
    totalLen += Math.sqrt(dx*dx + dy*dy);
    arcLUT.push({ t, dist: totalLen });
    prev = pt;
  }
}

/** Return the Bézier parameter t that corresponds to arc-length `d`. */
function getT(d) {
  if (d <= 0) { return 0; }
  if (d >= totalLen) { return 1; }
  let lo = 0, hi = arcLUT.length - 1;
  while (hi - lo > 1) {
    const mid = (lo + hi) >> 1;
    if (arcLUT[mid].dist < d) { lo = mid; } else { hi = mid; }
  }
  const s = (d - arcLUT[lo].dist) / (arcLUT[hi].dist - arcLUT[lo].dist);
  return arcLUT[lo].t + s * (arcLUT[hi].t - arcLUT[lo].t);
}

// ── Homography math ──────────────────────────────────────────────────────────

/** Solve Ax = b for an 8×8 system via Gauss-Jordan elimination. */
function gaussElim(A, b) {
  const n = 8;
  const M = A.map((row, i) => [...row, b[i]]);
  for (let col = 0; col < n; col++) {
    // Partial pivoting for numerical stability
    let maxRow = col;
    for (let row = col + 1; row < n; row++) {
      if (Math.abs(M[row][col]) > Math.abs(M[maxRow][col])) { maxRow = row; }
    }
    const tmp = M[col]; M[col] = M[maxRow]; M[maxRow] = tmp;
    for (let row = 0; row < n; row++) {
      if (row === col) { continue; }
      const f = M[row][col] / M[col][col];
      for (let j = col; j <= n; j++) { M[row][j] -= f * M[col][j]; }
    }
  }
  return M.map((row, i) => row[n] / row[i]);
}

/** 3×3 matrix multiply (flat 9-element row-major). */
function mat3mul(A, B) {
  const [a0,a1,a2, a3,a4,a5, a6,a7,a8] = A;
  const [b0,b1,b2, b3,b4,b5, b6,b7,b8] = B;
  return [
    a0*b0+a1*b3+a2*b6,  a0*b1+a1*b4+a2*b7,  a0*b2+a1*b5+a2*b8,
    a3*b0+a4*b3+a5*b6,  a3*b1+a4*b4+a5*b7,  a3*b2+a4*b5+a5*b8,
    a6*b0+a7*b3+a8*b6,  a6*b1+a7*b4+a8*b7,  a6*b2+a7*b5+a8*b8,
  ];
}

/**
 * Compute the src→dst homography that `warpPerspective` expects.
 * Kornia internally inverts the forward (src→dst) matrix to do backward-mapping.
 *
 * Maps each portrait corner (srcPts[i]) to the corresponding canvas quad
 * corner (dstPts[i]).
 *
 * Uses coordinate normalization to keep the 8×8 linear system well-conditioned
 * even for very large source images (e.g. 2048×2048 portraits).
 */
function computeH_fwd(srcPts, dstPts) {
  // srcPts: portrait corners (0,0),(sw,0),(sw,sh),(0,sh)  — coords up to 2047
  // dstPts: canvas quad corners                           — coords up to 1200
  const sw = srcPts[1].x;   // source width  - 1
  const sh = srcPts[2].y;   // source height - 1

  // Normalize src to [0,1] and dst to [0,1] for numerical stability.
  // Without normalization, matrix entries reach ~2047×1200 ≈ 2.5M causing
  // catastrophic cancellation in Gaussian elimination.
  const sn = srcPts.map(p => ({ x: p.x / sw,       y: p.y / sh       }));
  const dn = dstPts.map(p => ({ x: p.x / CANVAS_W,  y: p.y / CANVAS_H }));

  // Solve H_n: normalized-portrait → normalized-canvas (forward transform)
  // For each correspondence (xs, ys) → (xd, yd) [all in [0,1]]:
  //   xs*h1 + ys*h2 + h3 - xd*xs*h7 - xd*ys*h8 = xd   (h9=1)
  //   xs*h4 + ys*h5 + h6 - yd*xs*h7 - yd*ys*h8 = yd
  const A = [], b = [];
  for (let i = 0; i < 4; i++) {
    const { x: xs, y: ys } = sn[i];   // normalized src
    const { x: xd, y: yd } = dn[i];   // normalized dst
    A.push([xs, ys, 1, 0,  0,  0, -xd*xs, -xd*ys]); b.push(xd);
    A.push([0,  0,  0, xs, ys, 1, -yd*xs, -yd*ys]); b.push(yd);
  }
  const h = gaussElim(A, b);
  // H_n: [0,1]-portrait → [0,1]-canvas
  const H_n = [h[0],h[1],h[2], h[3],h[4],h[5], h[6],h[7], 1.0];

  // Compose to get actual-portrait → actual-canvas:
  //   H_fwd = denorm_dst @ H_n @ norm_src
  //   norm_src   = diag(1/sw, 1/sh, 1)    (scale actual-portrait → [0,1])
  //   denorm_dst = diag(CW, CH, 1)         (scale [0,1]-canvas → actual-canvas)
  const normSrc   = [1/sw, 0, 0,          0, 1/sh, 0,          0, 0, 1];
  const denormDst = [CANVAS_W, 0, 0,      0, CANVAS_H, 0,      0, 0, 1];
  return mat3mul(denormDst, mat3mul(H_n, normSrc));
}

// ── Rendering ────────────────────────────────────────────────────────────────

const src = Engine.loadImage(INPUT);
const srcW = src.width;
const srcH = src.height;

// White canvas background (create 1×1 white, resize to full canvas size)
const canvas = Engine.createImage(1, 1);
canvas.setPixel(px(0, 0), new Pixel(1, 1, 1, 1));
canvas.resize(CANVAS_W, CANVAS_H);

// Leave 15 % margin at each end so images don't fall off the canvas edges
const margin    = 0.15 * totalLen;
const usableLen = totalLen - 2 * margin;

for (let i = 0; i < NUM_IMAGES; i++) {
  const dist = margin + (i / (NUM_IMAGES - 1)) * usableLen;
  const t    = getT(dist);
  const pos  = bezierPoint(t);
  const tang = bezierTangent(t);

  // Normal vector pointing "down" on screen (90° CW from tangent):
  //   if tang=(1,0)  → N_cw=(0,1)  = screen-down  ✓
  //   if tang=(0,1)  → N_cw=(1,0)  = screen-right
  const N = { x: tang.y, y: -tang.x };

  const hw = TILE_W / 2; // half-width  along tangent
  const hh = TILE_H / 2; // half-height along normal

  // Destination quad on the canvas:
  //   Image X-axis (left→right) aligns with tangent T
  //   Image Y-axis (top→bottom) aligns with -N (i.e. T rotated 90° CCW = upward)
  //   The quad is centered on `pos`.
  const TL = { x: pos.x - hw*tang.x + hh*N.x,  y: pos.y - hw*tang.y + hh*N.y };
  const TR = { x: pos.x + hw*tang.x + hh*N.x,  y: pos.y + hw*tang.y + hh*N.y };
  const BR = { x: pos.x + hw*tang.x - hh*N.x,  y: pos.y + hw*tang.y - hh*N.y };
  const BL = { x: pos.x - hw*tang.x - hh*N.x,  y: pos.y - hw*tang.y - hh*N.y };

  // Source corners (actual portrait pixel coordinates)
  const srcPts = [
    { x: 0,       y: 0       },  // → TL
    { x: srcW-1,  y: 0       },  // → TR
    { x: srcW-1,  y: srcH-1  },  // → BR
    { x: 0,       y: srcH-1  },  // → BL
  ];
  const dstPts = [TL, TR, BR, BL];

  // Numerically stable forward (src→dst) matrix; Kornia inverts it internally
  const H_fwd = computeH_fwd(srcPts, dstPts);

  // Warp source image into its canvas-sized slot, then composite Over
  const tile = src.clone().warpPerspective(CANVAS_W, CANVAS_H, H_fwd);
  canvas.blendAt(tile, px(0, 0), 1.0, BLEND_OVER);
  tile.free();
}

// ── Overlay: Bézier curve + control points (tiny-skia vector rendering) ─────

const ov = Engine.createCanvas(CANVAS_W, CANVAS_H);

// Helper: fill + stroke a circle — creates+frees the path internally
function circleMarker(cv, cx, cy, r, color, strokeW) {
  const p = Engine.createPath().circle(cx, cy, r);
  cv.fill(color).drawPath(p);
  if (strokeW > 0) { cv.pen("#ffffffff", strokeW).drawPath(p, true); }
  p.free();
}

// 1. The S-curve — single cubic Bézier, 3 px dark stroke, anti-aliased
ov.pen(new Pixel(0.10, 0.10, 0.10, 0.85), 3.0);
ov.drawPathStr(`M ${P0.x},${P0.y} C ${P1.x},${P1.y} ${P2.x},${P2.y} ${P3.x},${P3.y}`, true);

// 2. Handle tangent lines P0→P1 and P3→P2 (thin orange, 1.5 px)
ov.pen(new Pixel(0.92, 0.42, 0.12, 0.75), 1.5);
ov.drawPathStr(`M ${P0.x},${P0.y} L ${P1.x},${P1.y}`, true);
ov.drawPathStr(`M ${P3.x},${P3.y} L ${P2.x},${P2.y}`, true);

// 3. Anchor endpoints (P0, P3) — blue filled circles, white ring
for (const p of [P0, P3]) { circleMarker(ov, p.x, p.y, 7, new Pixel(0.18, 0.40, 0.92, 1.0), 1.5); }

// 4. Control handles (P1, P2) — orange filled circles, white ring
for (const p of [P1, P2]) { circleMarker(ov, p.x, p.y, 5, new Pixel(0.92, 0.42, 0.12, 0.90), 1.5); }

// 5. Image placement positions — green dots, white ring
for (let i = 0; i < NUM_IMAGES; i++) {
  const dist = margin + (i / (NUM_IMAGES - 1)) * usableLen;
  const pos  = bezierPoint(getT(dist));
  circleMarker(ov, pos.x, pos.y, 4, new Pixel(0.15, 0.72, 0.32, 1.0), 1.0);
}

// Composite vector overlay onto the photo canvas (Porter-Duff Over)
const ovImg = ov.toImage();
canvas.blendAt(ovImg, px(0, 0), 1.0, BLEND_OVER);
ovImg.free();
ov.free();

canvas.save(OUTPUT);

// © 2026 Michael Lechner · mlc OpticScript · https://mlcgo.eu · Elastic License 2.0