VoxelForge/Standalone/Research/WorldAssemblerM2/RenderSharpMesherCaptures.py

#!/usr/bin/env python3
"""Render the standalone spike's exported meshes from matched orthographic views."""
from __future__ import annotations

import math
from pathlib import Path

from PIL import Image, ImageDraw


ROOT = Path(__file__).resolve().parent
RESULTS = ROOT / "results"
CAPTURES = RESULTS / "captures"
CAPTURES.mkdir(parents=True, exist_ok=True)
METHODS = ("MC", "Hybrid", "DC")
PIECES = (
    ("corner90", "90 degree corner"),
    ("doorway", "doorway"),
    ("riser1m", "one metre riser"),
    ("cliff_overhang", "cliff ledge / overhang"),
    ("organic_cave_chamber", "organic cave chamber wall"),
    ("native_surface_cliff", "sampled SurfaceWorld cliff"),
)
WIDTH, HEIGHT = 288, 216
BG = (239, 242, 244)
INK = (36, 44, 51)
SEAM_FINE = (73, 133, 174)
SEAM_COARSE = (188, 112, 68)


def normalize(v):
    n = math.sqrt(sum(c * c for c in v))
    return tuple(c / n for c in v) if n > 1.0e-12 else (0.0, 0.0, 1.0)


def dot(a, b):
    return sum(x * y for x, y in zip(a, b))


def cross(a, b):
    return (a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0])


def read_obj(path: Path):
    vertices, normals, faces = [], [], []
    for line in path.read_text(encoding="ascii").splitlines():
        if line.startswith("v "):
            vertices.append(tuple(map(float, line.split()[1:4])))
        elif line.startswith("vn "):
            normals.append(tuple(map(float, line.split()[1:4])))
        elif line.startswith("f "):
            face = []
            for item in line.split()[1:4]:
                fields = item.split("//")
                face.append((int(fields[0]) - 1, int(fields[1]) - 1 if len(fields) > 1 else int(fields[0]) - 1))
            faces.append(tuple(face))
    return vertices, normals, faces


def project_mesh(vertices, faces, seam=False, target=None):
    # One view is reused for every candidate in a row: azimuth 38°, elevation 27°.
    if target is None:
        if vertices:
            bounds_min = tuple(min(v[i] for v in vertices) for i in range(3))
            bounds_max = tuple(max(v[i] for v in vertices) for i in range(3))
            target = tuple((bounds_min[i] + bounds_max[i]) * 0.5 for i in range(3))
        else:
            target = (0.0, 0.0, 0.0)
    view_dir = normalize((8.0, -10.0, 7.0))
    forward = tuple(-c for c in view_dir)
    right = normalize(cross(forward, (0.0, 0.0, 1.0)))
    up = normalize(cross(right, forward))
    projected_bounds = []
    for point in vertices:
        rel = tuple(point[i] - target[i] for i in range(3))
        projected_bounds.append((dot(rel, right), dot(rel, up)))
    span_x = max((point[0] for point in projected_bounds), default=0.0) - min((point[0] for point in projected_bounds), default=0.0)
    span_y = max((point[1] for point in projected_bounds), default=0.0) - min((point[1] for point in projected_bounds), default=0.0)
    scale = min(20.2, 0.84 * min((WIDTH - 24) / max(span_x, 1.0e-6), (HEIGHT - 24) / max(span_y, 1.0e-6)))
    projected = []
    for p in vertices:
        rel = tuple(p[i] - target[i] for i in range(3))
        projected.append((WIDTH * 0.5 + dot(rel, right) * scale,
                          HEIGHT * 0.51 - dot(rel, up) * scale,
                          dot(rel, forward)))
    return projected


