from __future__ import annotations
import math
from numpy import float64
import numpy as np
import pyclipper
from typing import Any, Literal, cast
import copy
from shapely.geometry.base import BaseGeometry
from shapely.geometry import Polygon, MultiPolygon, LineString, LinearRing, MultiLineString, GeometryCollection

from ..infra.blender.blender_mesh_capture import PdfProperties

from .geometry import ndarray_lru_cache

# Scaling factor for pyclipper (it works with integers)
SCALE = 1_000_000

def _transfer_path_attr(src_attr: np.ndarray, orig_pts: np.ndarray, new_pts: np.ndarray) -> np.ndarray:
    """
    Transfers 1D per-vertex attribute from orig_pts (shape N, 3) to new_pts (shape M, 3).
    If src_attr is uniform (all values equal), returns array of shape (M,) filled with that value.
    Otherwise, interpolates along normalized cumulative length.
    """
    N = len(orig_pts)
    M = len(new_pts)
    if M == 0:
        return np.zeros(0, dtype=src_attr.dtype if hasattr(src_attr, 'dtype') else float)
    if N == 0:
        return np.zeros(M, dtype=src_attr.dtype if hasattr(src_attr, 'dtype') else float)
    if N == M and np.allclose(orig_pts, new_pts):
        return src_attr.copy()
    if len(src_attr) > 0 and np.all(src_attr == src_attr[0]):
        return np.full(M, src_attr[0], dtype=src_attr.dtype if hasattr(src_attr, 'dtype') else float)

    orig_diffs = np.diff(orig_pts[:, :3], axis=0)
    orig_dists = np.sqrt((orig_diffs ** 2).sum(axis=1))
    orig_cum = np.concatenate(([0.0], np.cumsum(orig_dists)))
    total_orig = orig_cum[-1]
    if total_orig > 1e-7:
        orig_t = orig_cum / total_orig
    else:
        orig_t = np.linspace(0.0, 1.0, N)

    new_diffs = np.diff(new_pts[:, :3], axis=0)
    new_dists = np.sqrt((new_diffs ** 2).sum(axis=1))
    new_cum = np.concatenate(([0.0], np.cumsum(new_dists)))
    total_new = new_cum[-1]
    if total_new > 1e-7:
        new_t = new_cum / total_new
    else:
        new_t = np.linspace(0.0, 1.0, M)

    return np.interp(new_t, orig_t, src_attr)


def _transfer_group_attributes(src_sp: SlicePaths, res_paths: list[np.ndarray], res_groups: list[int]) -> tuple[list[np.ndarray], list[np.ndarray], dict[str, list[np.ndarray]]]:
    """
    Given input SlicePaths src_sp and output res_paths/res_groups,
    constructs extrusion_rate, feedrate, and extra_attributes for res_paths
    by matching groups or nearest paths in src_sp.
    """
    out_ext: list[np.ndarray] = []
    out_feed: list[np.ndarray] = []
    out_extra: dict[str, list[np.ndarray]] = {k: [] for k in getattr(src_sp, "extra_attributes", {})}

    src_groups = getattr(src_sp, "groups", [])
    src_group_map: dict[int, list[int]] = {}
    for idx, g in enumerate(src_groups):
        src_group_map.setdefault(g, []).append(idx)

    for idx, (p, g) in enumerate(zip(res_paths, res_groups)):
        M = len(p)
        src_indices = src_group_map.get(g, [])
        if not src_indices and len(src_sp) > 0:
            z = p[0, 2] if len(p) > 0 else 0.0
            src_indices = [min(range(len(src_sp)), key=lambda i: abs(src_sp[i][0, 2] - z) if len(src_sp[i]) > 0 else float('inf'))]

        if src_indices:
            src_i = src_indices[0]
            ext = _transfer_path_attr(src_sp.extrusion_rate[src_i], src_sp[src_i], p)
            feed = _transfer_path_attr(src_sp.feedrate[src_i], src_sp[src_i], p)
            out_ext.append(ext)
            out_feed.append(feed)
            for k in src_sp.extra_attributes:
                attr_arr = _transfer_path_attr(src_sp.extra_attributes[k][src_i], src_sp[src_i], p)
                out_extra[k].append(attr_arr)
        else:
            out_ext.append(np.ones(M))
            out_feed.append(np.zeros(M))
            for k in src_sp.extra_attributes:
                out_extra[k].append(np.zeros(M))

    return out_ext, out_feed, out_extra


def align_path_seam_sharpest(path: np.ndarray) -> np.ndarray:
    """
    Reorders a closed 2D/3D polygon path (N, 3) or (N, 2) so that its starting vertex (seam)
    is located at the sharpest concave angle (inner corner). If no concave angles exist,
    falls back to the sharpest angle overall.
    
    If path is open or has fewer than 3 distinct vertices, returns path unchanged.
    Uses vectorized NumPy operations for high performance.
    """
    if len(path) < 4:
        return path

    if not np.allclose(path[0], path[-1]):
        return path

    pts = path[:-1]
    M = len(pts)
    if M < 3:
        return path

    pts2d = pts[:, :2]
    x = pts2d[:, 0]
    y = pts2d[:, 1]
    x_next = np.roll(x, -1)
    y_next = np.roll(y, -1)
    signed_area = 0.5 * np.sum(x * y_next - x_next * y)
    if abs(signed_area) < 1e-9:
        return path

    prev_pts = np.roll(pts2d, 1, axis=0)
    next_pts = np.roll(pts2d, -1, axis=0)

    v_in = pts2d - prev_pts
    v_out = next_pts - pts2d

    l_in = np.linalg.norm(v_in, axis=1)
    l_out = np.linalg.norm(v_out, axis=1)

    valid = (l_in > 1e-6) & (l_out > 1e-6)
    if not np.any(valid):
        return path

    cross_2d = v_in[:, 0] * v_out[:, 1] - v_in[:, 1] * v_out[:, 0]
    if signed_area > 0:
        concave_mask = valid & (cross_2d < -1e-6)
    else:
        concave_mask = valid & (cross_2d > 1e-6)

    cos_angles = np.full(M, 1.0)
    dot_prods = np.sum(v_in[valid] * v_out[valid], axis=1)
    cos_angles[valid] = dot_prods / (l_in[valid] * l_out[valid])

    if np.any(concave_mask):
        concave_indices = np.where(concave_mask)[0]
        best_idx = int(concave_indices[np.argmin(cos_angles[concave_mask])])
    else:
        best_idx = int(np.argmin(cos_angles))

    if best_idx == 0:
        return path

    return np.vstack([path[best_idx:-1], path[:best_idx + 1]])


