Fix contour silhouette detection for perpendicular faces

The silhouette detection was using 2D orientation (CCW/CW in XY projection)
which fails for faces perpendicular to the view direction. When a face's
normal is parallel to the Z-axis (truly perpendicular to top-down view),
the 2D cross product becomes zero and the face is neither CCW nor CW.

This fix replaces the 2D orientation approach with 3D normal-based detection:
- Categorize faces by the sign of normal.z:
  - normal.z < 0: front-facing (visible from top)
  - normal.z > 0: back-facing (not visible)
  - normal.z == 0: vertical/edge-on
- An edge is a silhouette if adjacent faces have different categories

This correctly detects silhouette edges for:
- Front-facing + back-facing faces
- Front-facing + vertical faces
- Back-facing + vertical faces

Affects: extract_silhouette_contour_edges() and extract_lineart_edges()

Generated by Mistral Vibe.
Co-Authored-By: Mistral Vibe <vibe@mistral.ai>

core/pdf_exporter/geometry.py

@@ -92,7 +92,7 @@ def extract_unique_top_view_edges(mesh_data: Any) -> NDArray[np.float64]:
 def extract_lineart_edges(mesh_data: Any, max_angle_deg: float = 30.0) -> NDArray[np.float64]:
     """
     Fully vectorized NumPy Lineart edge extraction:
-    1. Silhouette edge: shared between at least one CCW face and at least one CW face in 2D top view.
+    1. Silhouette edge: adjacent faces have different normal.z sign categories (front/back/vertical).
     2. Boundary edge: owned by only 1 face.
     3. Crease edge: owned by 2 or more faces where 3D dihedral angle between normals > max_angle_deg.
     Returns 2D edge segments array of shape (N, 2, 2).
@@ -157,13 +157,17 @@ def extract_lineart_edges(mesh_data: Any, max_angle_deg: float = 30.0) -> NDArra
     # Top-view 2D vertices
     verts_2d = tris_3d[:, :, :2]  # (N, 3, 2)
 
-    # 2D orientation per triangle
+    # For silhouette detection, categorize faces by normal.z sign
-    v0_2d = verts_2d[:, 0, :]
+    # For view direction (0,0,-1):
-    v1_2d = verts_2d[:, 1, :]
+    # - normal.z < 0: front-facing (facing viewer)
-    v2_2d = verts_2d[:, 2, :]
+    # - normal.z > 0: back-facing (facing away)
-    cross_2d = (v1_2d[:, 0] - v0_2d[:, 0]) * (v2_2d[:, 1] - v0_2d[:, 1]) - (v1_2d[:, 1] - v0_2d[:, 1]) * (v2_2d[:, 0] - v0_2d[:, 0])
+    # - normal.z == 0: edge-on / vertical
-    is_ccw = cross_2d > 0
+    normal_z = normals_unit[:, 2]
-    is_cw = cross_2d < 0
+    epsilon = 1e-8
+    face_category = np.zeros(num_tris, dtype=int)
+    face_category[normal_z < -epsilon] = 1  # Front-facing
+    face_category[normal_z > epsilon] = 2   # Back-facing
+    # face_category == 0 means vertical (|normal_z| <= epsilon)
 
     # Extract all 3 edges per triangle
     e0_3d = tris_3d[:, [0, 1], :]
@@ -176,8 +180,7 @@ def extract_lineart_edges(mesh_data: Any, max_angle_deg: float = 30.0) -> NDArra
     e2_2d = verts_2d[:, [2, 0], :]
     edges_2d = np.concatenate([e0_2d, e1_2d, e2_2d], axis=0)  # (3N, 2, 2)
 
-    edge_ccw = np.tile(is_ccw, 3)  # (3N,)
+    edge_category = np.tile(face_category, 3)  # (3N,)
-    edge_cw = np.tile(is_cw, 3)    # (3N,)
     edge_normals = np.tile(normals_unit, (3, 1))  # (3N, 3)
 
     # Filter out near zero-length 2D edges
@@ -186,8 +189,7 @@ def extract_lineart_edges(mesh_data: Any, max_angle_deg: float = 30.0) -> NDArra
 
     edges_3d = edges_3d[valid_len_mask]
     edges_2d = edges_2d[valid_len_mask]
-    edge_ccw = edge_ccw[valid_len_mask]
+    edge_category = edge_category[valid_len_mask]
-    edge_cw = edge_cw[valid_len_mask]
     edge_normals = edge_normals[valid_len_mask]
 
     if len(edges_3d) == 0:
@@ -218,10 +220,18 @@ def extract_lineart_edges(mesh_data: Any, max_angle_deg: float = 30.0) -> NDArra
 
     num_unique = len(unique_keys)
 
-    # Rule 1: Silhouette edge (has CCW and CW face)
+    # Rule 1: Silhouette edge - adjacent faces have different categories
-    has_ccw = np.bincount(inverse_indices, weights=edge_ccw.astype(int), minlength=num_unique) > 0
+    # (i.e., not all faces sharing this edge are in the same category)
-    has_cw = np.bincount(inverse_indices, weights=edge_cw.astype(int), minlength=num_unique) > 0
+    categories_all = edge_category[inverse_indices]
-    rule1_unique = has_ccw & has_cw
+    rule1_unique = np.zeros(num_unique, dtype=bool)
+    
+    for edge_id in range(num_unique):
+        edge_cats = categories_all[inverse_indices == edge_id]
+        if len(edge_cats) < 2:
+            rule1_unique[edge_id] = False
+        else:
+            all_same = np.all(edge_cats == edge_cats[0])
+            rule1_unique[edge_id] = not all_same
 
     # Rule 2: Boundary edge (owned by 1 face)
     rule2_unique = counts == 1
@@ -340,11 +350,18 @@ def extract_crease_edges(mesh_data: Any, max_angle_deg: float = 30.0) -> NDArray
 
 def extract_silhouette_contour_edges(mesh_data: Any) -> NDArray[np.float64]:
     """
