"""
2D Vector Edge Extraction Module for PDF Exporter.
Extracts top-view wireframe, lineart, boundary, contour, and crease edges from 3D SlicingObjects.
"""
import math
from typing import Any
import numpy as np
from numpy.typing import NDArray

def deduplicate_2d_edges(edges_2d: NDArray[np.float64]) -> NDArray[np.float64]:
    """
    Deduplicates 2D edge segments regardless of direction (A->B or B->A).
    """
    if len(edges_2d) == 0:
        return edges_2d

    p0 = edges_2d[:, 0, :]
    p1 = edges_2d[:, 1, :]

    swap_mask = (p0[:, 0] > p1[:, 0]) | ((p0[:, 0] == p1[:, 0]) & (p0[:, 1] > p1[:, 1]))

    canonical_p0 = np.where(swap_mask[:, None], p1, p0)
    canonical_p1 = np.where(swap_mask[:, None], p0, p1)

    quantized_p0 = np.round(canonical_p0, decimals=4)
    quantized_p1 = np.round(canonical_p1, decimals=4)

    stacked = np.column_stack([quantized_p0, quantized_p1])

    _, unique_indices = np.unique(stacked, axis=0, return_index=True)
    return edges_2d[unique_indices]


def transform_slice_paths(
    sp: Any,
    t_arr: NDArray[np.float64],
    r_arr: NDArray[np.float64],
    s_arr: NDArray[np.float64],
) -> Any:
    """
    Transforms 3D coordinates in a SlicePaths object using translation, Euler rotation (radians), and scale.
    """
    from ..slice import SlicePaths

    if not isinstance(sp, SlicePaths):
        return sp

    has_scale = any(s != 1.0 for s in s_arr)
    has_rot = any(r != 0.0 for r in r_arr)
    has_trans = any(t != 0.0 for t in t_arr)

    if not (has_scale or has_rot or has_trans):
        return sp.copy_with(paths=list(sp))

    R = None
    if has_rot:
        rx, ry, rz = r_arr
        cx, sx = np.cos(rx), np.sin(rx)
        cy, sy = np.cos(ry), np.sin(ry)
        cz, sz = np.cos(rz), np.sin(rz)

        Rx = np.array([
            [1.0, 0.0, 0.0],
            [0.0, cx, -sx],
            [0.0, sx, cx]
        ], dtype=np.float64)
        Ry = np.array([
            [cy, 0.0, sy],
            [0.0, 1.0, 0.0],
            [-sy, 0.0, cy]
        ], dtype=np.float64)
        Rz = np.array([
            [cz, -sz, 0.0],
            [sz, cz, 0.0],
            [0.0, 0.0, 1.0]
        ], dtype=np.float64)
        R = Rz @ (Ry @ Rx)

    new_paths = []
    for p in sp:
        if len(p) == 0:
            new_paths.append(p.copy())
            continue

        new_p = p.astype(np.float64, copy=True)
        dim = new_p.shape[1] if new_p.ndim > 1 else 0

        if has_scale and dim > 0:
            if dim >= 3:
                new_p[:, :3] *= s_arr
            else:
                new_p[:, :dim] *= s_arr[:dim]

        if has_rot and R is not None and dim > 0:
            if dim >= 3:
                new_p[:, :3] = new_p[:, :3] @ R.T
            elif dim == 2:
                rot_2d = R[:2, :2]
                new_p[:, :2] = new_p[:, :2] @ rot_2d.T

        if has_trans and dim > 0:
            if dim >= 3:
                new_p[:, :3] += t_arr
            else:
                new_p[:, :dim] += t_arr[:dim]

        new_paths.append(new_p)

    return sp.copy_with(paths=new_paths)


