#!/usr/bin/env catnip
# Aires orientées, bivecteurs et détection de faces inversées
#
# Le produit extérieur (wedge) de deux arêtes d'un triangle forme un BIVECTEUR :
# un élément d'aire orientée. Pour un triangle (a, b, c), le bivecteur
# (a - c) ∧ (b - c) a pour magnitude 2×l'aire du triangle et pour « plan »
# l'orientation de la face. En 3D son dual est un vecteur, la normale :
# dual((a-c) ∧ (b-c)) = cross(a-c, b-c) = cross(b-a, c-a).
#
# On mesure l'aire orientée de chaque face d'un tétraèdre dont UNE face a ses
# sommets dans le mauvais ordre : sa normale pointe alors vers l'intérieur, ce
# qui la marque comme inversée. Contrôle numérique indépendant : l'aire tirée du
# bivecteur est confrontée à la formule de Héron (longueurs d'arêtes seules).
# La sortie est un rapport structuré et un PNG coloré par orientation.
#
# DEPS: numpy, pillow

numpy = import('numpy')
math = import('math')
Image = import('PIL.Image')
ImageDraw = import('PIL.ImageDraw')
builtins = import('builtins')
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 = 520.0
VIEW_DISTANCE = 5.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])) }

struct Face {
    indices: list[int]; name: str;
}

# Rapport d'une face : aire (via bivecteur), magnitude du bivecteur (= 2×aire),
# écart avec l'aire de Héron, drapeau d'inversion, intensité lambertienne, normale.
struct FaceReport {
    name: str; area: float; wedge_mag: float; area_error: float;
    inverted: bool; intensity: float; normal: Vec3;
}

struct MeshReport {
    total_area: float; global_orientation: str; suspect: list[str]; max_area_error: float;
}

# Matrice de rotation de vue (Rodrigues) : ne sert qu'au rendu. Une rotation de
# déterminant +1 préserve aires, magnitudes de bivecteurs et orientation, donc
# tout le calcul géométrique peut se faire dans l'espace vue.
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)
}

# Tétraèdre régulier centré à l'origine (le centroïde est donc l'origine).
tetra_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),
)

# Faces triangulaires en enroulement CCW vu de l'extérieur : cross(b-a, c-a)
# pointe alors vers l'extérieur. La face « base » est volontairement bruitée —
# ses sommets sont dans l'ordre (0, 2, 1) au lieu de (0, 1, 2), sa normale rentre.
tetra_faces = list(
    Face(list(1, 3, 2), 'front'),
    Face(list(0, 2, 3), 'left'),
    Face(list(0, 3, 1), 'right'),
    Face(list(0, 2, 1), 'base'),
)

VIEW_AXIS = Vec3(1.0, 0.0, 1.0)
VIEW_ANGLE = 2.2
view = rotation_matrix(VIEW_AXIS, VIEW_ANGLE)

# Sommets dans l'espace vue. Tout le reste travaille sur ces coordonnées.
rotated = tetra_vertices.[(v) => { from_arr(view.dot(v.arr())) }]
mesh_centroid = reduce(rotated.[(v) => { v.arr() }], (a, b) => { a + b }) / len(rotated)

# Lumière fixe dans l'espace vue (haut-gauche-avant), pour l'éclairage lambertien.
light_dir = Vec3(-0.4, 0.7, 1.0).arr() / numpy.linalg.norm(Vec3(-0.4, 0.7, 1.0).arr())

# Aire orientée d'une face via le bivecteur des arêtes (a-c) et (b-c). Son dual
# est la normale non unitaire ; sa magnitude vaut 2×l'aire. On confronte l'aire
# ainsi obtenue à celle de Héron (chemin indépendant : longueurs d'arêtes seules).
make_report = (face: Face): FaceReport => {
    a = rotated[face.indices[0]].arr()
    b = rotated[face.indices[1]].arr()
    c = rotated[face.indices[2]].arr()
    dual = numpy.cross(a - c, b - c)  # dual du bivecteur = normale sortante attendue
    wedge_mag = float(numpy.linalg.norm(dual))  # |bivecteur| = 2 × aire
    area = wedge_mag / 2.0
    la = float(numpy.linalg.norm(a - b))
    lb = float(numpy.linalg.norm(b - c))
    lc = float(numpy.linalg.norm(c - a))
    s = (la + lb + lc) / 2.0
    heron = math.sqrt(s * (s - la) * (s - lb) * (s - lc))
    area_error = float(numpy.abs(area - heron))
    centroid_face = (a + b + c) / 3.0
    # Face inversée : la normale du bivecteur pointe vers le centroïde du maillage
    # au lieu de s'en éloigner (produit scalaire négatif).
    dp = float(numpy.dot(dual, centroid_face - mesh_centroid))
    inverted = dp < 0.0
    normal_unit = dual / wedge_mag
    ndl = float(numpy.dot(normal_unit, light_dir))
    intensity = if ndl > 0.0 { ndl } else { 0.0 }  # Lambert, borné à 0 (face rentrante = sombre)
    FaceReport(face.name, area, wedge_mag, area_error, inverted, intensity, from_arr(normal_unit))
}