class SlicePaths(list[np.ndarray]):
    groups: list[int]
    extrusion_rate: list[np.ndarray]
    feedrate: list[np.ndarray]
    extra_attributes: dict[str, list[np.ndarray]]
    layer_height: float
    first_layer_height: float | None
    transforms: list[dict[str, Any]]
    pdf_object_properties: list[tuple[str, str]]
    pdf_properties: PdfProperties
    name: str
    objects_info: dict[str, Any] | None

    def __init__(self, paths: list[np.ndarray] = [], groups: list[int] = [], extrusion_rate: list[np.ndarray] | None = None, feedrate: list[np.ndarray] | None = None, extra_attributes: dict[str, list[np.ndarray]] | None = None, layer_height: float = 0.2, first_layer_height: float | None = None, transforms: list[dict[str, Any]] | None = None, name: str = "", objects_info: dict[str, Any] | None = None):
        super().__init__(paths)
        
        self.groups = []
        if groups: self.groups = list(groups)
        else: self.groups = list(range(0,len(paths)))

        self.extrusion_rate = []
        if extrusion_rate: 
            self.extrusion_rate = list(extrusion_rate)
        else:
            self.extrusion_rate = [np.ones(len(p)) for p in paths]

        self.feedrate = []
        if feedrate:
            self.feedrate = list(feedrate)
        else:
            self.feedrate = [np.zeros(len(p)) for p in paths]
        
        self.extra_attributes = extra_attributes or {}
        self.layer_height = layer_height
        self.first_layer_height = first_layer_height
        self.transforms = list(transforms) if transforms is not None else []
        self.name = name
        self.objects_info = copy.deepcopy(objects_info) if objects_info is not None else None
    
    def copy_with(self, paths: list[np.ndarray], groups: list[int] | None = None, extrusion_rate: list[np.ndarray] | None = None, feedrate: list[np.ndarray] | None = None, extra_attributes: dict[str, list[np.ndarray]] | None = None) -> SlicePaths:
        """Returns a new SlicePaths object with new paths/groups, preserving all metadata and instance attributes."""
        target_groups = groups if groups is not None else (self.groups if len(self.groups) == len(paths) else list(range(len(paths))))

        if extrusion_rate is None or feedrate is None or extra_attributes is None:
            calc_ext, calc_feed, calc_extra = _transfer_group_attributes(self, paths, target_groups)
            if extrusion_rate is None:
                extrusion_rate = calc_ext
            if feedrate is None:
                feedrate = calc_feed
            if extra_attributes is None:
                extra_attributes = calc_extra

        new_obj = SlicePaths(paths, target_groups, extrusion_rate, feedrate, extra_attributes)
        for k, v in self.__dict__.items():
            if k not in ["groups", "extrusion_rate", "feedrate", "extra_attributes"]:
                setattr(new_obj, k, v)
        return new_obj

    def offset(self, offset_vec: np.ndarray) -> SlicePaths:
        """Returns a new SlicePaths object shifted by the given vector."""
        return self.copy_with([p + offset_vec for p in self])

    def transform(self, matrix: np.ndarray) -> SlicePaths:
        """Applies a transformation matrix to all paths."""
        return self.copy_with([p @ matrix.T for p in self])

    def buffer(self, distance: float, miter: str = 'round', rotation: Any = None) -> SlicePaths:
        """Buffers (offsets) the paths. Handles optional rotation for non-planar buffering."""
        p = self
        if rotation and any(rotation):
            from mathutils import Euler
            mat = Euler(rotation).to_matrix().inverted()
            p = self.transform(np.array(mat))

        res_paths, res_groups = buffer_slice(p, distance, miter=miter)
        new_paths = self.copy_with(res_paths, res_groups)

        if rotation and any(rotation):
            from mathutils import Euler
            mat_inv = Euler(rotation).to_matrix()
            new_paths = new_paths.transform(np.array(mat_inv))
        
        return new_paths

    def buffer_multiple(self, base_distance: float, count: int, width_step: float, miter: str = 'round', rotation: Any = None) -> SlicePaths:
        """Generates multiple concentric buffers and interleaves them by group (Z-layer), preserving all metadata."""
        if not self or count <= 0:
            return self.copy_with([])

        all_buffers: list[SlicePaths] = []
        for i in range(count):
            dist = base_distance + (i * width_step)
            all_buffers.append(self.buffer(dist, miter=miter, rotation=rotation))

        return SlicePaths.zip_paths(all_buffers, primary_source=self)

    def buffer_curve(self, distance: float, cap_style: str = 'round', join_style: str = 'round', single_sided: bool = False, rotation: Any = None) -> SlicePaths:
        """Buffers open curves using shapely. Handles optional rotation for non-planar buffering."""
        p = self
        if rotation and any(rotation):
            from mathutils import Euler
            mat = Euler(rotation).to_matrix().inverted()
            p = self.transform(np.array(mat))

        res_paths, res_groups = buffer_curve_slice(p, distance, cap_style=cap_style, join_style=join_style, single_sided=single_sided)
        new_paths = self.copy_with(res_paths, res_groups)

        if rotation and any(rotation):
            from mathutils import Euler
            mat_inv = Euler(rotation).to_matrix()
            new_paths = new_paths.transform(np.array(mat_inv))
        
        return new_paths

    def offset_curve(self, distance: float, join_style: Literal['round', 'mitre'] = 'round', rotation: Any = None) -> SlicePaths:
        """Offsets open/closed curves by a distance using shapely. Handles optional rotation."""
        p = self
        if rotation and any(rotation):
            from mathutils import Euler
            mat = Euler(rotation).to_matrix().inverted()
            p = self.transform(np.array(mat))

        res_paths, res_groups = offset_curve_slice(p, distance, join_style=join_style)
        new_paths = self.copy_with(res_paths, res_groups)

        if rotation and any(rotation):
            from mathutils import Euler
            mat_inv = Euler(rotation).to_matrix()
            new_paths = new_paths.transform(np.array(mat_inv))
        
        return new_paths

    def simplify(self, tolerance: float = 0.1, preserve_topology: bool = True) -> SlicePaths:
        """Simplifies paths using Shapely while conserving all per-path metadata and attributes."""
        if not self or tolerance <= 0:
            return self

        from shapely.geometry import LineString

        simplified_paths = []
        for p in self:
            if len(p) < 3:
                simplified_paths.append(p)
                continue
            line = LineString(p)
            simplified_line = line.simplify(tolerance, preserve_topology=preserve_topology)
            if simplified_line.geom_type == 'LineString' and len(simplified_line.coords) >= 2:
                simplified_paths.append(np.array(simplified_line.coords, dtype=p.dtype))
            else:
                simplified_paths.append(p)

        return self.copy_with(simplified_paths)

    @classmethod
    def extend_paths(cls, paths_list: list[SlicePaths], primary_source: SlicePaths | None = None) -> SlicePaths:
        """Combines multiple SlicePaths sequentially into one SlicePaths, conserving all metadata and attributes."""
        if not paths_list:
            return SlicePaths()

        base = primary_source if primary_source is not None else paths_list[0]
        merged = base.copy_with([])

        merged_transforms = []
        for sp in paths_list:
            merged.extend_with_groups(sp)
            for t in getattr(sp, "transforms", []):
                if t not in merged_transforms:
                    merged_transforms.append(t)
        merged.transforms = merged_transforms

        return merged

    @classmethod
    def zip_paths(cls, paths_list: list[SlicePaths], primary_source: SlicePaths | None = None) -> SlicePaths:
        """Interleaves multiple SlicePaths layer by layer (by group ID), conserving all metadata and attributes."""
        if not paths_list:
            return SlicePaths()

        base = primary_source if primary_source is not None else paths_list[0]
        
        all_groups = set()
        grouped_inputs = []
        for p_list in paths_list:
            d = {}
            for i, g in enumerate(p_list.groups):
                all_groups.add(g)
                d.setdefault(g, []).append((p_list[i], p_list.extrusion_rate[i], p_list.feedrate[i], {k: v[i] for k, v in p_list.extra_attributes.items()}))
            grouped_inputs.append(d)

        sorted_groups = sorted(list(all_groups))

        merged = base.copy_with([])
        for g in sorted_groups:
            for d in grouped_inputs:
                items = d.get(g, [])
                for path, ext, feed, extra in items:
                    merged.append_with_group(path, group=g, extrusion_rate=ext, feedrate=feed, extra_attributes=extra)

        merged_transforms = []
        merged_names = []
        merged_objs = []
        for sp in paths_list:
            obj_info = getattr(sp, "objects_info", None)
            if obj_info:
                for obj_entry in obj_info.get("objects", []):
                    if obj_entry not in merged_objs:
                        merged_objs.append(dict(obj_entry))
                        n_str = obj_entry.get("name", "")
                        if n_str and n_str not in merged_names:
                            merged_names.append(n_str)
            elif hasattr(sp, "name") and getattr(sp, "name"):
                n_str = str(getattr(sp, "name"))
                if n_str and n_str not in merged_names:
                    merged_names.append(n_str)
                    merged_objs.append({"name": n_str, "polygon": []})

            for t in getattr(sp, "transforms", []):
                if t not in merged_transforms:
                    merged_transforms.append(t)

        merged.transforms = merged_transforms
        if merged_objs:
            merged.objects_info = {"objects": merged_objs}
            merged.name = "-".join(merged_names)

        return merged

    def get_infill_mask(self) -> SlicePaths:
        """Computes the intersection of all paths in each group using shapely."""
        res_paths, res_groups = get_infill_mask_slice(self)
        return self.copy_with(res_paths, res_groups)

    def get_overhang_mask(self, tolerance: float = 0.05) -> SlicePaths:
        """Computes the overhang mask for each layer, propagating from top to bottom."""
        res_paths, res_groups = get_overhang_mask_slice(self, tuple(self.groups), tolerance=tolerance)
        return self.copy_with(res_paths, res_groups)

    def infill(self, density: float, rotation_deg: float = 0.0, infill_type: str = 'GRID', layer_rotation: Any = None) -> SlicePaths:
        """Generates infill for the paths."""
        p = self
        if layer_rotation and any(layer_rotation):
            from mathutils import Euler
            mat = Euler(layer_rotation).to_matrix().inverted()
            p = self.transform(np.array(mat))

        res_paths, res_groups = infill_slice(p, density, rotation_deg=rotation_deg, infill_type=infill_type)
        new_paths = self.copy_with(res_paths, res_groups)

        if layer_rotation and any(layer_rotation):
            from mathutils import Euler
            mat_inv = Euler(layer_rotation).to_matrix()
            new_paths = new_paths.transform(np.array(mat_inv))
        
        return new_paths

    def rotate(self, angle_deg: float, center: np.ndarray | None = None) -> SlicePaths:
        """Rotates the paths around Z-axis."""
        if angle_deg == 0:
            return self
            
        if center is None:
            center = self.get_first_layer_midpoint()
            
        angle_rad = np.radians(angle_deg)
        c, s = np.cos(angle_rad), np.sin(angle_rad)
        rot_mat = np.array([
            [c, -s, 0],
            [s,  c, 0],
            [0,  0, 1]
        ], dtype=np.float32)

        new_paths = []
        for p in self:
            # Shift to center, rotate, shift back
            new_p = (p - center) @ rot_mat.T + center
            new_paths.append(new_p)
            
        return self.copy_with(new_paths)

        new_paths = []
        for p in self:
            # Shift to center, rotate, shift back
            new_p = (p - center) @ rot_mat.T + center
            new_paths.append(new_p)
            
        return self.copy_with(new_paths)

    def get_first_layer_midpoint(self) -> np.ndarray:
        if not self:
            return np.zeros(3)
        
        # Find minimum group index
        min_g = min(self.groups)
        first_layer_indices = [i for i, g in enumerate(self.groups) if g == min_g]
        
        all_pts = []
        for i in first_layer_indices:
            all_pts.append(self[i])
            
        if not all_pts:
            return np.zeros(3)
            
        pts = np.concatenate(all_pts)
        min_pts = np.min(pts, axis=0)
        max_pts = np.max(pts, axis=0)
        return (min_pts + max_pts) / 2.0

    def boolean(self, other: SlicePaths, operation: str = 'UNION', connect_along_boundary: bool = False) -> SlicePaths:
        return boolean_slice(self, other, operation, connect_along_boundary=connect_along_boundary)

    def extend_with_groups(self, other: SlicePaths | list[np.ndarray]):
        if isinstance(other, SlicePaths):
            self.extend(other)
            self.groups.extend(other.groups)
            self.extrusion_rate.extend(other.extrusion_rate)
            self.feedrate.extend(other.feedrate)
            for k, v in other.extra_attributes.items():
                if k in self.extra_attributes:
                    self.extra_attributes[k].extend(v)
                else:
                    # Pad with zeros for existing paths if we didn't have this attribute
                    self.extra_attributes[k] = [np.zeros(len(p)) for p in self[:-len(other)]] + list(v)
            # Pad existing attributes if 'other' didn't have them
            for k in self.extra_attributes:
                if k not in other.extra_attributes:
                    self.extra_attributes[k].extend([np.zeros(len(p)) for p in other])
        else:
            self.extend(other)
            self.groups.extend(list(range(0,len(other))))
            self.extrusion_rate.extend([np.ones(len(p)) for p in other])
            self.feedrate.extend([np.zeros(len(p)) for p in other])
            for k in self.extra_attributes:
                self.extra_attributes[k].extend([np.zeros(len(p)) for p in other])

    def append_with_group(self, path: np.ndarray, group: int = 0, extrusion_rate: np.ndarray | None = None, feedrate: np.ndarray | None = None, extra_attributes: dict[str, np.ndarray] | None = None):
        self.append(path)
        self.groups.append(group)
        self.extrusion_rate.append(extrusion_rate if extrusion_rate is not None else np.ones(len(path)))
        self.feedrate.append(feedrate if feedrate is not None else np.zeros(len(path)))
        
        extra_attributes = extra_attributes or {}
        for k, v in extra_attributes.items():
            if k in self.extra_attributes:
                self.extra_attributes[k].append(v)
            else:
                self.extra_attributes[k] = [np.zeros(len(p)) for p in self[:-1]] + [v]
        
        for k in self.extra_attributes:
            if k not in extra_attributes:
                self.extra_attributes[k].append(np.zeros(len(path)))

    def normalize_groups(self) -> None:
        self.groups = [m.setdefault(x, len(m)) for m in [{}] for x in self.groups]

    def align_seams(self, mode: str = "sharpest") -> SlicePaths:
        """
        Aligns the seams (starting point of closed paths) in SlicePaths.
        Supported modes: 'sharpest'.
        """
        if not self or mode != "sharpest":
            return self

        new_paths: list[np.ndarray] = []
        new_ext: list[np.ndarray] = []
        new_feed: list[np.ndarray] = []
        new_extra: dict[str, list[np.ndarray]] = {k: [] for k in self.extra_attributes}

        for i, p in enumerate(self):
            if len(p) < 4 or not np.allclose(p[0], p[-1]):
                new_paths.append(p)
                new_ext.append(self.extrusion_rate[i])
                new_feed.append(self.feedrate[i])
                for k in self.extra_attributes:
                    new_extra[k].append(self.extra_attributes[k][i])
                continue

            pts = p[:-1]
            M = len(pts)
            if M < 3:
                new_paths.append(p)
                new_ext.append(self.extrusion_rate[i])
                new_feed.append(self.feedrate[i])
                for k in self.extra_attributes:
                    new_extra[k].append(self.extra_attributes[k][i])
                continue

            pts2d = pts[:, :2]
            x = pts2d[:, 0]
            y = pts2d[:, 1]
            x_next = np.roll(x, -1)
            y_next = np.roll(y, -1)
            signed_area = 0.5 * np.sum(x * y_next - x_next * y)

            prev_pts = np.roll(pts2d, 1, axis=0)
            next_pts = np.roll(pts2d, -1, axis=0)

            v_in = pts2d - prev_pts
            v_out = next_pts - pts2d

            l_in = np.linalg.norm(v_in, axis=1)
            l_out = np.linalg.norm(v_out, axis=1)

            valid = (l_in > 1e-6) & (l_out > 1e-6)
            if not np.any(valid) or abs(signed_area) < 1e-9:
                best_idx = 0
            else:
                cross_2d = v_in[:, 0] * v_out[:, 1] - v_in[:, 1] * v_out[:, 0]
                if signed_area > 0:
                    concave_mask = valid & (cross_2d < -1e-6)
                else:
                    concave_mask = valid & (cross_2d > 1e-6)

                cos_angles = np.full(M, 1.0)
                dot_prods = np.sum(v_in[valid] * v_out[valid], axis=1)
                cos_angles[valid] = dot_prods / (l_in[valid] * l_out[valid])

                if np.any(concave_mask):
                    concave_indices = np.where(concave_mask)[0]
                    best_idx = int(concave_indices[np.argmin(cos_angles[concave_mask])])
                else:
                    best_idx = int(np.argmin(cos_angles))

            if best_idx == 0:
                new_paths.append(p)
                new_ext.append(self.extrusion_rate[i])
                new_feed.append(self.feedrate[i])
                for k in self.extra_attributes:
                    new_extra[k].append(self.extra_attributes[k][i])
            else:
                new_p = np.vstack([p[best_idx:-1], p[:best_idx + 1]])
                new_paths.append(new_p)

                ext = self.extrusion_rate[i]
                if len(ext) == len(p):
                    new_ext.append(np.concatenate([ext[best_idx:-1], ext[:best_idx + 1]]))
                else:
                    new_ext.append(ext)

                feed = self.feedrate[i]
                if len(feed) == len(p):
                    new_feed.append(np.concatenate([feed[best_idx:-1], feed[:best_idx + 1]]))
                else:
                    new_feed.append(feed)

                for k in self.extra_attributes:
                    attr = self.extra_attributes[k][i]
                    if len(attr) == len(p):
                        new_extra[k].append(np.concatenate([attr[best_idx:-1], attr[:best_idx + 1]]))
                    else:
                        new_extra[k].append(attr)

        return self.copy_with(new_paths, extrusion_rate=new_ext, feedrate=new_feed, extra_attributes=new_extra)

