pdm/creoparams/mesh.py

"""Lettura di mesh triangolari (STL) e calcolo delle grandezze derivate.

La geometria dei `.prt` e' scritta nel formato del kernel Granite di PTC e
non e' leggibile senza la libreria di PTC: verificato cercando nel file i 22
valori esatti del modello (raggi e coordinate presi da un export STEP) in
cinque codifiche diverse, senza un solo riscontro.

La forma del pezzo arriva quindi da una mesh esportata da Creo. Questo
modulo la legge; la fonte e' sostituibile senza toccare il resto (per
esempio con il Granite Interoperability Kernel, che legge i .prt nativi
senza Creo installato).
"""

from __future__ import annotations

import math
import struct
from dataclasses import dataclass

BINARY_HEADER = 84
TRIANGLE_SIZE = 50


@dataclass
class Mesh:
    """Triangoli come terne di vertici (x, y, z)."""

    triangles: list
    source: str = ""

    def __len__(self):
        return len(self.triangles)

    @property
    def bounds(self):
        xs = [v[0] for t in self.triangles for v in t]
        ys = [v[1] for t in self.triangles for v in t]
        zs = [v[2] for t in self.triangles for v in t]
        return (min(xs), min(ys), min(zs)), (max(xs), max(ys), max(zs))

    @property
    def size(self):
        lo, hi = self.bounds
        return tuple(round(hi[i] - lo[i], 4) for i in range(3))

    @property
    def center(self):
        lo, hi = self.bounds
        return tuple((hi[i] + lo[i]) / 2 for i in range(3))

    def area(self) -> float:
        """Area della superficie, in unita' del modello al quadrato."""
        total = 0.0
        for a, b, c in self.triangles:
            ux, uy, uz = b[0]-a[0], b[1]-a[1], b[2]-a[2]
            vx, vy, vz = c[0]-a[0], c[1]-a[1], c[2]-a[2]
            nx, ny, nz = uy*vz-uz*vy, uz*vx-ux*vz, ux*vy-uy*vx
            total += math.sqrt(nx*nx + ny*ny + nz*nz) / 2
        return total

    def volume(self) -> float:
        """Volume col metodo dei tetraedri con segno. Vale solo se la mesh
        e' chiusa: su una mesh aperta il risultato non ha senso."""
        total = 0.0
        for a, b, c in self.triangles:
            total += (a[0]*(b[1]*c[2]-b[2]*c[1])
                      - a[1]*(b[0]*c[2]-b[2]*c[0])
                      + a[2]*(b[0]*c[1]-b[1]*c[0])) / 6
        return abs(total)

    def indexed(self, tolerance: int = 5):
        """Vertici unici, facce e spigoli, con le facce adiacenti a ciascuno.

        Un STL ripete i vertici a ogni triangolo: fondendoli possiamo sapere
        quali facce condividono uno spigolo, e quindi distinguere gli spigoli
        vivi del pezzo dalle suddivisioni interne della tassellatura. Sono
        quelli che rendono leggibile un modello meccanico.
        """
        lookup, vertices, faces = {}, [], []
        for triangle in self.triangles:
            face = []
            for vertex in triangle:
                key = tuple(round(c, tolerance) for c in vertex)
                index = lookup.get(key)
                if index is None:
                    index = len(vertices)
                    lookup[key] = index
                    vertices.append([round(c, 4) for c in vertex])
                face.append(index)
            if len(set(face)) == 3:          # scarta i triangoli degeneri
                faces.append(face)

        edges = {}
        for number, (a, b, c) in enumerate(faces):
            for pair in ((a, b), (b, c), (c, a)):
                key = (min(pair), max(pair))
                edges.setdefault(key, []).append(number)

        edge_list = [[a, b, adj[0], adj[1] if len(adj) > 1 else -1]
                     for (a, b), adj in edges.items()]
        return vertices, faces, edge_list

    def to_dict(self) -> dict:
        lo, hi = self.bounds
        return {
            "source": self.source,
            "triangles": len(self.triangles),
            "size": self.size,
            "bounds": {"min": [round(v, 4) for v in lo],
                       "max": [round(v, 4) for v in hi]},
            "area": round(self.area(), 3),
            "volume": round(self.volume(), 3),
        }


def _is_binary(data: bytes) -> bool:
    """Un STL ASCII inizia con 'solid', ma anche molti binari lo fanno:
    l'unico controllo affidabile e' la dimensione attesa."""
    if len(data) < BINARY_HEADER:
        return False
    count = struct.unpack("<I", data[80:84])[0]
    return len(data) == BINARY_HEADER + count * TRIANGLE_SIZE


def _read_binary(data: bytes) -> list:
    count = struct.unpack("<I", data[80:84])[0]
    triangles = []
    offset = BINARY_HEADER
    for _ in range(count):
        values = struct.unpack("<12f", data[offset:offset + 48])
        triangles.append((values[3:6], values[6:9], values[9:12]))
        offset += TRIANGLE_SIZE
    return triangles


def _read_ascii(data: bytes) -> list:
    triangles, current = [], []
    for line in data.decode("ascii", errors="replace").splitlines():
        parts = line.split()
        if len(parts) == 4 and parts[0] == "vertex":
            try:
                current.append(tuple(float(v) for v in parts[1:]))
            except ValueError:
                current = []
            if len(current) == 3:
                triangles.append(tuple(current))
                current = []
        elif parts and parts[0] == "facet":
            current = []
    return triangles


def read_stl(path: str) -> Mesh:
    with open(path, "rb") as handle:
        data = handle.read()
    triangles = _read_binary(data) if _is_binary(data) else _read_ascii(data)
    if not triangles:
        raise ValueError("nessun triangolo trovato: il file non e' un STL valido")
    return Mesh(triangles, source=path)