"""
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

from ..geometry import euler_rotation_matrix

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

    rot_mat: np.ndarray | None = None
    if has_rot:
        rot_mat = euler_rotation_matrix(r_arr[0], r_arr[1], r_arr[2])

    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 rot_mat is not None and dim > 0:
            if dim >= 3:
                new_p[:, :3] = new_p[:, :3] @ rot_mat.T
            elif dim == 2:
                rot_2d = rot_mat[: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,
    u0: float | None = None,
    v0: float | None = None,
    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).
    Major grid lines stem from origin u0, v0.
    """
    from ..slice import SlicePaths

    u_0 = u0 if u0 is not None else grid_min_x
    v_0 = v0 if v0 is not None else grid_min_y

    u_dir = 1.0 if grid_max_x >= u_0 else -1.0
    num_u_steps = int(np.ceil(abs(grid_max_x - u_0) / minor_step))
    if u_min_step := int(np.ceil(abs(grid_min_x - u_0) / minor_step)):
        num_u_steps = max(num_u_steps, u_min_step)

    v_dir = 1.0 if grid_max_y >= v_0 else -1.0
    num_v_steps = int(np.ceil(abs(grid_max_y - v_0) / minor_step))
    if v_min_step := int(np.ceil(abs(grid_min_y - v_0) / minor_step)):
        num_v_steps = max(num_v_steps, v_min_step)

    if num_u_steps > max_grid_steps or num_v_steps > max_grid_steps:
        return SlicePaths(), SlicePaths()

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

    # Lines parallel to V (constant U): gu_i = u_0 + u_dir * i * minor_step
    # Generate all u positions within [min(grid_min_x, grid_max_x), max(grid_min_x, grid_max_x)]
    min_u_val = min(grid_min_x, grid_max_x)
    max_u_val = max(grid_min_x, grid_max_x)
    min_v_val = min(grid_min_y, grid_max_y)
    max_v_val = max(grid_min_y, grid_max_y)

    k_u_start = int(np.floor((min_u_val - u_0) / minor_step))
    k_u_end = int(np.ceil((max_u_val - u_0) / minor_step))

    for i in range(k_u_start, k_u_end + 1):
        gu = u_0 + i * minor_step
        if gu < min_u_val - 1e-6 or gu > max_u_val + 1e-6:
            continue
        p1 = map_2d_to_3d(gu, min_v_val, z_level, plane)
        p2 = map_2d_to_3d(gu, max_v_val, z_level, plane)
        line = np.array([p1, p2], dtype=np.float64)
        if i % 10 == 0:
            main_paths.append(line)
        else:
            minor_paths.append(line)

    # Lines parallel to U (constant V): gv_j = v_0 + j * minor_step
    k_v_start = int(np.floor((min_v_val - v_0) / minor_step))
    k_v_end = int(np.ceil((max_v_val - v_0) / minor_step))

    for j in range(k_v_start, k_v_end + 1):
        gv = v_0 + j * minor_step
        if gv < min_v_val - 1e-6 or gv > max_v_val + 1e-6:
            continue
        p1 = map_2d_to_3d(min_u_val, gv, z_level, plane)
        p2 = map_2d_to_3d(max_u_val, gv, z_level, plane)
        line = np.array([p1, p2], dtype=np.float64)
        if j % 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]
                v0_3d = verts_3d[:, 0, :]
                v1_3d = verts_3d[:, 1, :]
                v2_3d = verts_3d[:, 2, :]
                if item.mesh.shape[1] >= 4:
                    normals = item.mesh[:, 3, :]
                else:
                    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)