def slice_at_range(models, z_values, ignore_holes: bool = False) -> list[list[np.ndarray]]:
    if not models:
        return []

    mmps: list[list[np.ndarray]] = []

    for model in models:
        model_edges: np.ndarray = edges_from_model(model)
        layered_edges: np.ndarray = layered_model(model_edges, z_values)
        sliced_edges: list[np.ndarray] = slice_at_z(layered_edges)
        
        all_model_paths: list[np.ndarray] = []

        for z, e in zip(z_values, sliced_edges):
            if e is None or len(e) == 0:
                continue

            zf = float(z)

            # Build 3D segments directly
            seg3d = e.copy()
            seg3d[..., 2] = zf

            # Build geometry for this slice
            slice_paths = segments_to_slice(seg3d, ignore_holes=ignore_holes)
            all_model_paths.extend(slice_paths)

        mmps.append(all_model_paths)
    return mmps

def edges_from_model(model): # returns a numpy: (facet, 3, 2, xyz)
    pairs_indices = np.array([[0, 1], [1, 2], [2, 0]]) # Create an array of indices to select the pairs
    return model[:, pairs_indices, :] # Use advanced indexing to create the array

def layered_model(models: np.ndarray, z_values: np.ndarray) -> np.ndarray:
    models_array: np.ndarray = np.array(models)[None, :]  # (1, F, 3, 2, 3)
    m = len(z_values)
    # array of zeros with shape (m, 1, 1, 1, 3)
    z_offsets_full = np.zeros((m, 1, 1, 1, 3))
    z_offsets_full[..., 2] = z_values[:, None, None, None]
    
    layered = models_array - z_offsets_full  # (m, F, 3, 2, 3)
    return layered

def slice_at_z(edges, normal_vector=np.array([0,0,1]), tolerance=1e-15):
    edges_offs = edges  # Shape (m, n, 3, 2, 3)

    # Obtain intersection points
    line_vector = edges_offs[..., 1, :] - edges_offs[..., 0, :]

    numerator = np.sum(normal_vector * -edges_offs[..., 0, :], axis=-1)
    denominator = np.sum(normal_vector[np.newaxis, np.newaxis, :] * line_vector, axis=-1)

    # Avoid division by zero
    with np.errstate(divide='ignore', invalid='ignore'):
        t = numerator / denominator
        t[denominator == 0] = np.nan  # Assign NaN where there's no valid intersection
    t[(t < 0) | (t > 1)] = np.nan

    intersections = edges_offs[..., 0, :] + t[..., np.newaxis] * line_vector
    intersections[np.abs(intersections) < tolerance] = 0  # Zero out small values

    all_segments = []
    for itx in intersections:
        mask = ~np.any(np.isnan(itx), axis=-1)
        filtered_arr = itx[mask]

        segment_count = filtered_arr.shape[-2] // 2
        segments = filtered_arr.reshape(segment_count, 2, 3)
        segments = np.around(segments, 6)
        all_segments.append(segments)

    return all_segments