def raster(path: Path, color_by_x=False):
    vertices, normals, faces = read_obj(path)
    if color_by_x:
        clipped_vertices, clipped_normals, clipped_faces = [], [], []
        for face in faces:
            center_x = sum(vertices[i][0] for i, _ in face) / 3.0
            left = center_x < 0.0
            polygon = [(vertices[vi], normals[ni] if ni < len(normals) else (0.0, 0.0, 1.0)) for vi, ni in face]
            clipped = []
            for i, (pa, na) in enumerate(polygon):
                pb, nb = polygon[(i + 1) % len(polygon)]
                inside_a = pa[0] <= 0.0 if left else pa[0] >= 0.0
                inside_b = pb[0] <= 0.0 if left else pb[0] >= 0.0
                if inside_a:
                    clipped.append((pa, na))
                if inside_a != inside_b:
                    t = (0.0 - pa[0]) / (pb[0] - pa[0])
                    point = tuple(pa[k] + t * (pb[k] - pa[k]) for k in range(3))
                    normal = normalize(tuple(na[k] + t * (nb[k] - na[k]) for k in range(3)))
                    clipped.append((point, normal))
            if len(clipped) < 3:
                continue
            first = len(clipped_vertices)
            for point, normal in clipped:
                clipped_vertices.append(point)
                clipped_normals.append(normal)
            for i in range(1, len(clipped) - 1):
                clipped_faces.append(((first, first), (first + i, first + i), (first + i + 1, first + i + 1)))
        vertices, normals, faces = clipped_vertices, clipped_normals, clipped_faces
    bounds_min = tuple(min(v[i] for v in vertices) for i in range(3))
    bounds_max = tuple(max(v[i] for v in vertices) for i in range(3))
    target = tuple((bounds_min[i] + bounds_max[i]) * 0.5 for i in range(3))
    projected = project_mesh(vertices, faces, color_by_x, target)
    pixels = Image.new("RGB", (WIDTH, HEIGHT), BG)
    zbuf = [float("inf")] * (WIDTH * HEIGHT)
    px = pixels.load()
    light = normalize((-0.45, -0.55, 0.75))
    for face in faces:
        ids = [v for v, _ in face]
        points = [projected[i] for i in ids]
        x0 = max(0, int(math.floor(min(p[0] for p in points))))
        x1 = min(WIDTH - 1, int(math.ceil(max(p[0] for p in points))))
        y0 = max(0, int(math.floor(min(p[1] for p in points))))
        y1 = min(HEIGHT - 1, int(math.ceil(max(p[1] for p in points))))
        if x0 > x1 or y0 > y1:
            continue
        (ax, ay, az), (bx, by, bz), (cx, cy, cz) = points
        denom = (by - cy) * (ax - cx) + (cx - bx) * (ay - cy)
        if abs(denom) < 1.0e-9:
            continue
        face_normals = [normals[ni] if ni < len(normals) else (0.0, 0.0, 1.0) for _, ni in face]
        for y in range(y0, y1 + 1):
            for x in range(x0, x1 + 1):
                sx, sy = x + 0.5, y + 0.5
                u = ((by - cy) * (sx - cx) + (cx - bx) * (sy - cy)) / denom
                v = ((cy - ay) * (sx - cx) + (ax - cx) * (sy - cy)) / denom
                w = 1.0 - u - v
                if u < -1.0e-8 or v < -1.0e-8 or w < -1.0e-8:
                    continue
                depth = u * az + v * bz + w * cz
                idx = y * WIDTH + x
                if depth >= zbuf[idx]:
                    continue
                zbuf[idx] = depth
                n = normalize(tuple(u * face_normals[0][k] + v * face_normals[1][k] + w * face_normals[2][k] for k in range(3)))
                shade = 0.36 + 0.64 * max(0.0, dot(n, light))
                if color_by_x:
                    center_x = sum(vertices[i][0] for i in ids) / 3.0
                    base = SEAM_FINE if center_x < 0.0 else SEAM_COARSE
                else:
                    base = (130, 145, 151)
                px[x, y] = tuple(max(0, min(255, int(c * shade))) for c in base)
    if color_by_x:
        draw = ImageDraw.Draw(pixels)
        # Project a visible segment on the LOD boundary plane, x=0.
        seam_world = [(-0.0, y, z) for y, z in ((-4, -4), (4, -4), (4, 4), (-4, 4), (-4, -4))]
        seam_points = project_mesh(seam_world, [], True, target)
        draw.line([(p[0], p[1]) for p in seam_points], fill=(181, 55, 49), width=2)
    return pixels


def capture_title(image: Image.Image, title: str, subtitle: str):
    canvas = Image.new("RGB", (WIDTH, HEIGHT + 30), BG)
    canvas.paste(image, (0, 30))
    draw = ImageDraw.Draw(canvas)
    draw.rectangle((0, 0, WIDTH, 29), fill=(31, 42, 50))
    draw.text((8, 4), title, fill=(255, 255, 255))
    draw.text((8, 16), subtitle, fill=(196, 210, 218))
    return canvas


def main():
    captures = []
    for key, label in PIECES:
        row = []
        for method in METHODS:
            obj = RESULTS / "meshes" / f"{key}_{method}.obj"
            image = capture_title(raster(obj), label, f"{method} | matched orthographic view")
            path = CAPTURES / f"{key}_{method}.jpg"
            image.save(path, quality=86, optimize=True)
            row.append(image)
        captures.append(row)
    seam_row = []
    for method in METHODS:
        obj = RESULTS / "meshes" / f"lodseam_cliff_crossing_{method}.obj"
        image = capture_title(raster(obj, color_by_x=True), "Sharp cliff across LOD border", f"{method} | LOD0 left, LOD1 right, x=0")
        path = CAPTURES / f"lodseam_cliff_crossing_{method}.jpg"
        image.save(path, quality=86, optimize=True)
        seam_row.append(image)
    captures.append(seam_row)

    sheet = Image.new("RGB", (WIDTH * 3, (HEIGHT + 30) * len(captures)), (221, 227, 230))
    draw = ImageDraw.Draw(sheet)
    for row_index, row in enumerate(captures):
        for col_index, image in enumerate(row):
            x, y = col_index * WIDTH, row_index * (HEIGHT + 30)
            sheet.paste(image, (x, y))
            draw.rectangle((x, y, x + WIDTH - 1, y + HEIGHT + 29), outline=(73, 83, 89), width=1)
    sheet_path = CAPTURES / "contact_sheet.jpg"
    sheet.save(sheet_path, quality=86, optimize=True)
    oversized = []
    for path in sorted(CAPTURES.glob("*.jpg")):
        size = path.stat().st_size
        print(f"capture file={path.name} bytes={size}")
        if size >= 300_000:
            oversized.append((path.name, size))
    if oversized:
        raise RuntimeError(f"captures exceed 300 KB: {oversized}")
    print(f"contact_sheet={sheet_path} bytes={sheet_path.stat().st_size} image_count={len(list(CAPTURES.glob('*.jpg')))}")


if __name__ == "__main__":
    main()