def get_first_object_rotation(
    slicing_objs: list[Any],
    raw_geometry: Any = None,
    slice_paths: Any = None,
) -> tuple[float, float, float]:
    """
    Extracts the accumulated 3D Euler rotation angles (rx, ry, rz) in radians
    from the first object or raw_geometry / slice_paths.
    """
    candidates: list[Any] = []
    if slicing_objs:
        candidates.extend(slicing_objs)
    if raw_geometry is not None and raw_geometry not in candidates:
        candidates.append(raw_geometry)
    if slice_paths is not None and slice_paths not in candidates:
        candidates.append(slice_paths)

    for item in candidates:
        transforms = getattr(item, "transforms", None)
        if transforms:
            rx, ry, rz = 0.0, 0.0, 0.0
            has_rot = False
            for t in transforms:
                if isinstance(t, dict) and "rotation" in t:
                    r = t["rotation"]
                    rx += float(r[0])
                    ry += float(r[1])
                    rz += float(r[2])
                    has_rot = True
            if has_rot:
                return (rx, ry, rz)

    return (0.0, 0.0, 0.0)


def map_2d_to_3d(u: float, v: float, d: float, plane: str = "XY") -> NDArray[np.float64]:
    p = plane.upper()
    if p == "XZ":
        return np.array([u, d, v], dtype=np.float64)
    elif p == "YZ":
        return np.array([d, u, v], dtype=np.float64)
    else:  # XY
        return np.array([u, v, d], dtype=np.float64)


def get_plane_2d_indices(plane: str = "XY") -> tuple[list[int], int]:
    p = plane.upper()
    if p == "XZ":
        return [0, 2], 1
    elif p == "YZ":
        return [1, 2], 0
    else:  # XY
        return [0, 1], 2


def generate_grid_slice_paths(
    grid_min_x: float,
    grid_max_x: float,
    grid_min_y: float,
    grid_max_y: float,
    minor_step: float,
    max_grid_steps: int = 2000,
    z_level: float = 0.0,
    plane: str = "XY",
) -> tuple[Any, Any]:
    """
    Generates 3D grid line geometry sitting at local depth d=z_level on specified local plane (XY, XZ, YZ)
    as two SlicePaths objects: (minor_grid_slice_paths, main_grid_slice_paths).
    """
    from ..slice import SlicePaths

    k_x_start = int(np.round(grid_min_x / minor_step))
    k_x_end = int(np.round(grid_max_x / minor_step))
    k_x = np.arange(k_x_start, k_x_end + 1)

    k_y_start = int(np.round(grid_min_y / minor_step))
    k_y_end = int(np.round(grid_max_y / minor_step))
    k_y = np.arange(k_y_start, k_y_end + 1)

    if len(k_x) > max_grid_steps or len(k_y) > max_grid_steps:
        return SlicePaths(), SlicePaths()

    minor_paths: list[NDArray[np.float64]] = []
    main_paths: list[NDArray[np.float64]] = []

    # Vertical grid lines (constant U)
    for k in k_x:
        gu = float(k * minor_step)
        p1 = map_2d_to_3d(gu, grid_min_y, z_level, plane)
        p2 = map_2d_to_3d(gu, grid_max_y, z_level, plane)
        line = np.array([p1, p2], dtype=np.float64)
        if k % 10 == 0:
            main_paths.append(line)
        else:
            minor_paths.append(line)

    # Horizontal grid lines (constant V)
    for k in k_y:
        gv = float(k * minor_step)
        p1 = map_2d_to_3d(grid_min_x, gv, z_level, plane)
        p2 = map_2d_to_3d(grid_max_x, gv, z_level, plane)
        line = np.array([p1, p2], dtype=np.float64)
        if k % 10 == 0:
            main_paths.append(line)
        else:
            minor_paths.append(line)

    return SlicePaths(minor_paths), SlicePaths(main_paths)