def link_segments(segments: np.ndarray, tolerance=1e-6) -> list[np.ndarray]:
    """Links segments (N, 2, 3) into continuous paths."""
    if len(segments) == 0: return []
    z = segments[0, 0, 2]
    pts = segments[:, :, :2]
    
    def to_key(p): return tuple((p / tolerance).round().astype(np.int64))
    
    adj = {}
    for i in range(len(segments)):
        p1, p2 = pts[i, 0], pts[i, 1]
        k1, k2 = to_key(p1), to_key(p2)
        if k1 == k2: continue
        adj.setdefault(k1, []).append((k2, i))
        adj.setdefault(k2, []).append((k1, i))
        
    paths = []
    used = set()
    for i in range(len(segments)):
        if i in used: continue
        p1, p2 = pts[i, 0], pts[i, 1]
        k1, k2 = to_key(p1), to_key(p2)
        if k1 == k2: continue
        
        path_pts = [p1, p2]
        used.add(i)
        
        # Forward
        curr_k = k2
        while True:
            found = False
            for next_k, idx in adj.get(curr_k, []):
                if idx not in used:
                    pk1, pk2 = to_key(pts[idx, 0]), to_key(pts[idx, 1])
                    other_pt = pts[idx, 1] if pk1 == curr_k else pts[idx, 0]
                    path_pts.append(other_pt)
                    used.add(idx)
                    curr_k = to_key(other_pt)
                    found = True
                    break
            if not found: break
            
        # Backward
        curr_k = k1
        while True:
            found = False
            for next_k, idx in adj.get(curr_k, []):
                if idx not in used:
                    pk1, pk2 = to_key(pts[idx, 0]), to_key(pts[idx, 1])
                    other_pt = pts[idx, 1] if pk1 == curr_k else pts[idx, 0]
                    path_pts.insert(0, other_pt)
                    used.add(idx)
                    curr_k = to_key(other_pt)
                    found = True
                    break
            if not found: break
            
        res = np.column_stack([np.array(path_pts), np.full(len(path_pts), z)])
        paths.append(res)
    return paths

@ndarray_lru_cache(maxsize=128)
def segments_to_slice(segments: np.ndarray, ignore_holes: bool = False) -> list[np.ndarray]:

    linked = link_segments(segments)
    if not linked: return []
    
    closed_paths = []
    open_paths = []
    for p in linked:
        path_pts = p[:, :2]
        if len(path_pts) > 3 and np.allclose(path_pts[0], path_pts[-1]):
            closed_paths.append(p)
        else:
            open_paths.append(p)
            
    res = []
    z = segments[0, 0, 2]
    
    if closed_paths:
        pc = pyclipper.Pyclipper() # pyright: ignore[reportAttributeAccessIssue]
        for p in closed_paths:
            path_pts = p[:, :2]
            path_pts = path_pts[:-1]
            scaled = (path_pts * SCALE).astype(np.int64).tolist()
            pc.AddPath(scaled, pyclipper.PT_SUBJECT, closed = True) # pyright: ignore[reportAttributeAccessIssue]
        
        if ignore_holes:
            tree = pc.Execute2(pyclipper.CT_UNION, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD) # pyright: ignore[reportAttributeAccessIssue]
            def collect_outer_contours(node) -> list[list[int]]:
                contours = []
                if hasattr(node, 'Contour') and node.Contour and not getattr(node, 'IsHole', False):
                    contours.append(node.Contour)
                if hasattr(node, 'Childs') and node.Childs:
                    for child in node.Childs:
                        contours.extend(collect_outer_contours(child))
                return contours
            
            solution = collect_outer_contours(tree)
        else:
            # Use EVENODD for robust hole detection.
            solution = pc.Execute(pyclipper.CT_UNION, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD) # pyright: ignore[reportAttributeAccessIssue]
        
        for path in solution:
            if not path: continue
            fpath = np.array(path) / SCALE
            fpath_3d = np.column_stack([fpath, np.full(len(fpath), z)])
            # Close loop by appending the first point for internal consistency
            fpath_3d = np.vstack([fpath_3d, fpath_3d[0]])
            fpath_3d = align_path_seam_sharpest(fpath_3d)
            res.append(fpath_3d)
            
    res.extend(open_paths)
    return res

def find_point_on_loop(pt: np.ndarray, loop: np.ndarray, tol: float = 1.0) -> tuple[int, float]:
    """Finds the edge index and parameter u [0,1] for a point on a loop. loop is (N,2) and closed."""
    p1 = loop[:-1]
    p2 = loop[1:]
    v = p2 - p1
    w = pt - p1
    
    l2 = np.sum(v**2, axis=1)
    mask = l2 > 0
    u = np.zeros(len(l2))
    u[mask] = np.sum(w[mask] * v[mask], axis=1) / l2[mask]
    u = np.clip(u, 0, 1)
    
    dist_sq = np.sum((p1 + u[:, np.newaxis] * v - pt)**2, axis=1)
    min_idx = np.argmin(dist_sq)
    if dist_sq[min_idx] < tol**2:
        return int(min_idx), float(u[min_idx])
    return -1, 0.0

def get_boundary_path(pt1: np.ndarray, pt2: np.ndarray, loops: list[np.ndarray]) -> np.ndarray | None:
    """Returns the shortest path along one of the loops between pt1 and pt2."""
    best_path = None
    min_len = float('inf')
    
    for loop in loops:
        # Pyclipper paths are usually not closed in the array (last != first) if ET_CLOSEDPOLYGON
        # but NodeSlicer usually keeps them closed. Ensure it's closed for logic.
        working_loop = loop
        if not np.allclose(loop[0], loop[-1]):
            working_loop = np.vstack([loop, loop[0]])
            
        idx1, u1 = find_point_on_loop(pt1, working_loop)
        if idx1 == -1: continue
        idx2, u2 = find_point_on_loop(pt2, working_loop)
        if idx2 == -1: continue
        
        # Found on same loop. Try forward and backward.
        num_verts = len(working_loop) - 1
        
        # Forward path
        fwd_pts = [pt1]
        curr = idx1
        while curr != idx2:
            curr = (curr + 1) % num_verts
            fwd_pts.append(working_loop[curr])
        fwd_pts.append(pt2)
        fwd_path = np.array(fwd_pts)
        fwd_len = np.sum(np.sqrt(np.sum(np.diff(fwd_path, axis=0)**2, axis=1)))
        
        # Backward path
        bwd_pts = [pt1]
        curr = idx1
        while curr != (idx2 + 1) % num_verts:
            bwd_pts.append(working_loop[curr])
            curr = (curr - 1 + num_verts) % num_verts
        bwd_pts.append(pt2)
        bwd_path = np.array(bwd_pts)
        bwd_len = np.sum(np.sqrt(np.sum(np.diff(bwd_path, axis=0)**2, axis=1)))
        
        if fwd_len < bwd_len:
            if fwd_len < min_len:
                min_len = fwd_len
                best_path = fwd_path
        else:
            if bwd_len < min_len:
                min_len = bwd_len
                best_path = bwd_path
                
    return best_path