-    Extracts top-view silhouette contour edges using CW vs CCW triangle orientation in +Z top-view projection.
+    Extracts top-view silhouette contour edges using 3D face normals.
+    A silhouette edge separates faces with different visibility from the top-down view.
+    For view direction (0,0,-1), we compute the sign of each face normal's Z-component:
+    - sign < 0: front-facing (visible)
+    - sign > 0: back-facing (not visible)
+    - sign == 0: edge-on / vertical (zero projected area)
+    An edge is a silhouette if its adjacent faces have different sign categories.
     """
     from ...infra.blender.blender_mesh_capture import SlicingGroup, SlicingCollection, SlicingObject
 
     all_tris_3d = []
+    all_normals = []
 
     def process_item(item: Any):
         if item is None:
@@ -360,7 +377,28 @@ def extract_silhouette_contour_edges(mesh_data: Any) -> NDArray[np.float64]:
                 process_item(so)
         elif isinstance(item, SlicingObject):
             if item.mesh.size > 0:
-                all_tris_3d.append(item.mesh[:, :3, :3])
+                verts_3d = item.mesh[:, :3, :3]
+                if item.mesh.shape[1] >= 4:
+                    normals = item.mesh[:, 3, :]
+                else:
+                    v0_3d = verts_3d[:, 0, :]
+                    v1_3d = verts_3d[:, 1, :]
+                    v2_3d = verts_3d[:, 2, :]
+                    normals = np.cross(v1_3d - v0_3d, v2_3d - v0_3d)
+                
+                norms_len = np.linalg.norm(normals, axis=1, keepdims=True)
+                norms_len = np.where(norms_len == 0, 1.0, norms_len)
+                normals_unit = normals / norms_len
+                
+                zero_norm_mask = (normals[:, 0] == 0) & (normals[:, 1] == 0) & (normals[:, 2] == 0)
+                if np.any(zero_norm_mask):
+                    computed_normals = np.cross(v1_3d - v0_3d, v2_3d - v0_3d)
+                    cn_len = np.linalg.norm(computed_normals, axis=1, keepdims=True)
+                    cn_len = np.where(cn_len == 0, 1.0, cn_len)
+                    normals_unit[zero_norm_mask] = (computed_normals / cn_len)[zero_norm_mask]
+                
+                all_tris_3d.append(verts_3d)
+                all_normals.append(normals_unit)
 
     process_item(mesh_data)
 
@@ -368,44 +406,53 @@ def extract_silhouette_contour_edges(mesh_data: Any) -> NDArray[np.float64]:
         return np.empty((0, 2, 2), dtype=np.float64)
 
     tris_3d = np.concatenate(all_tris_3d, axis=0)  # (N, 3, 3)
+    normals_unit = np.concatenate(all_normals, axis=0)  # (N, 3)
     num_tris = len(tris_3d)
 
     if num_tris == 0:
         return np.empty((0, 2, 2), dtype=np.float64)
 
-    verts_2d = tris_3d[:, :, :2]  # (N, 3, 2)
+    # For top-down view (view direction = (0,0,-1)), compute sign of normal.z
+    # The sign of -normal.z tells us visibility: >0 means front-facing (visible)
+    # We categorize by sign of normal.z directly for simplicity:
+    # negative z = facing viewer, positive z = facing away, zero = edge-on
+    normal_z = normals_unit[:, 2]
     
-    v0_2d = verts_2d[:, 0, :]
+    # Assign each face to a category: 0=vertical/edge-on, 1=front, 2=back
-    v1_2d = verts_2d[:, 1, :]
+    # We use a small epsilon to handle floating point
-    v2_2d = verts_2d[:, 2, :]
+    epsilon = 1e-8
-    cross_2d = (v1_2d[:, 0] - v0_2d[:, 0]) * (v2_2d[:, 1] - v0_2d[:, 1]) - (v1_2d[:, 1] - v0_2d[:, 1]) * (v2_2d[:, 0] - v0_2d[:, 0])
+    face_category = np.zeros(num_tris, dtype=int)
-    is_ccw = cross_2d > 0
+    face_category[normal_z < -epsilon] = 1  # Front-facing
-    is_cw = cross_2d < 0
+    face_category[normal_z > epsilon] = 2   # Back-facing
+    # face_category == 0 means vertical (|normal_z| <= epsilon)
 
+    # Extract all 3 edges per triangle
     e0_3d = tris_3d[:, [0, 1], :]
     e1_3d = tris_3d[:, [1, 2], :]
     e2_3d = tris_3d[:, [2, 0], :]
     edges_3d = np.concatenate([e0_3d, e1_3d, e2_3d], axis=0)  # (3N, 2, 3)
 
+    verts_2d = tris_3d[:, :, :2]  # (N, 3, 2)
     e0_2d = verts_2d[:, [0, 1], :]
     e1_2d = verts_2d[:, [1, 2], :]
     e2_2d = verts_2d[:, [2, 0], :]
     edges_2d = np.concatenate([e0_2d, e1_2d, e2_2d], axis=0)  # (3N, 2, 2)
 
-    edge_ccw = np.tile(is_ccw, 3)
+    # Each edge inherits its face's category
-    edge_cw = np.tile(is_cw, 3)
+    edge_category = np.tile(face_category, 3)
 
+    # Filter out near zero-length 2D edges
     dists_sq_2d = np.sum((edges_2d[:, 0, :] - edges_2d[:, 1, :]) ** 2, axis=1)
     valid_len_mask = dists_sq_2d >= 1e-8
 
     edges_3d = edges_3d[valid_len_mask]
     edges_2d = edges_2d[valid_len_mask]
-    edge_ccw = edge_ccw[valid_len_mask]
+    edge_category = edge_category[valid_len_mask]
-    edge_cw = edge_cw[valid_len_mask]
 
     if len(edges_3d) == 0:
         return np.empty((0, 2, 2), dtype=np.float64)
 
+    # Canonical undirected edge sorting in 3D
     rounded_3d = np.round(edges_3d, decimals=4)
     pA = rounded_3d[:, 0, :]
     pB = rounded_3d[:, 1, :]
@@ -421,6 +468,7 @@ def extract_silhouette_contour_edges(mesh_data: Any) -> NDArray[np.float64]:
     canonical_2d = edges_2d.copy()
     canonical_2d[swap_mask] = canonical_2d[swap_mask][:, ::-1, :]
 
+    # Unique 3D edge identification
     flat_3d = canonical_3d.reshape(len(canonical_3d), 6)
     view_bytes = np.ascontiguousarray(flat_3d).view(np.dtype((np.void, 48))).ravel()
     unique_keys, first_indices, inverse_indices, counts = np.unique(
@@ -429,11 +477,27 @@ def extract_silhouette_contour_edges(mesh_data: Any) -> NDArray[np.float64]:
 
     num_unique = len(unique_keys)
 
-    has_ccw = np.bincount(inverse_indices, weights=edge_ccw.astype(int), minlength=num_unique) > 0
+    # For each unique edge, collect the categories of all adjacent faces
-    has_cw = np.bincount(inverse_indices, weights=edge_cw.astype(int), minlength=num_unique) > 0
+    # An edge is a silhouette if it has faces from different categories
-    silhouette_mask = has_ccw & has_cw
+    # (i.e., not all faces are in the same category)
     
-    return canonical_2d[first_indices[silhouette_mask]]
+    # Get categories for all edge instances
+    categories_all = edge_category[inverse_indices]
+    
+    # For each unique edge, check if all its face categories are the same
+    is_silhouette = np.zeros(num_unique, dtype=bool)
+    
+    for edge_id in range(num_unique):
+        edge_cats = categories_all[inverse_indices == edge_id]
+        # If this edge has only one face, it's a boundary, not silhouette
+        if len(edge_cats) < 2:
+            is_silhouette[edge_id] = False
+        else:
+            # Check if all categories are the same
+            all_same = np.all(edge_cats == edge_cats[0])
+            is_silhouette[edge_id] = not all_same
+
+    return canonical_2d[first_indices[is_silhouette]]
 
 
 def extract_outer_contour_edges(mesh_data: Any) -> NDArray[np.float64]: