"""Check an exported stand with Blender's independent triangle BVH.

blender --background --python-exit-code 1 --python script/inspect_stand.py -- stand.3mf report.json
"""
import json
import sys
import zipfile
import xml.etree.ElementTree as ET

from mathutils import Vector
from mathutils.bvhtree import BVHTree

source, output = sys.argv[sys.argv.index("--") + 1:]
with zipfile.ZipFile(source) as package:
    root = ET.fromstring(package.read("3D/3dmodel.model"))
ns = {"m": "http://schemas.microsoft.com/3dmanufacturing/core/2015/02"}
meshes = root.findall("m:resources/m:object/m:mesh", ns)
assert len(meshes) == 1, "Stand must contain exactly one mesh"
mesh = meshes[0]
vertices = [Vector(tuple(float(v.attrib[a]) for a in ("x", "y", "z")))
            for v in mesh.findall("m:vertices/m:vertex", ns)]
faces = [tuple(int(t.attrib[a]) for a in ("v1", "v2", "v3"))
         for t in mesh.findall("m:triangles/m:triangle", ns)]
edges = {}
volume = 0.0
flat_area = 0.0
minimum_area = float("inf")
for i, face in enumerate(faces):
    a, b, c = (vertices[j] for j in face)
    normal = (b-a).cross(c-a)
    area = normal.length / 2
    minimum_area = min(minimum_area, area)
    volume += a.dot(normal) / 6
    if max(abs(p.z) for p in (a, b, c)) < 1e-6:
        flat_area += area
    for j in range(3):
        u, v = face[j], face[(j+1) % 3]
        edges.setdefault(tuple(sorted((u, v))), []).append((i, u < v))
assert all(len(uses) == 2 and uses[0][1] != uses[1][1] for uses in edges.values()), "Open edges or inconsistent winding"
adjacency = [[] for _ in faces]
for uses in edges.values():
    a, b = uses[0][0], uses[1][0]
    adjacency[a].append(b)
    adjacency[b].append(a)
seen = {0}
queue = [0]
for face in queue:
    for other in adjacency[face]:
        if other not in seen:
            seen.add(other)
            queue.append(other)
assert len(seen) == len(faces), "Disconnected shells"
assert volume > 0 and minimum_area > 1e-10 and flat_area > 0
tree = BVHTree.FromPolygons(vertices, faces, all_triangles=True, epsilon=0.0)
intersections = []
for a, b in tree.overlap(tree):
    if a >= b or set(faces[a]).intersection(faces[b]):
        continue
    intersections.append((a, b))
report = {"vertices": len(vertices), "triangles": len(faces), "connected_components": 1,
          "watertight_consistent_winding": True, "signed_volume_mm3": volume,
          "table_contact_area_mm2": flat_area, "minimum_triangle_area_mm2": minimum_area,
          "nonadjacent_triangle_intersections": len(intersections), "intersection_examples": intersections[:10],
          "slicer_validation": "This script checks the mesh only; consult separate slicer results when available."}
with open(output, "w") as stream:
    json.dump(report, stream, indent=2)
print(json.dumps(report, indent=2))
assert not intersections, "Nonadjacent triangles intersect"