@ndarray_lru_cache(maxsize=128)
def boolean_slice(subject: SlicePaths, object: SlicePaths, operation: str = 'UNION', connect_along_boundary: bool = False) -> SlicePaths:
    all_z_keys: set[int] = set()
    subject_dict: dict[int, list[np.ndarray]] = {}
    object_dict: dict[int, list[np.ndarray]] = {}

    for s in subject:
        z = int(round(s[0, 2] * 1000))
        all_z_keys.add(z)
        subject_dict.setdefault(z, []).append(s)

    for s in object:
        z = int(round(s[0, 2] * 1000))
        all_z_keys.add(z)
        object_dict.setdefault(z, []).append(s)

    results: SlicePaths = SlicePaths()

    group = 0
    for k in sorted(all_z_keys):
        subject_slices = subject_dict.get(k, [])
        object_slices = object_dict.get(k, [])

        if not subject_slices:
            if object_slices and operation == 'UNION':
                results.extend(object_slices)
                results.groups.extend([group] * len(object_slices))
            continue

        if not object_slices:
            if operation != 'INTERSECTION':
                results.extend(subject_slices)
                results.groups.extend([group] * len(subject_slices))
            continue

        z = subject_slices[0][0, 2]

        pc = pyclipper.Pyclipper() # pyright: ignore[reportAttributeAccessIssue]

        # Function to process and add paths to Clipper
        def add_paths_to_clipper(paths, role):
            for s in paths:
                path_pts = s[:, :2]
                is_closed = len(path_pts) > 2 and np.allclose(path_pts[0], path_pts[-1])
                if is_closed:
                    path_pts = path_pts[:-1]
                scaled = (path_pts * SCALE).astype(np.int64).tolist()
                pc.AddPath(scaled, role, is_closed) # pyright: ignore[reportAttributeAccessIssue]

        layer_results = []

        def process_op(subs, clips_scaled, clips_pts, op_type):
            for s in subs:
                path_pts = s[:, :2]
                is_closed = len(path_pts) > 2 and np.allclose(path_pts[0], path_pts[-1])
                
                if not connect_along_boundary or is_closed:
                    pc.Clear()
                    if is_closed: path_pts = path_pts[:-1]
                    pc.AddPath((path_pts * SCALE).astype(np.int64).tolist(), pyclipper.PT_SUBJECT, is_closed)  # pyright: ignore[reportAttributeAccessIssue]
                    for c in clips_scaled:
                        pc.AddPath(c, pyclipper.PT_CLIP, True) # pyright: ignore[reportAttributeAccessIssue]
                    
                    solution = pc.Execute2(op_type, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD) # pyright: ignore[reportAttributeAccessIssue]
                    for path in pyclipper.ClosedPathsFromPolyTree(solution): # pyright: ignore[reportAttributeAccessIssue]
                        fpath = np.array(path) / SCALE
                        fpath_3d = np.column_stack([fpath, np.full(len(fpath), z)])
                        layer_results.append(np.vstack([fpath_3d, fpath_3d[0]]))
                    for path in pyclipper.OpenPathsFromPolyTree(solution): # pyright: ignore[reportAttributeAccessIssue]
                        fpath = np.array(path) / SCALE
                        layer_results.append(np.column_stack([fpath, np.full(len(fpath), z)]))
                else:
                    pc.Clear()
                    pc.AddPath((path_pts * SCALE).astype(np.int64).tolist(), pyclipper.PT_SUBJECT, False) # pyright: ignore[reportAttributeAccessIssue]
                    for c in clips_scaled:
                        pc.AddPath(c, pyclipper.PT_CLIP, True) # pyright: ignore[reportAttributeAccessIssue]
                    
                    solution = pc.Execute2(op_type, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD) # pyright: ignore[reportAttributeAccessIssue]
                    open_results = pyclipper.OpenPathsFromPolyTree(solution) # pyright: ignore[reportAttributeAccessIssue]
                    if not open_results: continue
                    
                    def get_p(pt):
                        p1 = path_pts[:-1]
                        p2 = path_pts[1:]
                        v = p2 - p1
                        w = pt - p1
                        l2 = np.sum(v**2, axis=1)
                        mask = l2 > 0
                        u = np.zeros(len(l2))
                        u[mask] = np.sum(w[mask] * v[mask], axis=1) / l2[mask]
                        u = np.clip(u, 0, 1)
                        dist_sq = np.sum((p1 + u[:, np.newaxis] * v - pt)**2, axis=1)
                        i = np.argmin(dist_sq)
                        return i + u[i]

                    segments = []
                    for path in open_results:
                        fpath = np.array(path) / SCALE
                        p_start = get_p(fpath[0])
                        p_end = get_p(fpath[-1])
                        if p_start > p_end:
                            fpath = fpath[::-1]
                            p_start, p_end = p_end, p_start
                        segments.append((p_start, fpath))
                    
                    segments.sort(key=lambda x: x[0])
                    
                    current_group = [segments[0][1]]
                    for i in range(1, len(segments)):
                        prev_end = current_group[-1][-1] * SCALE
                        curr_start = segments[i][1][0] * SCALE
                        scaled_loops = [(pts * SCALE) for pts in clips_pts]
                        b_path = get_boundary_path(prev_end, curr_start, scaled_loops)
                        
                        if b_path is not None:
                            if len(b_path) > 2:
                                current_group.append(b_path[1:-1] / SCALE)
                            current_group.append(segments[i][1])
                        else:
                            # Path is not continuous along boundary, finalize current group
                            path_to_add = np.vstack(current_group)
                            layer_results.append(np.column_stack([path_to_add, np.full(len(path_to_add), z)]))
                            current_group = [segments[i][1]]
                    
                    path_to_add = np.vstack(current_group)
                    layer_results.append(np.column_stack([path_to_add, np.full(len(path_to_add), z)]))

        if operation == 'INTERSECTION':
            obj_closed_pts = [s[:, :2] for s in object_slices if len(s) > 2 and np.allclose(s[0], s[-1])]
            obj_closed_scaled = [(pts[:-1] * SCALE).astype(np.int64).tolist() for pts in obj_closed_pts]
            process_op(subject_slices, obj_closed_scaled, obj_closed_pts, pyclipper.CT_INTERSECTION) # pyright: ignore[reportAttributeAccessIssue]
            
            # Pass 2 for B_open ∩ A_closed
            b_open = [s for s in object_slices if not (len(s) > 2 and np.allclose(s[0], s[-1]))]
            if b_open:
                sub_closed_pts = [s[:, :2] for s in subject_slices if len(s) > 2 and np.allclose(s[0], s[-1])]
                sub_closed_scaled = [(pts[:-1] * SCALE).astype(np.int64).tolist() for pts in sub_closed_pts]
                process_op(b_open, sub_closed_scaled, sub_closed_pts, pyclipper.CT_INTERSECTION) # pyright: ignore[reportAttributeAccessIssue]
        
        elif operation == 'DIFFERENCE':
            obj_closed_pts = [s[:, :2] for s in object_slices if len(s) > 2 and np.allclose(s[0], s[-1])]
            obj_closed_scaled = [(pts[:-1] * SCALE).astype(np.int64).tolist() for pts in obj_closed_pts]
            process_op(subject_slices, obj_closed_scaled, obj_closed_pts, pyclipper.CT_DIFFERENCE) # pyright: ignore[reportAttributeAccessIssue]
            # DIFFERENCE is not symmetric for open paths in the same way, usually A - B means A is the main focus.

        else:
            # UNION, XOR
            add_paths_to_clipper(subject_slices, pyclipper.PT_SUBJECT) # pyright: ignore[reportAttributeAccessIssue]
            obj_closed = []
            obj_open = []
            for s in object_slices:
                if len(s) > 2 and np.allclose(s[0], s[-1]): obj_closed.append(s)
                else: obj_open.append(s)

            if operation in ['UNION', 'XOR']:
                add_paths_to_clipper(obj_open, pyclipper.PT_SUBJECT) # pyright: ignore[reportAttributeAccessIssue]
            
            for s in obj_closed:
                path_pts = s[:, :2]
                scaled = (path_pts[:-1] * SCALE).astype(np.int64).tolist()
                pc.AddPath(scaled, pyclipper.PT_CLIP, True) # pyright: ignore[reportAttributeAccessIssue]
            
            op_type = pyclipper.CT_UNION # pyright: ignore[reportAttributeAccessIssue]
            if operation == 'UNION': op_type = pyclipper.CT_UNION # pyright: ignore[reportAttributeAccessIssue]
            elif operation == 'XOR': op_type = pyclipper.CT_XOR # pyright: ignore[reportAttributeAccessIssue]

            solution = pc.Execute2(op_type, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD) # pyright: ignore[reportAttributeAccessIssue]
            for path in pyclipper.ClosedPathsFromPolyTree(solution): # pyright: ignore[reportAttributeAccessIssue]
                fpath = np.array(path) / SCALE
                fpath_3d = np.column_stack([fpath, np.full(len(fpath), z)])
                layer_results.append(np.vstack([fpath_3d, fpath_3d[0]]))
            for path in pyclipper.OpenPathsFromPolyTree(solution): # pyright: ignore[reportAttributeAccessIssue]
                fpath = np.array(path) / SCALE
                layer_results.append(np.column_stack([fpath, np.full(len(fpath), z)]))

        results.extend(layer_results)
        results.groups.extend([group] * len(layer_results))
        group += 1
            
    # Transfer per-group/path metadata (extrusion_rate, feedrate, extra_attributes) from subject
    calc_ext, calc_feed, calc_extra = _transfer_group_attributes(subject, results, results.groups)
    results.extrusion_rate = calc_ext
    results.feedrate = calc_feed
    results.extra_attributes = calc_extra

    for attr, val in subject.__dict__.items():
        if attr not in ["groups", "extrusion_rate", "feedrate", "extra_attributes"]:
            setattr(results, attr, val)
    for attr, val in object.__dict__.items():
        if attr not in ["groups", "extrusion_rate", "feedrate", "extra_attributes"]:
            if not hasattr(results, attr):
                setattr(results, attr, val)

    return results

def geometry_to_paths(geom: Any, z: float) -> list[np.ndarray]:
    """Helper to convert shapely geometries to a list of 3D numpy arrays with a given Z height."""
    paths = []
    if geom.is_empty:
        return paths
    
    if isinstance(geom, Polygon):
        ext_coords = np.array(geom.exterior.coords)
        ext_coords_3d = np.column_stack([ext_coords, np.full(len(ext_coords), z)])
        paths.append(ext_coords_3d)
        for interior in geom.interiors:
            int_coords = np.array(interior.coords)
            int_coords_3d = np.column_stack([int_coords, np.full(len(int_coords), z)])
            paths.append(int_coords_3d)
    elif isinstance(geom, (MultiPolygon, GeometryCollection)):
        for poly in geom.geoms:
            paths.extend(geometry_to_paths(poly, z))
    
    return paths

@ndarray_lru_cache(maxsize=64)
def buffer_curve_slice(
    slice_paths: list[np.ndarray],
    distance: float,
    cap_style: str = 'round',
    join_style: str = 'round',
    single_sided: bool = False
) -> tuple[list[np.ndarray], list[int]]:
    if not slice_paths or distance <= 0:
        return [], []

    from shapely.geometry import LineString
    
    # Group paths by Z height to avoid flattening the entire model into 2D
    groups: dict[int, list[np.ndarray]] = {}
    for p in slice_paths:
        z = p[0, 2]
        # Use scaled integer for stable dictionary grouping of floats
        z_key = int(round(z * 1000))
        groups.setdefault(z_key, []).append(p)

    paths_list: list[np.ndarray] = []
    groups_list: list[int] = []

    group_id = 0
    for z_key in sorted(groups.keys()):
        layer_paths = groups[z_key]
        z = layer_paths[0][0, 2]

        for p in layer_paths:
            if len(p) < 2:
                continue

            pts_2d = p[:, :2]
            line = LineString(pts_2d)
            
            poly = line.buffer(
                distance,
                cap_style=cast(Literal['round', 'square', 'flat'], cap_style),
                join_style=cast(Literal['round', 'mitre', 'bevel'], join_style),
                single_sided=single_sided,
                quad_segs=16
            )
            
            layer_result_paths = geometry_to_paths(poly, z)
            
            for path in layer_result_paths:
                paths_list.append(path)
                groups_list.append(group_id)

        group_id += 1

    return paths_list, groups_list