def generate_legend_slice_paths(
    grid_min_x: float,
    grid_min_y: float,
    grid_max_x: float,
    main_step: float,
    minor_step: float,
    line_width_pt: float,
    font_size_pt: float,
    z_level: float = 0.0,
    plane: str = "XY",
    draw_x: bool = True,
    draw_y: bool = True,
    draw_z: bool = False,
    local_origin_3d: tuple[float, float, float] | None = None,
    local_max_3d: tuple[float, float, float] | None = None,
) -> tuple[Any, list[dict[str, Any]]]:
    """
    Generates 3D legend line geometry (axes, arrowheads, tick marks) starting at origin (min x, max y, min z)
    as a SlicePaths object, alongside tick and label metadata for text positioning.
    Returns: (legend_slice_paths, label_metadata_list)
    """
    from ..slice import SlicePaths
    from .sizing import LegendSizing, MM_TO_PT

    leg_sizing = LegendSizing(main_step, minor_step, line_width_pt, font_size_pt)

    if local_origin_3d is not None:
        x0, y0, z0 = local_origin_3d
    else:
        x0 = grid_min_x
        y0 = grid_min_y
        z0 = z_level

    max_x_val = local_max_3d[0] if local_max_3d else grid_max_x
    min_y_val = local_max_3d[1] if local_max_3d else grid_min_y
    max_z_val = local_max_3d[2] if local_max_3d else z_level

    arrow_len_mm = leg_sizing.arrow_len_pt / MM_TO_PT
    arrow_half_w_mm = leg_sizing.arrow_half_w_pt / MM_TO_PT
    label_margin_mm = leg_sizing.label_margin_h_pt / MM_TO_PT

    legend_paths: list[NDArray[np.float64]] = []
    label_meta: list[dict[str, Any]] = []

    # 1. Local X Axis (extends in +X direction from x0)
    if draw_x:
        x_grid_end_mm = min(max_x_val, x0 + 10.0 * main_step) if max_x_val > x0 else x0 + 10.0 * main_step
        num_main_steps = max(1, int(np.round(abs(x_grid_end_mm - x0) / main_step)))

        x_arrow_base = x_grid_end_mm + leg_sizing.x_axis_extension_mm
        x_arrow_tip = x_arrow_base + arrow_len_mm

        legend_paths.append(np.array([[x0, y0, z0], [x_arrow_tip, y0, z0]], dtype=np.float64))
        legend_paths.append(np.array([[x_arrow_base, y0 + arrow_half_w_mm, z0], [x_arrow_tip, y0, z0]], dtype=np.float64))
        legend_paths.append(np.array([[x_arrow_tip, y0, z0], [x_arrow_base, y0 - arrow_half_w_mm, z0]], dtype=np.float64))

        for k in range(num_main_steps + 1):
            xk = x0 + k * main_step
            tick_h = leg_sizing.major_tick_mm if k in (0, 5, num_main_steps) else leg_sizing.minor_tick_mm
            legend_paths.append(np.array([[xk, y0, z0], [xk, y0 - tick_h, z0]], dtype=np.float64))

            if k > 0:
                val_mm = k * main_step
                unit = " m" if (main_step >= 1000.0 and k == num_main_steps) else (" mm" if (main_step < 1000.0 and k == num_main_steps) else "")
                val_str = f"{val_mm/1000.0:g}{unit}" if main_step >= 1000.0 else f"{val_mm:g}{unit}"
                label_meta.append({
                    "type": "tick",
                    "text": val_str,
                    "local_pt": (xk, y0 - tick_h - label_margin_mm, z0),
                })

        label_meta.append({
            "type": "x_axis",
            "text": "X",
            "local_pt": (x_arrow_tip + label_margin_mm, y0, z0),
        })

    # 2. Local Y Axis (extends in -Y direction from y0 towards min_y_val)
    if draw_y:
        y_grid_end_mm = max(min_y_val, y0 - 10.0 * main_step) if min_y_val < y0 else y0 - 10.0 * main_step
        num_main_steps = max(1, int(np.round(abs(y0 - y_grid_end_mm) / main_step)))

        y_arrow_base = y_grid_end_mm - leg_sizing.y_axis_height_mm
        y_arrow_tip = y_arrow_base - arrow_len_mm

        legend_paths.append(np.array([[x0, y0, z0], [x0, y_arrow_tip, z0]], dtype=np.float64))
        legend_paths.append(np.array([[x0 - arrow_half_w_mm, y_arrow_base, z0], [x0, y_arrow_tip, z0]], dtype=np.float64))
        legend_paths.append(np.array([[x0, y_arrow_tip, z0], [x0 + arrow_half_w_mm, y_arrow_base, z0]], dtype=np.float64))

        if not draw_x:
            for k in range(num_main_steps + 1):
                yk = y0 - k * main_step
                tick_h = leg_sizing.major_tick_mm if k in (0, 5, num_main_steps) else leg_sizing.minor_tick_mm
                legend_paths.append(np.array([[x0, yk, z0], [x0 + tick_h, yk, z0]], dtype=np.float64))

                if k > 0:
                    val_mm = k * main_step
                    unit = " m" if (main_step >= 1000.0 and k == num_main_steps) else (" mm" if (main_step < 1000.0 and k == num_main_steps) else "")
                    val_str = f"{val_mm/1000.0:g}{unit}" if main_step >= 1000.0 else f"{val_mm:g}{unit}"
                    label_meta.append({
                        "type": "tick",
                        "text": val_str,
                        "local_pt": (x0 + tick_h + label_margin_mm, yk, z0),
                    })

        label_meta.append({
            "type": "y_axis",
            "text": "Y",
            "local_pt": (x0, y_arrow_tip - label_margin_mm, z0),
        })

    # 3. Local Z Axis (extends in +Z direction from z0 towards max_z_val)
    if draw_z:
        z_grid_end_mm = min(max_z_val, z0 + 10.0 * main_step) if max_z_val > z0 else z0 + 10.0 * main_step
        num_main_steps = max(1, int(np.round(abs(z_grid_end_mm - z0) / main_step)))

        z_arrow_base = z_grid_end_mm + leg_sizing.y_axis_height_mm
        z_arrow_tip = z_arrow_base + arrow_len_mm

        legend_paths.append(np.array([[x0, y0, z0], [x0, y0, z_arrow_tip]], dtype=np.float64))
        legend_paths.append(np.array([[x0 - arrow_half_w_mm, y0, z_arrow_base], [x0, y0, z_arrow_tip]], dtype=np.float64))
        legend_paths.append(np.array([[x0, y0, z_arrow_tip], [x0 + arrow_half_w_mm, y0, z_arrow_base]], dtype=np.float64))

        if not draw_x:
            for k in range(num_main_steps + 1):
                zk = z0 + k * main_step
                tick_h = leg_sizing.major_tick_mm if k in (0, 5, num_main_steps) else leg_sizing.minor_tick_mm
                legend_paths.append(np.array([[x0, y0, zk], [x0 + tick_h, y0, zk]], dtype=np.float64))

                if k > 0:
                    val_mm = k * main_step
                    unit = " m" if (main_step >= 1000.0 and k == num_main_steps) else (" mm" if (main_step < 1000.0 and k == num_main_steps) else "")
                    val_str = f"{val_mm/1000.0:g}{unit}" if main_step >= 1000.0 else f"{val_mm:g}{unit}"
                    label_meta.append({
                        "type": "tick",
                        "text": val_str,
                        "local_pt": (x0 + tick_h + label_margin_mm, y0, zk),
                    })

        label_meta.append({
            "type": "z_axis",
            "text": "Z",
            "local_pt": (x0, y0, z_arrow_tip + label_margin_mm),
        })

    return SlicePaths(legend_paths), label_meta


