#!/usr/bin/env python3
# $Id: StlExtrude 181 2026-07-13 10:56:05Z peter $
#
# StlExtrude - extrude a potrace GeoJSON silhouette to a millimetre STL.
#
# Copyright (C) 2026  Peter Felecan
# Copyright (C) 2026  Anthropic AI (Claude Sonnet 5)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
#
# DESCRIPTION:
#   Reads a potrace GeoJSON silhouette (one or more polygons, each an
#   exterior ring plus zero or more hole rings), scales it to millimetres,
#   and writes a binary STL of the silhouette extruded to a flat plate.
#   The cap faces are triangulated with the vendored earcut module; the
#   side walls and the two caps form a closed (watertight) solid. This is
#   the default backend for dvi3v's --stl option and needs only Python 3,
#   so an end user does not have to install OpenSCAD.
#
# USAGE:
#   StlExtrude INPUT_GEOJSON OUTPUT_STL MM_PER_PX THICKNESS_MM
import json
import os
import struct
import sys

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import earcut as ec  # vendored; see earcut.py


def ring_area(lRing):
    lSum = 0.0
    lN = len(lRing)
    for iC in range(lN):
        lX1, lY1 = lRing[iC]
        lX2, lY2 = lRing[(iC + 1) % lN]
        lSum += lX1 * lY2 - lX2 * lY1
    return lSum / 2.0


def as_ccw(lRing):
    return lRing if ring_area(lRing) > 0 else lRing[::-1]


def as_cw(lRing):
    return lRing if ring_area(lRing) < 0 else lRing[::-1]


def load_polygons(sPath):
    lData = json.load(open(sPath))
    lFeatures = lData.get("features", [lData])
    lPolygons = []
    for lFeature in lFeatures:
        lGeom = lFeature.get("geometry", lFeature)
        sType = lGeom["type"]
        if sType == "Polygon":
            lPolygons.append(lGeom["coordinates"])
        elif sType == "MultiPolygon":
            lPolygons.extend(lGeom["coordinates"])
    return lPolygons


def main():
    if len(sys.argv) != 5:
        sys.stderr.write(
            "Usage: StlExtrude INPUT_GEOJSON OUTPUT_STL MM_PER_PX THICKNESS_MM\n")
        return 2
    sIn, sOut = sys.argv[1], sys.argv[2]
    try:
        fMmPerPx = float(sys.argv[3])
        fThick = float(sys.argv[4])
    except ValueError:
        sys.stderr.write("Error: MM_PER_PX and THICKNESS_MM must be numbers.\n")
        return 2
    if fThick <= 0:
        sys.stderr.write("Error: THICKNESS_MM must be positive.\n")
        return 2

    lTriangles = []

    def emit(lA, lB, lC):
        lTriangles.append((lA, lB, lC))

    for lRings in load_polygons(sIn):
        # Scale every ring to millimetres and drop a repeated closing vertex.
        lScaled = []
        for lRing in lRings:
            lR = [(lX * fMmPerPx, lY * fMmPerPx) for lX, lY in lRing]
            if len(lR) > 1 and lR[0] == lR[-1]:
                lR = lR[:-1]
            lScaled.append(lR)
        if not lScaled:
            continue
        lOuter = as_ccw(lScaled[0])
        lHoles = [as_cw(lH) for lH in lScaled[1:]]

        # Triangulate the cap (polygon with holes) with earcut.
        lFlat = ec.flatten([lOuter] + lHoles)
        lIndices = ec.earcut(lFlat["vertices"], lFlat["holes"],
                             lFlat["dimensions"])
        lVerts = lFlat["vertices"]

        def vert(iIndex):
            return (lVerts[2 * iIndex], lVerts[2 * iIndex + 1])

        # Top cap at z = thickness (CCW, normal +Z) and bottom at z = 0
        # (reversed winding, normal -Z).
        for iC in range(0, len(lIndices), 3):
            lA = vert(lIndices[iC])
            lB = vert(lIndices[iC + 1])
            lC = vert(lIndices[iC + 2])
            emit((lA[0], lA[1], fThick), (lB[0], lB[1], fThick),
                 (lC[0], lC[1], fThick))
            emit((lA[0], lA[1], 0.0), (lC[0], lC[1], 0.0),
                 (lB[0], lB[1], 0.0))

        # Side walls for every ring (exterior CCW, holes CW), two triangles
        # per edge, connecting the bottom loop to the top loop.
        for lRing in [lOuter] + lHoles:
            lN = len(lRing)
            for iC in range(lN):
                lX1, lY1 = lRing[iC]
                lX2, lY2 = lRing[(iC + 1) % lN]
                emit((lX1, lY1, 0.0), (lX2, lY2, 0.0), (lX2, lY2, fThick))
                emit((lX1, lY1, 0.0), (lX2, lY2, fThick), (lX1, lY1, fThick))

    if not lTriangles:
        sys.stderr.write("Error: no geometry produced from the GeoJSON.\n")
        return 1

    with open(sOut, "wb") as oFile:
        oFile.write(b"\0" * 80)
        oFile.write(struct.pack("<I", len(lTriangles)))
        for lA, lB, lC in lTriangles:
            oFile.write(struct.pack("<3f", 0.0, 0.0, 0.0))  # normal (recomputed)
            for lV in (lA, lB, lC):
                oFile.write(struct.pack("<3f", lV[0], lV[1], lV[2]))
            oFile.write(struct.pack("<H", 0))
    return 0


if __name__ == "__main__":
    sys.exit(main())