def offset_geometry_to_paths(geom: Any, z: float) -> list[np.ndarray]:
    """Helper to convert shapely line geometries to a list of 3D numpy arrays with a given Z height."""
    paths = []
    if geom.is_empty:
        return paths
    
    if isinstance(geom, (LineString, LinearRing)):
        coords = np.array(geom.coords)
        coords_3d = np.column_stack([coords, np.full(len(coords), z)])
        paths.append(coords_3d)
    elif isinstance(geom, MultiLineString):
        for line in geom.geoms:
            coords = np.array(line.coords)
            coords_3d = np.column_stack([coords, np.full(len(coords), z)])
            paths.append(coords_3d)
    elif isinstance(geom, GeometryCollection):
        for sub_geom in geom.geoms:
            paths.extend(offset_geometry_to_paths(sub_geom, z))
    elif isinstance(geom, (Polygon, MultiPolygon)):
        boundary = geom.boundary
        paths.extend(offset_geometry_to_paths(boundary, z))
    return paths

@ndarray_lru_cache(maxsize=64)
def offset_curve_slice(
    slice_paths: list[np.ndarray],
    distance: float,
    join_style: Literal['round', 'mitre'] = 'round'
) -> tuple[list[np.ndarray], list[int]]:
    if not slice_paths:
        return [], []

    from shapely.geometry import LineString, Polygon, LinearRing
    
    # Group paths by Z height
    groups: dict[int, list[np.ndarray]] = {}
    for p in slice_paths:
        z = p[0, 2]
        z_key = int(round(z * 1000))
        groups.setdefault(z_key, []).append(p)

    paths_list: list[np.ndarray] = []
    groups_list: list[int] = []

    group_id = 0
    for z_key in sorted(groups.keys()):
        layer_paths = groups[z_key]
        z = layer_paths[0][0, 2]

        for p in layer_paths:
            if len(p) < 2:
                continue

            pts_2d = p[:, :2]
            
            # Check if path is closed
            is_closed = len(pts_2d) > 2 and np.allclose(pts_2d[0], pts_2d[-1])
            
            if is_closed:
                try:
                    # Check direction/orientation
                    ccw = LinearRing(pts_2d).is_ccw
                    poly = Polygon(pts_2d)
                    if poly.is_valid:
                        # Use Polygon buffer which keeps boundary closed
                        buf_dist = distance if ccw else -distance
                        offset_geom = poly.buffer(
                            buf_dist,
                            join_style=join_style,
                            quad_segs=16
                        )
                    else:
                        line = LineString(pts_2d)
                        offset_geom = line.offset_curve(
                            distance,
                            join_style=join_style,
                            quad_segs=16
                        )
                except Exception:
                    line = LineString(pts_2d)
                    offset_geom = line.offset_curve(
                        distance,
                        join_style=join_style,
                        quad_segs=16
                    )
            else:
                line = LineString(pts_2d)
                offset_geom = line.offset_curve(
                    distance,
                    join_style=join_style,
                    quad_segs=16
                )
            
            layer_result_paths = offset_geometry_to_paths(offset_geom, z)
            
            for path in layer_result_paths:
                paths_list.append(path)
                groups_list.append(group_id)

        group_id += 1

    return paths_list, groups_list

@ndarray_lru_cache(maxsize=64)
def get_infill_mask_slice(slice_paths: list[np.ndarray]) -> tuple[list[np.ndarray], list[int]]:
    if not slice_paths:
        return [], []

    from shapely.geometry import Polygon
    
    # Group paths by their Z height
    groups: dict[int, list[np.ndarray]] = {}
    for p in slice_paths:
        z = p[0, 2]
        z_key = int(round(z * 1000))
        groups.setdefault(z_key, []).append(p)

    paths_list: list[np.ndarray] = []
    groups_list: list[int] = []

    group_id = 0
    for z_key in sorted(groups.keys()):
        layer_paths = groups[z_key]
        if not layer_paths:
            continue
            
        z = layer_paths[0][0, 2]
        
        # Convert each closed path to a Shapely Polygon
        polygons = []
        for p in layer_paths:
            if len(p) < 3:
                continue
            pts_2d = p[:, :2]
            poly = Polygon(pts_2d)
            if poly.is_valid and not poly.is_empty:
                polygons.append(poly)
                
        if not polygons:
            continue
            
        # Compute the intersection of all polygons in the layer
        intersection_poly = polygons[0]
        for poly in polygons[1:]:
            intersection_poly = intersection_poly.intersection(poly)
            
        # Convert resulting polygon back to 3D paths
        result_paths = geometry_to_paths(intersection_poly, z)
        
        for path in result_paths:
            paths_list.append(path)
            groups_list.append(group_id)
            
        group_id += 1
            
    return paths_list, groups_list

def paths_to_shapely_polygon(paths: list[np.ndarray]) -> BaseGeometry:
    """Converts a list of paths on a layer to a single Shapely geometry using pyclipper to resolve holes."""
    if not paths:
        return Polygon()
        
    import pyclipper
    from shapely.ops import unary_union
    
    pc = pyclipper.Pyclipper()  # pyright: ignore[reportAttributeAccessIssue]
    for p in paths:
        if len(p) < 3:
            continue
        pts_2d = p[:, :2]
        if np.allclose(pts_2d[0], pts_2d[-1]):
            pts_2d = pts_2d[:-1]
        scaled = (pts_2d * SCALE).astype(np.int64).tolist()
        pc.AddPath(scaled, pyclipper.PT_SUBJECT, True) # pyright: ignore[reportAttributeAccessIssue]
        
    tree = pc.Execute2(pyclipper.CT_UNION, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD) # pyright: ignore[reportAttributeAccessIssue]
    
    def polytree_to_shapely(node) -> list[Polygon]:
        polygons = []
        if hasattr(node, 'Contour') and node.Contour and not node.IsHole:
            exterior = np.array(node.Contour) / SCALE
            interiors = []
            child_islands = []
            for child in node.Childs:
                if child.IsHole:
                    if child.Contour:
                        interiors.append(np.array(child.Contour) / SCALE)
                    for grandchild in child.Childs:
                        child_islands.append(grandchild)
                else:
                    child_islands.append(child)
            polygons.append(Polygon(exterior, interiors))
            for island in child_islands:
                polygons.extend(polytree_to_shapely(island))
        else:
            for child in node.Childs:
                polygons.extend(polytree_to_shapely(child))
        return polygons
        
    polys = polytree_to_shapely(tree)
    return unary_union(polys)

@ndarray_lru_cache(maxsize=32)
def get_overhang_mask_slice(slice_paths: list[np.ndarray], groups: tuple[int, ...], tolerance: float = 0.05) -> tuple[list[np.ndarray], list[int]]:
    if not slice_paths or not groups:
        return [], []

    from shapely.geometry import Polygon
    
    # Group paths by their group index
    group_map: dict[int, list[np.ndarray]] = {}
    for p, g in zip(slice_paths, groups):
        group_map.setdefault(g, []).append(p)
        
    sorted_groups = sorted(group_map.keys())
    if len(sorted_groups) < 2:
        # If there's only 1 layer, there are no overhangs under any layer
        return [], []

    # Pre-build solid geometry for each group
    layer_geoms = {g: paths_to_shapely_polygon(group_map[g]) for g in sorted_groups}

    paths_list: list[np.ndarray] = []
    groups_list: list[int] = []
    
    accumulated = Polygon()
    
    # Iterate from top group down to bottom group
    for i in range(len(sorted_groups) - 1, -1, -1):
        g = sorted_groups[i]
        z = group_map[g][0][0, 2] if group_map[g] and len(group_map[g][0]) > 0 else 0.0
        geom_i = layer_geoms[g]
        
        # Calculate raw support needed at layer g: union of all layers above g, minus geom_i.
        # This prevents boundary mismatches/slivers caused by subtracting and then unioning back the same layers.
        if i < len(sorted_groups) - 1:
            raw_support = accumulated.difference(geom_i)
            if not raw_support.is_empty:
                # Apply tolerance clearance only to the output of the current layer.
                geom_i_sub = geom_i.buffer(tolerance) if tolerance > 0 else geom_i
                output_geom = raw_support.difference(geom_i_sub)
                if not output_geom.is_empty:
                    result_paths = geometry_to_paths(output_geom, z)
                    for path in result_paths:
                        paths_list.append(path)
                        groups_list.append(g)
                        
        # Accumulate the union of model layers for the layers below
        accumulated = accumulated.union(geom_i)
                
    # Sort output by group index so that paths are ordered bottom-to-top
    sorted_indices = np.argsort(groups_list)
    paths_list = [paths_list[idx] for idx in sorted_indices]
    groups_list = [int(groups_list[idx]) for idx in sorted_indices]
    
    return paths_list, groups_list