def extract_slice_paths_edges(slice_paths: Any) -> NDArray[np.float64]:
    """
    Extracts 2D line segments (shape N x 2 x 2) from SlicePaths or a list/collection of SlicePaths.
    """
    from ..slice import SlicePaths

    if slice_paths is None:
        return np.empty((0, 2, 2), dtype=np.float64)

    all_segments: list[np.ndarray] = []

    def process_item(item: Any):
        if item is None:
            return
        elif isinstance(item, SlicePaths):
            for path in item:
                if path is not None and len(path) >= 2:
                    p_arr = np.asarray(path, dtype=np.float64)
                    pts = p_arr[:, :2]
                    segs = np.stack([pts[:-1], pts[1:]], axis=1)
                    all_segments.append(segs)
        elif isinstance(item, (list, tuple)):
            for sub in item:
                process_item(sub)

    process_item(slice_paths)

    if not all_segments:
        return np.empty((0, 2, 2), dtype=np.float64)

    return np.concatenate(all_segments, axis=0)


def extract_unique_top_view_edges(mesh_data: Any) -> NDArray[np.float64]:
    """
    Extracts 2D top-view unique edges (X, Y projection) from mesh_data (SlicingObjects or SlicePaths).
    """
    from ...infra.blender.blender_mesh_capture import SlicingGroup, SlicingCollection, SlicingObject
    from ..slice import SlicePaths

    all_edge_segments: list[NDArray[np.float64]] = []

    def process_item(item: Any):
        if item is None:
            return
        elif isinstance(item, SlicePaths):
            sp_edges = extract_slice_paths_edges(item)
            if len(sp_edges) > 0:
                all_edge_segments.append(sp_edges)
        elif isinstance(item, (list, tuple)):
            for sub in item:
                process_item(sub)
        elif isinstance(item, SlicingGroup):
            for coll in item.collections.values():
                process_item(coll)
        elif isinstance(item, SlicingCollection):
            all_verts, all_idxs = item.unique_verts
            if len(all_idxs) > 0:
                edges = np.vstack([
                    all_idxs[:, [0, 1]],
                    all_idxs[:, [1, 2]],
                    all_idxs[:, [2, 0]]
                ])
                sorted_edges = np.sort(edges, axis=1)
                unique_edge_indices = np.unique(sorted_edges, axis=0)
                edge_3d = all_verts[unique_edge_indices]
                all_edge_segments.append(edge_3d[:, :, :2])
        elif isinstance(item, SlicingObject):
            if item.mesh.size > 0:
                verts_3d = item.mesh[:, :3, :3]
                verts_2d = verts_3d[:, :, :2]
                e0 = verts_2d[:, [0, 1], :]
                e1 = verts_2d[:, [1, 2], :]
                e2 = verts_2d[:, [2, 0], :]
                edges_2d = np.concatenate([e0, e1, e2], axis=0)
                all_edge_segments.append(edges_2d)

    process_item(mesh_data)

    if not all_edge_segments:
        return np.empty((0, 2, 2), dtype=np.float64)

    combined_edges = np.concatenate(all_edge_segments, axis=0)
    if combined_edges.shape[0] == 0:
        return np.empty((0, 2, 2), dtype=np.float64)

    dists_sq = np.sum((combined_edges[:, 0, :] - combined_edges[:, 1, :]) ** 2, axis=1)
    valid_mask = dists_sq > 1e-8
    valid_edges = combined_edges[valid_mask]
    if valid_edges.shape[0] == 0:
        return np.empty((0, 2, 2), dtype=np.float64)

    return deduplicate_2d_edges(valid_edges)


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: 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).
    """
    from ...infra.blender.blender_mesh_capture import SlicingGroup, SlicingCollection, SlicingObject

    all_tris_3d = []
    all_normals = []

    def process_item(item: Any):
        if item is None:
            return
        elif isinstance(item, (list, tuple)):
            for sub in item:
                process_item(sub)
        elif isinstance(item, SlicingGroup):
            for coll in item.collections.values():
                process_item(coll)
        elif isinstance(item, SlicingCollection):
            for so in item.objects:
                process_item(so)
        elif isinstance(item, SlicingObject):
            if item.mesh.size > 0:
                verts_3d = item.mesh[:, :3, :3]  # (N, 3, 3)
                if item.mesh.shape[1] >= 4:
                    normals = item.mesh[:, 3, :]  # (N, 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):
                    v0_3d = verts_3d[:, 0, :]
                    v1_3d = verts_3d[:, 1, :]
                    v2_3d = verts_3d[:, 2, :]
                    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)

    if not all_tris_3d:
        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)

    # 2D vertices projected on plane
    verts_2d = tris_3d[:, :, :2]  # (N, 3, 2)

    normal_z = normals_unit[:, 2]
    epsilon = 1e-3
    face_category = np.zeros(num_tris, dtype=int)
    face_category[normal_z < -epsilon] = 1  # Front-facing
    face_category[normal_z > epsilon] = 2   # Back-facing

    # 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)

    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_category = np.tile(face_category, 3)  # (3N,)
    edge_normals = np.tile(normals_unit, (3, 1))  # (3N, 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_category = edge_category[valid_len_mask]
    edge_normals = edge_normals[valid_len_mask]

    if len(edges_3d) == 0:
        return np.empty((0, 2, 2), dtype=np.float64)

    rounded_3d = np.round(edges_3d, decimals=4)
    pA = rounded_3d[:, 0, :]
    pB = rounded_3d[:, 1, :]
    swap_mask = (pA[:, 0] > pB[:, 0]) | (
        (pA[:, 0] == pB[:, 0]) & (pA[:, 1] > pB[:, 1])
    ) | (
        (pA[:, 0] == pB[:, 0]) & (pA[:, 1] == pB[:, 1]) & (pA[:, 2] > pB[:, 2])
    )

    canonical_3d = rounded_3d.copy()
    canonical_3d[swap_mask] = canonical_3d[swap_mask][:, ::-1, :]

    canonical_2d = edges_2d.copy()
    canonical_2d[swap_mask] = canonical_2d[swap_mask][:, ::-1, :]

    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(
        view_bytes, return_index=True, return_inverse=True, return_counts=True
    )

    num_unique = len(unique_keys)

    max_cat = np.full(num_unique, -1, dtype=int)
    min_cat = np.full(num_unique, 999, dtype=int)
    np.maximum.at(max_cat, inverse_indices, edge_category)
    np.minimum.at(min_cat, inverse_indices, edge_category)

    rule1_unique = (counts >= 2) & (max_cat != min_cat)
    rule2_unique = counts == 1

    cos_threshold = math.cos(math.radians(max_angle_deg))
    rule3_unique = np.zeros(num_unique, dtype=bool)

    multi_edge_mask = counts >= 2
    if np.any(multi_edge_mask):
        sort_by_edge = np.argsort(inverse_indices)
        sorted_edge_ids = inverse_indices[sort_by_edge]
        sorted_normals = edge_normals[sort_by_edge]

        match = sorted_edge_ids[:-1] == sorted_edge_ids[1:]
        if np.any(match):
            n1 = sorted_normals[:-1][match]
            n2 = sorted_normals[1:][match]
            e_ids = sorted_edge_ids[:-1][match]

            dot_prods = np.sum(n1 * n2, axis=1)
            is_crease = dot_prods < (cos_threshold - 1e-6)
            rule3_unique[e_ids[is_crease]] = True

    marked_unique_mask = rule1_unique | rule2_unique | rule3_unique
    marked_first_indices = first_indices[marked_unique_mask]

    return canonical_2d[marked_first_indices]


def extract_boundary_edges(mesh_data: Any) -> NDArray[np.float64]:
    """
    Extracts 2D edge segments owned by only 1 face (single-face boundary edges).
    """
    from ...infra.blender.blender_mesh_capture import SlicingGroup, SlicingCollection, SlicingObject

    all_tris_3d = []

    def process_item(item: Any):
        if item is None:
            return
        elif isinstance(item, (list, tuple)):
            for sub in item:
                process_item(sub)
        elif isinstance(item, SlicingGroup):
            for coll in item.collections.values():
                process_item(coll)
        elif isinstance(item, SlicingCollection):
            for so in item.objects:
                process_item(so)
        elif isinstance(item, SlicingObject):
            if item.mesh.size > 0:
                all_tris_3d.append(item.mesh[:, :3, :3])

    process_item(mesh_data)

    if not all_tris_3d:
        return np.empty((0, 2, 2), dtype=np.float64)

    tris_3d = np.concatenate(all_tris_3d, axis=0)
    if len(tris_3d) == 0:
        return np.empty((0, 2, 2), dtype=np.float64)

    verts_2d = tris_3d[:, :, :2]
    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)

    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)

    dists_sq_2d = np.sum((edges_2d[:, 0, :] - edges_2d[:, 1, :]) ** 2, axis=1)
    valid_mask = dists_sq_2d >= 1e-8

    edges_3d = edges_3d[valid_mask]
    edges_2d = edges_2d[valid_mask]

    if len(edges_3d) == 0:
        return np.empty((0, 2, 2), dtype=np.float64)

    rounded_3d = np.round(edges_3d, decimals=4)
    pA = rounded_3d[:, 0, :]
    pB = rounded_3d[:, 1, :]
    swap_mask = (pA[:, 0] > pB[:, 0]) | (
        (pA[:, 0] == pB[:, 0]) & (pA[:, 1] > pB[:, 1])
    ) | (
        (pA[:, 0] == pB[:, 0]) & (pA[:, 1] == pB[:, 1]) & (pA[:, 2] > pB[:, 2])
    )

    canonical_3d = rounded_3d.copy()
    canonical_3d[swap_mask] = canonical_3d[swap_mask][:, ::-1, :]

    canonical_2d = edges_2d.copy()
    canonical_2d[swap_mask] = canonical_2d[swap_mask][:, ::-1, :]

    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(
        view_bytes, return_index=True, return_inverse=True, return_counts=True
    )

    boundary_mask = counts == 1
    return canonical_2d[first_indices[boundary_mask]]


def extract_crease_edges(mesh_data: Any, max_angle_deg: float = 30.0) -> NDArray[np.float64]:
    """
    Extracts 2D edge segments owned by 2 or more faces where the 3D dihedral angle between face normals > max_angle_deg.
    """
    return extract_lineart_edges(mesh_data, max_angle_deg=max_angle_deg)


def extract_silhouette_contour_edges(mesh_data: Any) -> NDArray[np.float64]:
    """
    Extracts top-view silhouette contour edges using 3D face normals.
    Filters front-facing faces (normal_z > epsilon) and extracts edges owned by exactly 1 front-facing face.
    """
    from ...infra.blender.blender_mesh_capture import SlicingGroup, SlicingCollection, SlicingObject

    all_tris_3d = []
    all_normals = []

    def process_item(item: Any):
        if item is None:
            return
        elif isinstance(item, (list, tuple)):
            for sub in item:
                process_item(sub)
        elif isinstance(item, SlicingGroup):
            for coll in item.collections.values():
                process_item(coll)
        elif isinstance(item, SlicingCollection):
            for so in item.objects:
                process_item(so)
        elif isinstance(item, SlicingObject):
            if item.mesh.size > 0:
                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)

    if not all_tris_3d:
        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)

    normal_z = normals_unit[:, 2]
    epsilon = 1e-3
    front_mask = normal_z > epsilon

    tris_3d = tris_3d[front_mask]
    if len(tris_3d) == 0:
        return np.empty((0, 2, 2), dtype=np.float64)

    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)  # (3M, 2, 3)

    verts_2d = tris_3d[:, :, :2]  # (M, 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)  # (3M, 2, 2)

    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]

    if len(edges_3d) == 0:
        return np.empty((0, 2, 2), dtype=np.float64)

    rounded_3d = np.round(edges_3d, decimals=4)
    pA = rounded_3d[:, 0, :]
    pB = rounded_3d[:, 1, :]
    swap_mask = (pA[:, 0] > pB[:, 0]) | (
        (pA[:, 0] == pB[:, 0]) & (pA[:, 1] > pB[:, 1])
    ) | (
        (pA[:, 0] == pB[:, 0]) & (pA[:, 1] == pB[:, 1]) & (pA[:, 2] > pB[:, 2])
    )

    canonical_3d = rounded_3d.copy()
    canonical_3d[swap_mask] = canonical_3d[swap_mask][:, ::-1, :]

    canonical_2d = edges_2d.copy()
    canonical_2d[swap_mask] = canonical_2d[swap_mask][:, ::-1, :]

    flat_3d = canonical_3d.reshape(len(canonical_3d), 6)
    view_bytes = np.ascontiguousarray(flat_3d).view(np.dtype((np.void, 48))).ravel()
    _, first_indices, counts = np.unique(
        view_bytes, return_index=True, return_counts=True
    )

    contour_edges = canonical_2d[first_indices[counts == 1]]
    return deduplicate_2d_edges(contour_edges)


def extract_outer_contour_edges(mesh_data: Any) -> NDArray[np.float64]:
    """
    Extracts top-view contour edges using backfacing CW vs CCW triangle orientation.
    """
    return extract_silhouette_contour_edges(mesh_data)