Source file Vec3.ml

(* Claude Code
 *
 * Copyright (C) 2026 Yoann Padioleau
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public License
 * (LGPL) as published by the Free Software Foundation; either version
 * 2 of the License, or (at your option) any later version.
 *)

(* See Vec3.mli *)

type t = float * float * float

let add ((ax, ay, az) : t) ((bx, by, bz) : t) : t = (ax +. bx, ay +. by, az +. bz)
let sub ((ax, ay, az) : t) ((bx, by, bz) : t) : t = (ax -. bx, ay -. by, az -. bz)
let scale (s : float) ((x, y, z) : t) : t = (s *. x, s *. y, s *. z)
let dot ((ax, ay, az) : t) ((bx, by, bz) : t) : float = (ax *. bx) +. (ay *. by) +. (az *. bz)

let cross ((ax, ay, az) : t) ((bx, by, bz) : t) : t =
  ((ay *. bz) -. (az *. by), (az *. bx) -. (ax *. bz), (ax *. by) -. (ay *. bx))

let length (v : t) : float = sqrt (dot v v)

let normalize (v : t) : t =
  let n = length v in
  if n = 0. then v else scale (1. /. n) v

let centroid (points : t list) : t =
  let n = float_of_int (List.length points) in
  scale (1. /. n) (List.fold_left add (0., 0., 0.) points)

(* claude: bugfix: spheres drawn "cut" at the top.
 *
 * The symptom: every sphere was missing its top cap, the ring of faces
 * around the north pole; you could see through it to the white
 * background. On the software and web backends; the bottom cap was
 * fine, the OpenGL backend was fine, and it didn't depend on the
 * shading mode.
 *
 * The cause: the face normal used to be computed from the polygon's
 * first three points, as the cross product of its first two edges:
 *
 *   normalize (cross (sub p1 p0) (sub p2 p0))
 *
 * That's right for any real triangle, but Playground3d.sphere builds
 * its faces as latitude/longitude quads, and at the pole a "quad" has
 * two corners at the very same point, the pole itself:
 *
 *           p0 = p1 = pole (0, 1, 0)       the top row of faces:
 *                /\                        each quad's first two
 *               /  \                       corners are both the
 *              /    \                      pole, so it's really a
 *            p3------p2                    triangle p0 p2 p3
 *
 * so p1 - p0 = (0, 0, 0), and the cross product is (0, 0, 0). And
 * [normalize] leaves the zero vector as it is (it checks for a zero
 * length, to not divide by 0): the "normal" is (0, 0, 0). The bottom
 * row is fine because there the pole comes *last* (p2 = p3), after
 * three distinct points.
 *
 * Why the faces vanished instead of, say, being drawn black: the
 * backface test is "draw if dot normal (eye - centroid) > 0.", and the
 * dot product of anything with (0, 0, 0) is exactly 0, which is not
 * > 0: every top face was culled, silently -- no exception, no warning,
 * just missing faces. The OpenGL backend didn't have the bug: the GPU
 * culls each *triangle* by the order its corners appear on screen,
 * clockwise or not, and never computes a face normal.
 *
 * The fix: Newell's method. Each coordinate of the normal is computed
 * from *all* the edges, going around the polygon:
 *
 *   nx = sum over the edges (p, q) of (p.y - q.y) * (p.z + q.z)
 *   ny = sum over the edges (p, q) of (p.z - q.z) * (p.x + q.x)
 *   nz = sum over the edges (p, q) of (p.x - q.x) * (p.y + q.y)
 *
 * Each sum is the shoelace formula (see graphics/2d/Stroke.signed_area):
 * twice the signed area of the polygon's shadow on one of the three
 * coordinate planes (nz: on the xy plane, looking down z), and those
 * three areas together are the plane's direction, scaled by the
 * polygon's area. An edge from a point to the same point adds exactly
 * 0 to every sum, so a repeated pole changes nothing, and a triangle
 * with a repeated corner gets the normal of the triangle it really is.
 * For a counterclockwise triangle it agrees with the cross product,
 * e.g. (0,0,0), (1,0,0), (0,1,0):
 *
 *   nz = (0-1)*(0+0) + (1-0)*(0+1) + (0-0)*(1+0) = 1, so (0, 0, 1)
 *
 * and for a quad that isn't quite flat (4 points need not lie in one
 * plane) it gives a sensible average normal, where the first three
 * points would only see one corner of it.
 *
 * The lesson, i.e. why such a bug survives: the old code is textbook
 * code, and it *looks* obviously right. Its hidden assumption is that
 * the first three points are distinct and not on one line. That's true
 * of every polygon written by hand (a cube's faces, a triangle), which
 * is how it got tested, and false for polygons generated by formulas,
 * where degenerate cases come naturally: a sphere's pole, a cylinder's
 * or cone's tip, or a quad with a corner in the middle of an edge (the
 * first three points in a straight line: their cross product is
 * (0, 0, 0) too). And the failure is silent, made so by a guard meant
 * to be "safe": without [normalize]'s zero check, the division by 0
 * would have given NaNs, just as silent here; a [failwith] on a zero
 * normal would have pointed right at it. So: a degenerate case that
 * "can't happen" deserves an error, not a quiet default. And when a
 * picture looks slightly off, find out why: the notch was visible in
 * screenshots for a while and was dismissed as "some rendering
 * artifact" before anyone asked. (Even this comment first said
 * the normal was NaN, until someone read [normalize] again.)
 *
 * Found by eye, on a screenshot of Spheres3d.exe.
 *
 * Reference: Martin Newell's method, as described in Ivan Sutherland,
 * Robert Sproull, Robert Schumacker, "A Characterization of Ten
 * Hidden-Surface Algorithms", ACM Computing Surveys 6(1):1-55, 1974. *)
let face_normal (points : t list) : t =
  match points with
  | [] | [ _ ] | [ _; _ ] -> failwith "a polygon needs at least 3 points"
  | first :: rest ->
      (* each point with the next one, the last with the first *)
      let edges = List.combine points (rest @ [ first ]) in
      normalize
        (List.fold_left
           (fun (nx, ny, nz) ((x0, y0, z0), (x1, y1, z1)) ->
             ( nx +. ((y0 -. y1) *. (z0 +. z1)),
               ny +. ((z0 -. z1) *. (x0 +. x1)),
               nz +. ((x0 -. x1) *. (y0 +. y1)) ))
           (0., 0., 0.) edges)