#!/usr/bin/env catnip
# Rotors 3D, caméra et rendu ombré (algèbre géométrique)
#
# Un rotor R = s + B (scalaire + bivecteur) tourne un vecteur par le produit
# sandwich v' = R v R̃, sans construire de matrice 4x4. Le bivecteur B encode le
# PLAN de rotation ; son dual est l'axe. Ici B = sin(θ/2)·k, s = cos(θ/2), ce qui
# est l'isomorphisme rotor ↔ quaternion unitaire — le sandwich se réduit alors à
# v' = v + 2s(B×v) + 2 B×(B×v).
#
# Catnip organise sommets, faces, preset caméra et export ; numpy effectue les
# produits croisés, les normalisations et la matrice de contrôle. La sortie est
# une image PNG et une erreur numérique max (rotor vs matrice de rotation).
#
# DEPS: numpy, pillow

numpy = import('numpy')
math = import('math')
Image = import('PIL.Image')
ImageDraw = import('PIL.ImageDraw')
import('pathlib', 'Path')

script_dir = Path(META.file).parent
output_dir = script_dir / 'output'
output_dir.mkdir(exist_ok=True)

WIDTH = 640
HEIGHT = 640
CENTER_X = WIDTH / 2.0
CENTER_Y = HEIGHT / 2.0
FOCAL = 620.0
NEAR = 0.1

struct Vec3 {
    x: float; y: float; z: float;

    arr(self) => { numpy.array(list(self.x, self.y, self.z)) }
}

from_arr = (a): Vec3 => { Vec3(float(a[0]), float(a[1]), float(a[2])) }

# Rotor = partie scalaire + bivecteur (les composantes du bivecteur sont duales
# de l'axe de rotation). Construit depuis un axe et un angle.
struct Rotor {
    s: float; b: Vec3;

    # Produit sandwich R v R̃, écrit sous forme fermée avec deux produits croisés.
    apply(self, v: Vec3): Vec3 => {
        bv = self.b.arr()
        t = 2.0 * numpy.cross(bv, v.arr())
        from_arr(v.arr() + self.s * t + numpy.cross(bv, t))
    }
}

rotor_from_axis_angle = (axis: Vec3, angle: float): Rotor => {
    k = axis.arr() / numpy.linalg.norm(axis.arr())
    half = angle / 2.0
    Rotor(math.cos(half), from_arr(math.sin(half) * k))
}

# Matrice de rotation de référence (Rodrigues) pour le contrôle numérique.
# Chemin totalement indépendant du rotor : même (axe, angle), formule classique.
rotation_matrix = (axis: Vec3, angle: float) => {
    k = axis.arr() / numpy.linalg.norm(axis.arr())
    cross_k = numpy.array(list(
            list(0.0, -k[2], k[1]),
            list(k[2], 0.0, -k[0]),
            list(-k[1], k[0], 0.0),
        ))
    numpy.eye(3) * math.cos(angle) + math.sin(angle) * cross_k + (1.0 - math.cos(angle)) * numpy.outer(k, k)
}

struct Face {
    indices: list[int]; color: list[int];
}

struct CameraPreset {
    name: str; axis: Vec3; angle: float; distance: float; light: Vec3;
}

# Cube centré à l'origine. Chaque face est orientée CCW vue de l'extérieur, si
# bien que cross(p1-p0, p2-p0) donne la normale sortante.
cube_vertices = list(
    Vec3(-1.0, -1.0, -1.0),
    Vec3(1.0, -1.0, -1.0),
    Vec3(1.0, 1.0, -1.0),
    Vec3(-1.0, 1.0, -1.0),
    Vec3(-1.0, -1.0, 1.0),
    Vec3(1.0, -1.0, 1.0),
    Vec3(1.0, 1.0, 1.0),
    Vec3(-1.0, 1.0, 1.0),
)

cube_faces = list(
    Face(list(4, 5, 6, 7), list(220, 60, 60)),  # +z rouge
    Face(list(0, 3, 2, 1), list(60, 180, 75)),  # -z vert
    Face(list(1, 2, 6, 5), list(60, 110, 220)),  # +x bleu
    Face(list(0, 4, 7, 3), list(230, 200, 50)),  # -x jaune
    Face(list(3, 7, 6, 2), list(200, 70, 200)),  # +y magenta
    Face(list(0, 1, 5, 4), list(60, 200, 210)),  # -y cyan
)