# Volume signé du maillage fermé : somme des volumes de tétraèdres (origine, a, b,
# c). Son signe donne l'orientation globale ; une face inversée le réduit sans en
# changer le signe tant qu'elle reste minoritaire.
face_signed_vol = (face: Face): float => {
    a = rotated[face.indices[0]].arr()
    b = rotated[face.indices[1]].arr()
    c = rotated[face.indices[2]].arr()
    float(numpy.dot(a, numpy.cross(b, c))) / 6.0
}

reports = tetra_faces.[(f) => { make_report(f) }]

report_for = (name: str): FaceReport => {
    i = 0
    found = reports[0]
    while i < len(reports) {
        if reports[i].name == name { found = reports[i] }
        i = i + 1
    }
    found
}

join = (items: list[str]): str => {
    out = ""
    i = 0
    while i < len(items) {
        sep = if i > 0 { ", " } else { "" }
        out = out + sep + items[i]
        i = i + 1
    }
    out
}

total_area = reduce(reports.[(r) => { r.area }], (a, b) => { a + b })
max_area_error = float(numpy.max(numpy.array(reports.[(r) => { r.area_error }])))
signed_volume = reduce(tetra_faces.[(f) => { face_signed_vol(f) }], (a, b) => { a + b })
orientation = if signed_volume > 0.0 { "sortante (CCW)" } else { "rentrante (CW)" }

suspect = list()
i = 0
while i < len(reports) {
    r = reports[i]
    if r.inverted { suspect.append(r.name) }
    i = i + 1
}

mesh_report = MeshReport(float(total_area), orientation, suspect, max_area_error)

print(f"⇒ Maillage : tétraèdre régulier, {len(tetra_faces)} faces triangulaires")
print()
print("⇒ Rapport par face (aire orientée via bivecteur)")
i = 0
while i < len(reports) {
    r = reports[i]
    flag = if r.inverted { "INVERSÉE" } else { "ok" }
    print(f"  {r.name:>6} : aire={round(r.area, 4):>7}  |bivecteur|={round(r.wedge_mag, 4):>7}  n·l={round(r.intensity, 3)}  {flag}")
    i = i + 1
}

suspect_str = if len(mesh_report.suspect) > 0 { join(mesh_report.suspect) } else { "aucune" }

print()
print("⇒ Rapport maillage")
print(f"  aire totale         : {round(mesh_report.total_area, 4)}")
print(f"  orientation globale : {mesh_report.global_orientation}")
print(f"  faces suspectes     : {suspect_str}")
print(f"  contrôle bivecteur vs Héron : erreur max = {mesh_report.max_area_error}")

# Projection perspective : caméra en (0, 0, VIEW_DISTANCE) regardant l'origine.
project = (v: Vec3) => {
    raw = VIEW_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) }]

# Profondeur d'une face (z moyen) pour l'algorithme du peintre : on trace du plus
# lointain au plus proche, sans culling, pour que la face inversée reste visible.
face_depth = (face: Face): float => {
    zs = face.indices.[(idx) => { float(rotated[idx].z) }]
    reduce(zs, (a, b) => { a + b }) / len(zs)
}

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

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

# Vert = face correctement orientée, rouge = face inversée ; l'intensité module.
ordered = builtins.sorted(tetra_faces, key=(f) => { face_depth(f) })
i = 0
while i < len(ordered) {
    face = ordered[i]
    rep = report_for(face.name)
    base = if rep.inverted { list(220, 70, 70) } else { list(70, 200, 110) }
    rgb = shade(base, 0.5 + 0.5 * rep.intensity)
    polygon = face.indices.[(idx) => { projected[idx] }]
    draw.polygon(polygon, fill=tuple(rgb[0], rgb[1], rgb[2]), outline=tuple(235, 235, 240))
    i = i + 1
}

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

# Aperçu navigateur : sert l'image une fois puis rend la main.
# --no-browser garde un chemin headless (le PNG reste écrit ci-dessus).
if '--no-browser' not in import('sys').argv {
    http = import('http')
    b64 = import('base64').b64encode(output_path.read_bytes()).decode('ascii')
    http.serve(f'<!doctype html><meta charset="utf-8"><title>Bivector mesh orientation</title><body style="margin:0;background:#0d1117"><img style="max-width:100%;display:block;margin:0 auto" src="data:image/png;base64,{b64}">',
        0, 'text/html; charset=utf-8', True)
}