@ndarray_lru_cache(maxsize=64)
def buffer_slice(slice_paths: list[np.ndarray], distance: float, miter='round', tolerance=.01) -> tuple[list[np.ndarray], list[int]]:

    if not slice_paths: return [], []
    
    # Group paths by Z height to avoid flattening the entire model into 2D
    groups: dict[int, list[np.ndarray]] = {}
    for p in slice_paths:
        z = p[0, 2]
        # Use scaled integer for stable dictionary grouping of floats
        z_key = int(round(z * 1000))
        groups.setdefault(z_key, []).append(p)
    
    join_types = {
        'round': pyclipper.JT_ROUND, # pyright: ignore[reportAttributeAccessIssue]
        'square': pyclipper.JT_SQUARE, # pyright: ignore[reportAttributeAccessIssue]
        'flat': pyclipper.JT_SQUARE, # pyright: ignore[reportAttributeAccessIssue]
        'miter': pyclipper.JT_MITER # pyright: ignore[reportAttributeAccessIssue]
    }
    jt = join_types.get(miter, pyclipper.JT_ROUND) # pyright: ignore[reportAttributeAccessIssue]
    
    paths_list: list[np.ndarray] = []
    groups_list: list[int] = []
    pco = pyclipper.PyclipperOffset(arc_tolerance=tolerance * SCALE) # pyright: ignore[reportAttributeAccessIssue]
    
    # Buffer each layer independently
    group_id = 0
    for z_key in sorted(groups.keys()):
        layer_paths = groups[z_key]
        z = layer_paths[0][0, 2]
        
        pco.Clear()
        scaled_paths = []
        for p in layer_paths:
            path_pts = p[:, :2]
            # Strip redundant closing points
            if len(path_pts) > 2 and np.allclose(path_pts[0], path_pts[-1]):
                path_pts = path_pts[:-1]
            
            if len(path_pts) < 3: continue
            scaled_paths.append((path_pts * SCALE).astype(np.int64).tolist())
        
        if not scaled_paths: continue

        # Add all paths for this specific Z layer
        pco.AddPaths(scaled_paths, jt, pyclipper.ET_CLOSEDPOLYGON) # pyright: ignore[reportAttributeAccessIssue]
        solution = pco.Execute(distance * SCALE)
        
        for path in solution:
            if not path: continue
            fpath = np.array(path) / SCALE
            fpath_3d = np.column_stack([fpath, np.full(len(fpath), z)])
            # Re-close loop
            fpath_3d = np.vstack([fpath_3d, fpath_3d[0]])
            fpath_3d = align_path_seam_sharpest(fpath_3d)
            paths_list.append(fpath_3d)
            groups_list.append(group_id)
        
        group_id += 1
            
    return paths_list, groups_list

@ndarray_lru_cache(maxsize=32)
def infill_slice(slice_paths: 'SlicePaths', density: float, rotation_deg: float = 0.0, infill_type: str = 'GRID') -> tuple[list[np.ndarray], list[int]]:
    """Generates a grid or rectilinear infill clipped by the provided paths. Consistent across layers."""
    if not slice_paths or density <= 0: return [], []
    
    rotation = np.radians(rotation_deg)
    
    # Calculate global bounding box for the entire model to ensure grid alignment
    all_pts_model = np.concatenate(slice_paths)
    g_min_x, g_min_y = all_pts_model[:, 0].min(), all_pts_model[:, 1].min()
    g_max_x, g_max_y = all_pts_model[:, 0].max(), all_pts_model[:, 1].max()
    global_center = np.array([(g_min_x + g_max_x) / 2, (g_min_y + g_max_y) / 2, 0])

    # Pre-calculate rotation matrices
    c_inv, s_inv = np.cos(-rotation), np.sin(-rotation)
    rot_mat_inv = np.array([[c_inv, -s_inv, 0], [s_inv, c_inv, 0], [0, 0, 1]])
    
    c, s = np.cos(rotation), np.sin(rotation)
    rot_mat = np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])

    # Determine grid bounds in the rotated space to ensure full coverage
    corners = np.array([
        [g_min_x, g_min_y, 0], [g_max_x, g_min_y, 0],
        [g_min_x, g_max_y, 0], [g_max_x, g_max_y, 0]
    ])
    rotated_corners = (corners - global_center) @ rot_mat_inv.T + global_center
    rg_min_x, rg_min_y = rotated_corners[:, 0].min(), rotated_corners[:, 1].min()
    rg_max_x, rg_max_y = rotated_corners[:, 0].max(), rotated_corners[:, 1].max()

    # Generate the global static grid lines in de-rotated space
    rg_min_x -= density
    rg_max_x += density
    rg_min_y -= density
    rg_max_y += density

    vertical_lines = []
    x = rg_min_x
    while x <= rg_max_x:
        vertical_lines.append([
            (int(x * SCALE), int(rg_min_y * SCALE)),
            (int(x * SCALE), int(rg_max_y * SCALE))
        ])
        x += density

    horizontal_lines = []
    y = rg_min_y
    while y <= rg_max_y:
        horizontal_lines.append([
            (int(rg_min_x * SCALE), int(y * SCALE)),
            (int(rg_max_x * SCALE), int(y * SCALE))
        ])
        y += density

    # Group by Z
    groups: dict[int, list[np.ndarray]] = {}
    if hasattr(slice_paths, 'groups'):
        for p, g in zip(slice_paths, slice_paths.groups):
            groups.setdefault(g, []).append(p)
    else:
        # Fallback if it's just a list
        for p in slice_paths:
            z_key = int(round(p[0, 2] * 1000))
            groups.setdefault(z_key, []).append(p)
        
    output_paths: list[np.ndarray] = []
    output_groups: list[int] = []

    # Sort groups by Z to determine layer index for alternating rectilinear infill
    unique_groups = sorted(groups.keys())
    group_to_layer_index = {g: i for i, g in enumerate(unique_groups)}

    for g in groups.keys():
        layer_idx = group_to_layer_index[g]
        
        if infill_type == 'GRID':
            current_grid_paths = vertical_lines + horizontal_lines
        elif infill_type == 'RECTILINEAR':
            current_grid_paths = vertical_lines if layer_idx % 2 == 0 else horizontal_lines
        elif infill_type == 'LINES':
            current_grid_paths = vertical_lines
        else:
            current_grid_paths = []

        layer_paths = groups[g]
        z = layer_paths[0][0, 2]
        
        # De-rotate layer paths
        offset_2d = np.array([global_center[0], global_center[1], 0])
        active_paths = [(p - offset_2d) @ rot_mat_inv.T + offset_2d for p in layer_paths]

        # Clip lines with active_paths
        pc = pyclipper.Pyclipper() # pyright: ignore[reportAttributeAccessIssue]
        
        for p in active_paths:
            path_pts = p[:, :2]
            if len(path_pts) > 2 and np.allclose(path_pts[0], path_pts[-1]):
                path_pts = path_pts[:-1]
            if len(path_pts) < 3: continue
            scaled = (path_pts * SCALE).astype(np.int64).tolist()
            pc.AddPath(scaled, pyclipper.PT_CLIP, True) # pyright: ignore[reportAttributeAccessIssue]
            
        if current_grid_paths:
            pc.AddPaths(current_grid_paths, pyclipper.PT_SUBJECT, False) # pyright: ignore[reportAttributeAccessIssue]
            
        solution = pc.Execute2(pyclipper.CT_INTERSECTION, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD) # pyright: ignore[reportAttributeAccessIssue]
        clipped_lines = pyclipper.OpenPathsFromPolyTree(solution)  # pyright: ignore[reportAttributeAccessIssue]
        
        # Rotate back
        for path in clipped_lines:
            if not path: continue
            fpath = np.array(path) / SCALE
            fpath_3d = np.column_stack([fpath, np.full(len(fpath), z)])
            rotated_fpath_3d = (fpath_3d - offset_2d) @ rot_mat.T + offset_2d
            output_paths.append(rotated_fpath_3d)
            output_groups.append(g)
        
    return output_paths, output_groups

def group_segments_by_z(segments: np.ndarray) -> dict[int, np.ndarray]:
    z_values: np.ndarray = segments[:, 0, 2]
    unique_z: np.ndarray = np.unique(z_values)
    groups: dict[int, np.ndarray] = {int(round(z*1000)): segments[z_values == z] for z in unique_z}
    return groups

############################################################