# Le nom est séparé des paramètres : le match sélectionne une configuration
# complète et empêche les presets partiels (même forme que les autres exemples).
preset_name = 'three-quarter'
cam = match preset_name {
    'three-quarter' => {
        CameraPreset('three-quarter', Vec3(1.0, 1.0, 0.3), 0.9, 5.0, Vec3(-0.4, 0.7, 1.0))
    }
    'tilt' => {
        CameraPreset('tilt', Vec3(1.0, 0.0, 0.0), 0.6, 5.0, Vec3(0.0, 0.5, 1.0))
    }
    'spin' => {
        CameraPreset('spin', Vec3(0.0, 1.0, 0.0), 1.2, 5.0, Vec3(-0.5, 0.4, 1.0))
    }
    _ => {
        CameraPreset('front', Vec3(0.0, 1.0, 0.0), 0.0, 5.0, Vec3(0.0, 0.0, 1.0))
    }
}

rotor = rotor_from_axis_angle(cam.axis, cam.angle)

print(f"⇒ Preset « {cam.name} »")
print(f"  rotor s={round(rotor.s, 4)}, bivecteur=({round(rotor.b.x, 4)}, {round(rotor.b.y, 4)}, {round(rotor.b.z, 4)})")

# Le rotor est diffusé sur tous les sommets : une seule expression, pas de boucle.
rotated = cube_vertices.[(v) => { rotor.apply(v) }]

# Contrôle numérique : la même rotation par matrice de référence, erreur max sur
# les 8 sommets. Le sandwich (via rotated) et Rodrigues sont deux chemins
# indépendants ; on réutilise rotated plutôt que de rappeler le rotor.
ref_matrix = rotation_matrix(cam.axis, cam.angle)
ref_points = numpy.array(cube_vertices.[(v) => { v.arr() }]).dot(ref_matrix.T)
rotor_points = numpy.array(rotated.[(v) => { v.arr() }])
max_error = float(numpy.max(numpy.linalg.norm(ref_points - rotor_points, axis=1)))
print(f"  contrôle rotor vs matrice : erreur max = {max_error}")

# Projection perspective : caméra en (0, 0, distance) regardant l'origine.
# depth est borné à NEAR pour éviter la division par un plan derrière la caméra
# si un preset rapproche la caméra sous l'étendue de l'objet.
project = (v: Vec3) => {
    raw = cam.distance - v.z
    depth = if raw > NEAR { raw } else { NEAR }
    tuple(CENTER_X + FOCAL * v.x / depth, CENTER_Y - FOCAL * v.y / depth)
}
projected = rotated.[(v) => { project(v) }]

cam_pos = numpy.array(list(0.0, 0.0, cam.distance))
light_dir = cam.light.arr() / numpy.linalg.norm(cam.light.arr())

# Normale sortante d'une face, dans l'espace après rotation.
face_normal = (face: Face) => {
    p0 = rotated[face.indices[0]].arr()
    p1 = rotated[face.indices[1]].arr()
    p2 = rotated[face.indices[2]].arr()
    n = numpy.cross(p1 - p0, p2 - p0)
    n / numpy.linalg.norm(n)
}

face_centroid = (face: Face) => {
    pts = face.indices.[(i) => { rotated[i].arr() }]
    reduce(pts, (a, b) => { a + b }) / len(pts)
}

# Éclairage lambertien : une composante ambiante fixe plus un terme diffus n·l.
shade = (base: list[int], intensity: float) => {
    base.[(c) => { int(numpy.clip(c * intensity, 0, 255)) }]
}

image = Image.new('RGB', tuple(WIDTH, HEIGHT), tuple(18, 18, 24))
draw = ImageDraw.Draw(image)

# Back-face culling au tracé : une face dont la normale ne pointe pas vers la
# caméra n'est pas dessinée. Le cube est convexe et les faces visibles ne se
# recouvrent pas en projection, donc l'ordre de tracé n'a pas d'effet (pas de
# tri peintre nécessaire).
draw_face = (face: Face): int => {
    normal = face_normal(face)
    if numpy.dot(normal, cam_pos - face_centroid(face)) > 0.0 {
        intensity = 0.25 + 0.75 * float(numpy.maximum(0.0, numpy.dot(normal, light_dir)))
        rgb = shade(face.color, intensity)
        polygon = face.indices.[(i) => { projected[i] }]
        draw.polygon(polygon, fill=tuple(rgb[0], rgb[1], rgb[2]), outline=tuple(15, 15, 20))
        1
    } else {
        0
    }
}
drawn = cube_faces.[(f) => { draw_face(f) }]
visible_count = int(numpy.sum(numpy.array(drawn)))

print(f"  faces visibles : {visible_count} / {len(cube_faces)}")

output_path = output_dir / 'rotor_camera.png'
image.save(str(output_path))
print()
print(f"⇒ Image → {output_path}")