class Pather(object):

    def __init__(self, paths: list[np.ndarray], **kwargs) -> None:
        self.paths = paths
        
        # Robust initialization: ensure metadata lists match paths length
        e_rate = getattr(paths, "extrusion_rate", None)
        if e_rate is not None and len(e_rate) == len(paths):
            self.extrusion_rates = e_rate
        else:
            self.extrusion_rates = [np.ones(len(p)) for p in paths]
            
        f_rate = getattr(paths, "feedrate", None)
        if f_rate is not None and len(f_rate) == len(paths):
            self.feedrates = f_rate
        else:
            self.feedrates = [np.zeros(len(p)) for p in paths]
            
        extra = getattr(paths, "extra_attributes", None)
        if extra is not None and all(len(v) == len(paths) for v in extra.values()):
            self.extra_attributes = extra
        else:
            self.extra_attributes = {}
            
        self.bed_center_x: float = 250.
        self.bed_center_y: float = 250.
        self.feedrate: float = 1200.
        self.travel_feedrate: float = 2400.
        self.extrusion_rate: float = 0.017
        self.layer_height: float = 0.2
        self.first_layer_height: float | None = None
        self._layer_width: float | None = None
        self.retraction: float = 0.0
        self.retraction_feedrate: float = 0.0
        self.lift: float = 0.0
        self.transforms: list[dict[str, Any]] = list(getattr(paths, "transforms", []))
        self.objects_info: dict[str, Any] | None = getattr(paths, "objects_info", None)
        self.name: str = str(getattr(paths, "name", ""))

    @property
    def layer_width(self) -> float:
        if self._layer_width is not None:
            return self._layer_width
        if self.layer_height > 0:
            import math
            filament_area = math.pi * (1.75 / 2.0) ** 2
            return (self.extrusion_rate * filament_area) / self.layer_height
        return 0.4

    @layer_width.setter
    def layer_width(self, val: float | None):
        self._layer_width = val

    def to_gcode(self, include_header: bool = True) -> str:
        if not self.paths:
            return ""

        # Centering logic using NumPy
        all_pts = np.concatenate(self.paths)
        min_xy = all_pts[:, :2].min(axis=0)
        max_xy = all_pts[:, :2].max(axis=0)
        center_xy = (min_xy + max_xy) / 2.0

        offset_x = self.bed_center_x - center_xy[0]
        offset_y = self.bed_center_y - center_xy[1]

        # Apply offset to all paths
        centered_paths = []
        for path in self.paths:
            new_path = path.copy()
            new_path[:, 0] += offset_x
            new_path[:, 1] += offset_y
            centered_paths.append(new_path)

        try:
            retraction_fr_val = float(self.retraction_feedrate) if self.retraction_feedrate is not None else 0.0
        except (TypeError, ValueError):
            retraction_fr_val = 0.0

        if retraction_fr_val <= 0.0:
            retraction_fr_val = self.feedrate

        gcode: list[str] = []

        if include_header:
            header: list[str] = [
                "; Generated from NodeSlicer",
                f";TYPE:Custom",
                f";DEFAULT_FEEDRATE={self.feedrate}",
                f";DEFAULT_TRAVEL_FEEDRATE={self.travel_feedrate}",
                f";DEFAULT_RETRACTION_FEEDRATE={retraction_fr_val}",
                f";DEFAULT_EXTRUSION_RATE={self.extrusion_rate}",
                f";DEFAULT_LAYER_HEIGHT={self.layer_height:.6f}",
            ]
            if self.first_layer_height is not None:
                header.append(f";DEFAULT_FIRST_LAYER_HEIGHT={self.first_layer_height:.6f}")
            header.extend([
                f";DEFAULT_LAYER_WIDTH={self.layer_width:.6f}",
                f";OFFSET_X={offset_x:.6f}",
                f";OFFSET_Y={offset_y:.6f}",
            ])

            if self.objects_info:
                import json
                header.append(f"; objects_info = {json.dumps(self.objects_info)}")
            elif self.name:
                import json
                header.append(f'; objects_info = {json.dumps({"objects": [{"name": self.name, "polygon": []}]})}')

            if hasattr(self, "transforms") and self.transforms:
                import json
                header.append(f";PREPARATION_TRANSFORMS={json.dumps(self.transforms)}")

                total_tx, total_ty, total_tz = 0.0, 0.0, 0.0
                total_rx, total_ry, total_rz = 0.0, 0.0, 0.0
                total_sx, total_sy, total_sz = 1.0, 1.0, 1.0

                for i, t in enumerate(self.transforms):
                    t_type = t.get("type", "transform")
                    tr = t.get("translation", (0.0, 0.0, 0.0))
                    rot = t.get("rotation", (0.0, 0.0, 0.0))
                    sc = t.get("scale", (1.0, 1.0, 1.0))

                    total_tx += tr[0]
                    total_ty += tr[1]
                    total_tz += tr[2]
                    total_rx += rot[0]
                    total_ry += rot[1]
                    total_rz += rot[2]
                    total_sx *= sc[0]
                    total_sy *= sc[1]
                    total_sz *= sc[2]

                    if t_type == "transform_geometry":
                        header.append(f";TRANSFORM_GEOMETRY_{i}_TRANSLATION={tr[0]:.6f},{tr[1]:.6f},{tr[2]:.6f}")
                        header.append(f";TRANSFORM_GEOMETRY_{i}_ROTATION={rot[0]:.6f},{rot[1]:.6f},{rot[2]:.6f}")
                        header.append(f";TRANSFORM_GEOMETRY_{i}_SCALE={sc[0]:.6f},{sc[1]:.6f},{sc[2]:.6f}")
                    elif t_type == "drop_geometry":
                        target_z = t.get("target_z", 0.0)
                        off_z = t.get("offset_z", 0.0)
                        header.append(f";DROP_GEOMETRY_{i}_TARGET_Z={target_z:.6f}")
                        header.append(f";DROP_GEOMETRY_{i}_OFFSET_Z={off_z:.6f}")

                header.append(f";PREPARATION_TOTAL_TRANSLATION={total_tx:.6f},{total_ty:.6f},{total_tz:.6f}")
                header.append(f";PREPARATION_TOTAL_ROTATION={total_rx:.6f},{total_ry:.6f},{total_rz:.6f}")
                header.append(f";PREPARATION_TOTAL_SCALE={total_sx:.6f},{total_sy:.6f},{total_sz:.6f}")

            header.extend([
                "G21 ; set units to mm",
                "G90 ; absolute XYZ positioning",
                "M83 ; relative extrusion mode",
                "G92 E0",
            ])
            gcode += header

        # Safely parse retraction and lift fields.
        # If they cannot be parsed or value is zero, we skip them.
        try:
            retraction_val = float(self.retraction) if self.retraction is not None else 0.0
        except (TypeError, ValueError):
            retraction_val = 0.0

        try:
            lift_val = float(self.lift) if self.lift is not None else 0.0
        except (TypeError, ValueError):
            lift_val = 0.0

        last_z = -1.0
        max_z = 0
        last_nozzle_z = 0.0

        for i, line in enumerate(centered_paths):
            coords = line.tolist()
            if len(coords) < 2:
                continue

            e_rate = self.extrusion_rates[i]
            f_rate = self.feedrates[i]
            extras = {prefix: vals[i] for prefix, vals in self.extra_attributes.items()}

            x0, y0, z0 = coords[0]

            active_feedrate = self.feedrate
            if i > 0:
                # 1. Retract if retraction is enabled
                if retraction_val > 0.0:
                    gcode.append(f"G1 E-{retraction_val:.5f} F{retraction_fr_val:.3f}")
                    active_feedrate = retraction_fr_val
                
                # 2. Lift Z if lift is enabled
                if lift_val > 0.0:
                    gcode.append(f"G0 Z{last_nozzle_z + lift_val:.3f} F720")
                    active_feedrate = 720.0

                # 3. Travel to the start of path `i` (x0, y0) in XY
                travel_fr = f_rate[0] if f_rate[0] > 0 else self.travel_feedrate
                gcode.append(f"G0 X{x0:.3f} Y{y0:.3f} F{travel_fr}")
                active_feedrate = travel_fr

                # 4. Handle Z alignment (layer change and/or un-lift)
                if abs(z0 - last_z) > 0.001:
                    gcode.append(f";BEFORE_LAYER_CHANGE")
                    gcode.append(f"G0 Z{z0:.3f} F720")
                    gcode.append(f";AFTER_LAYER_CHANGE")
                    last_z = z0
                    active_feedrate = 720.0
                elif lift_val > 0.0:
                    gcode.append(f"G0 Z{z0:.3f} F720")
                    active_feedrate = 720.0

                # 5. Un-retract if retraction was enabled
                if retraction_val > 0.0:
                    gcode.append(f"G1 E{retraction_val:.5f} F{retraction_fr_val:.3f}")
                    active_feedrate = retraction_fr_val
            else:
                if abs(z0 - last_z) > 0.001:
                    gcode.append(f";BEFORE_LAYER_CHANGE")
                    gcode.append(f"G0 Z{z0:.3f} F720")
                    gcode.append(f";AFTER_LAYER_CHANGE")
                    last_z = z0
                    active_feedrate = 720.0

                old_fr = f_rate[0] if f_rate[0] > 0 else self.travel_feedrate
                gcode.append(f"G0 X{x0:.3f} Y{y0:.3f} Z{z0:.3f} F{old_fr}")
                active_feedrate = old_fr

            old_fr = active_feedrate
            old_x, old_y, old_z = x0, y0, z0
            old_extras = {prefix: None for prefix in extras}
            for j, (p1, p2) in enumerate(zip(coords, coords[1:])):
                seg_len: float = math.dist(p1, p2)
                e_segment = seg_len * self.extrusion_rate * e_rate[j]
                
                # Per-segment feedrate value (not multiplier)
                fr = f_rate[j] if f_rate[j] > 0 else self.feedrate

                old_x, old_y, old_z = p1
                x, y, z = p2

                cur_gcode = "G1"

                if not old_x == x: cur_gcode += f" X{x:.3f}"
                if not old_y == y: cur_gcode += f" Y{y:.3f}"
                if not old_z == z: cur_gcode += f" Z{z:.3f}"
                if e_segment: cur_gcode += f" E{e_segment:.5f}"
                if not old_fr == fr: cur_gcode += f" F{fr:.3f}"

                for prefix, evals in extras.items():
                    val = evals[j]
                    if not old_extras[prefix] == val:
                        cur_gcode += f" {prefix}{val:.3f}"
                    old_extras[prefix] = val

                gcode += [cur_gcode]

                max_z = max(max_z, z)
                old_fr = fr

            last_nozzle_z = coords[-1][2]

        return "\n".join(gcode)