nodeslicer/core/slice.py

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 import Polygon, MultiPolygon, LineString, LinearRing, MultiLineString, GeometryCollection, Point
import manifold3d
import trimesh

from ..infra.blender.blender_mesh_capture import PdfProperties

from .geometry import ndarray_lru_cache, update_transform

# Scaling factor for 1 micron unit integer math (1 mm = 1000 microns)
SCALE = 1000

def is_closed_path(path: np.ndarray) -> bool:
    """Checks if a 2D or 3D path forms a closed loop in integer micron space (1 micron unit)."""
    if len(path) < 3:
        return False
    return (int(round(float(path[0, 0]) * SCALE)) == int(round(float(path[-1, 0]) * SCALE)) and
            int(round(float(path[0, 1]) * SCALE)) == int(round(float(path[-1, 1]) * SCALE)))

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 len(src_attr) > 0 and (len(src_attr) == 1 or 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)
    if N == M and is_closed_path(orig_pts) == is_closed_path(new_pts) and np.all(np.rint(orig_pts * SCALE) == np.rint(new_pts * SCALE)):
        return src_attr.copy()

    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)

    src_z_keys: list[int] | None = None
    if len(src_sp) > 0:
        src_z_keys = [int(round(float(src_sp[i][0, 2]) * SCALE)) if len(src_sp[i]) > 0 and src_sp[i].shape[1] > 2 else 0 for i in range(len(src_sp))]

    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 and src_z_keys is not None:
            z_um = int(round(float(p[0, 2]) * SCALE)) if len(p) > 0 and p.shape[1] > 2 else 0
            src_indices = [min(range(len(src_sp)), key=lambda i: abs(src_z_keys[i] - z_um))]

        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


import ctypes
import os

_c_seam_lib = None
try:
    _so_dir = os.path.dirname(os.path.abspath(__file__))
    _so_path = os.path.join(_so_dir, "libseam_opt.so")
    _c_src = os.path.join(_so_dir, "c_seam.c")
    rebuild = False
    if os.path.exists(_c_src):
        if not os.path.exists(_so_path) or os.path.getmtime(_so_path) < os.path.getmtime(_c_src):
            rebuild = True
    if rebuild:
        import subprocess
        for compiler in ("gcc", "clang"):
            for flags in (["-O3", "-fopenmp", "-fPIC", "-shared"], ["-O3", "-fPIC", "-shared"]):
                try:
                    subprocess.run([compiler, *flags, _c_src, "-o", _so_path, "-lm"], check=True, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
                    if os.path.exists(_so_path):
                        break
                except Exception:
                    pass
            if os.path.exists(_so_path):
                break
    if os.path.exists(_so_path):
        _c_seam_lib = ctypes.CDLL(_so_path)
        _c_seam_lib.find_sharpest_seam_c.argtypes = [
            ctypes.POINTER(ctypes.c_double),
            ctypes.c_int,
            ctypes.c_int,
            ctypes.c_double
        ]
        _c_seam_lib.find_sharpest_seam_c.restype = ctypes.c_int

        if hasattr(_c_seam_lib, "weld_vertices_c"):
            _c_seam_lib.weld_vertices_c.argtypes = [
                ctypes.POINTER(ctypes.c_float),
                ctypes.c_int,
                ctypes.POINTER(ctypes.c_float),
                ctypes.POINTER(ctypes.c_uint32)
            ]
            _c_seam_lib.weld_vertices_c.restype = ctypes.c_int

        if hasattr(_c_seam_lib, "polygon_signed_area_2d_c"):
            _c_seam_lib.polygon_signed_area_2d_c.argtypes = [
                ctypes.POINTER(ctypes.c_double),
                ctypes.c_int,
                ctypes.c_int
            ]
            _c_seam_lib.polygon_signed_area_2d_c.restype = ctypes.c_double

        if hasattr(_c_seam_lib, "find_path_mesh_intersection_c"):
            _c_seam_lib.find_path_mesh_intersection_c.argtypes = [
                ctypes.POINTER(ctypes.c_double),
                ctypes.c_int,
                ctypes.c_int,
                ctypes.POINTER(ctypes.c_double),
                ctypes.c_int,
                ctypes.c_int,
                ctypes.POINTER(ctypes.c_double),
                ctypes.POINTER(ctypes.c_int),
                ctypes.POINTER(ctypes.c_double),
                ctypes.POINTER(ctypes.c_double)
            ]
            _c_seam_lib.find_path_mesh_intersection_c.restype = ctypes.c_int

        if hasattr(_c_seam_lib, "reposition_path_seam_c"):
            _c_seam_lib.reposition_path_seam_c.argtypes = [
                ctypes.POINTER(ctypes.c_double),
                ctypes.c_int,
                ctypes.c_int,
                ctypes.c_int,
                ctypes.c_double,
                ctypes.POINTER(ctypes.c_double),
                ctypes.POINTER(ctypes.c_double)
            ]
            _c_seam_lib.reposition_path_seam_c.restype = ctypes.c_int

        if hasattr(_c_seam_lib, "slice_triangles_at_z_c"):
            _c_seam_lib.slice_triangles_at_z_c.argtypes = [
                ctypes.POINTER(ctypes.c_double),
                ctypes.c_int,
                ctypes.c_double,
                ctypes.POINTER(ctypes.c_double)
            ]
            _c_seam_lib.slice_triangles_at_z_c.restype = ctypes.c_int

        if hasattr(_c_seam_lib, "intersect_path_with_2d_segments_c"):
            _c_seam_lib.intersect_path_with_2d_segments_c.argtypes = [
                ctypes.POINTER(ctypes.c_double),
                ctypes.c_int,
                ctypes.c_int,
                ctypes.POINTER(ctypes.c_double),
                ctypes.c_int,
                ctypes.c_double,
                ctypes.c_int,
                ctypes.POINTER(ctypes.c_double),
                ctypes.POINTER(ctypes.c_int),
                ctypes.POINTER(ctypes.c_double),
                ctypes.POINTER(ctypes.c_double)
            ]
            _c_seam_lib.intersect_path_with_2d_segments_c.restype = ctypes.c_int

        if hasattr(_c_seam_lib, "shift_path_seam_by_distance_c"):
            _c_seam_lib.shift_path_seam_by_distance_c.argtypes = [
                ctypes.POINTER(ctypes.c_double),
                ctypes.c_int,
                ctypes.c_int,
                ctypes.c_double,
                ctypes.POINTER(ctypes.c_double)
            ]
            _c_seam_lib.shift_path_seam_by_distance_c.restype = ctypes.c_int

        if hasattr(_c_seam_lib, "rdp_simplify_batch_c"):
            _c_seam_lib.rdp_simplify_batch_c.argtypes = [
                ctypes.POINTER(ctypes.POINTER(ctypes.c_double)),
                ctypes.POINTER(ctypes.c_int),
                ctypes.POINTER(ctypes.c_int),
                ctypes.c_int,
                ctypes.c_double,
                ctypes.c_int,
                ctypes.POINTER(ctypes.POINTER(ctypes.c_double)),
                ctypes.POINTER(ctypes.c_int)
            ]
            _c_seam_lib.rdp_simplify_batch_c.restype = None

        if hasattr(_c_seam_lib, "rdp_simplify_polyline_c"):
            _c_seam_lib.rdp_simplify_polyline_c.argtypes = [
                ctypes.POINTER(ctypes.c_double),
                ctypes.c_int,
                ctypes.c_int,
                ctypes.c_double,
                ctypes.c_int,
                ctypes.POINTER(ctypes.c_double)
            ]
            _c_seam_lib.rdp_simplify_polyline_c.restype = ctypes.c_int
except Exception:
    _c_seam_lib = None


def weld_vertices(tris: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """
    Given an array of triangles of shape (F, 3, 3) or (F, 4, 3), welds identical vertices.
    Returns (unique_verts, faces) where unique_verts is (V, 3) float32 and faces is (F, 3) uint32.
    """
    flat = np.ascontiguousarray(tris[:, :3, :].reshape(-1, 3), dtype=np.float32)
    num_verts = len(flat)
    if num_verts == 0:
        return np.empty((0, 3), dtype=np.float32), np.empty((0, 3), dtype=np.uint32)

    if _c_seam_lib is not None and hasattr(_c_seam_lib, "weld_vertices_c"):
        try:
            out_unique = np.empty((num_verts, 3), dtype=np.float32)
            out_inv = np.empty(num_verts, dtype=np.uint32)
            unique_count = int(_c_seam_lib.weld_vertices_c(
                flat.ctypes.data_as(ctypes.POINTER(ctypes.c_float)),
                num_verts,
                out_unique.ctypes.data_as(ctypes.POINTER(ctypes.c_float)),
                out_inv.ctypes.data_as(ctypes.POINTER(ctypes.c_uint32))
            ))
            if unique_count >= 0:
                return out_unique[:unique_count], out_inv.reshape(-1, 3)
        except Exception:
            pass

    # NumPy fallback
    void_dtype = np.dtype((np.void, flat.dtype.itemsize * flat.shape[1]))
    flat_void = flat.view(void_dtype).ravel()
    unique_voids, inv = np.unique(flat_void, return_inverse=True)
    unique_verts = unique_voids.view(flat.dtype).reshape(-1, 3)
    faces = inv.reshape(-1, 3).astype(np.uint32)
    return unique_verts, faces


def find_path_seam(path: np.ndarray) -> int:
    """
    Finds the optimal starting vertex index (seam) positioned at the back of the object (+Y).
    Uses a weighted heuristic prioritizing maximum Y with a sharpness bonus for corners,
    with deterministic tie-breaking for consistent vertical layer alignment.
    """
    if len(path) < 4 or not is_closed_path(path):
        return 0

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

    if _c_seam_lib is not None:
        try:
            c_pts = np.ascontiguousarray(pts, dtype=np.float64)
            stride = int(c_pts.shape[1])
            ptr = c_pts.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
            return int(_c_seam_lib.find_sharpest_seam_c(ptr, M, stride, float(SCALE)))
        except Exception:
            pass

    # NumPy fallback
    x = pts[:, 0]
    y = pts[:, 1]
    y_min, y_max = float(np.min(y)), float(np.max(y))
    y_range = y_max - y_min
    bonus = min(2.0, 0.1 * y_range) if y_range > 0 else 2.0

    prev_pts = np.roll(pts[:, :2], 1, axis=0)
    next_pts = np.roll(pts[:, :2], -1, axis=0)

    v_in = pts[:, :2] - prev_pts
    v_out = next_pts - pts[:, :2]

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

    valid = (l_in > 1e-7) & (l_out > 1e-7)
    cos_ang = np.ones(M, dtype=np.float64)
    cos_ang[valid] = np.sum(v_in[valid] * v_out[valid], axis=1) / (l_in[valid] * l_out[valid])
    cos_ang = np.clip(cos_ang, -1.0, 1.0)
    sharpness = 1.0 - cos_ang

    scores = y + sharpness * bonus
    # Deterministic tie-break: highest score, smallest x
    return int(np.lexsort((x, -scores))[0])


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 back of the object (highest Y with corner preference).
    
    If path is open or has fewer than 3 distinct vertices, returns path unchanged.
    """
    if len(path) < 4 or not is_closed_path(path):
        return path

    M = len(path) - 1
    best_idx = find_path_seam(path)
    if best_idx == 0:
        return path

    res = np.empty_like(path)
    split = M - best_idx
    res[:split] = path[best_idx:-1]
    res[split:] = path[:best_idx + 1]
    return res


def extract_mesh_triangles(mesh_input: Any) -> np.ndarray:
    """
    Extracts 3D triangle vertices (F, 3, 3) from a SlicingObject, SlicingCollection,
    SlicingGroup, list of geometries, bpy Object, or raw numpy array.
    """
    if mesh_input is None:
        return np.empty((0, 3, 3), dtype=np.float64)

    if isinstance(mesh_input, np.ndarray):
        if mesh_input.ndim == 3 and mesh_input.shape[1] >= 3 and mesh_input.shape[2] == 3:
            return np.ascontiguousarray(mesh_input[:, :3, :], dtype=np.float64)
        return np.empty((0, 3, 3), dtype=np.float64)

    objects = []
    if hasattr(mesh_input, "collections"):  # SlicingGroup
        for c in mesh_input.collections.values():
            objects.extend(getattr(c, "objects", []))
    elif hasattr(mesh_input, "objects"):  # SlicingCollection
        objects.extend(mesh_input.objects)
    elif hasattr(mesh_input, "mesh"):  # SlicingObject
        objects.append(mesh_input)
    elif isinstance(mesh_input, list):
        for item in mesh_input:
            if hasattr(item, "collections"):
                for c in item.collections.values():
                    objects.extend(getattr(c, "objects", []))
            elif hasattr(item, "objects"):
                objects.extend(item.objects)
            elif hasattr(item, "mesh"):
                objects.append(item)
    elif hasattr(mesh_input, "type") and getattr(mesh_input, "type") == 'MESH':
        from ..infra.blender.blender_mesh_capture import SlicingObject
        objects.append(SlicingObject(mesh_input))

    if not objects:
        return np.empty((0, 3, 3), dtype=np.float64)

    valid_meshes = [o.mesh[:, :3, :] for o in objects if hasattr(o, "mesh") and o.mesh is not None and o.mesh.size > 0]
    if not valid_meshes:
        return np.empty((0, 3, 3), dtype=np.float64)

    return np.ascontiguousarray(np.concatenate(valid_meshes, axis=0), dtype=np.float64)


def extract_mesh_transforms(mesh_input: Any) -> list[dict[str, Any]]:
    """
    Extracts accumulated transform dicts from a SlicingObject, SlicingCollection,
    SlicingGroup, or list of objects.
    """
    if mesh_input is None:
        return []
    if hasattr(mesh_input, "transforms") and isinstance(mesh_input.transforms, list) and mesh_input.transforms:
        return [dict(t) for t in mesh_input.transforms if isinstance(t, dict)]
    if hasattr(mesh_input, "collections"):
        for c in mesh_input.collections.values():
            if hasattr(c, "transforms") and c.transforms:
                return [dict(t) for t in c.transforms if isinstance(t, dict)]
            for o in getattr(c, "objects", []):
                if hasattr(o, "transforms") and o.transforms:
                    return [dict(t) for t in o.transforms if isinstance(t, dict)]
    if hasattr(mesh_input, "objects"):
        for o in mesh_input.objects:
            if hasattr(o, "transforms") and o.transforms:
                return [dict(t) for t in o.transforms if isinstance(t, dict)]
    if hasattr(mesh_input, "location") and hasattr(mesh_input, "rotation_euler") and hasattr(mesh_input, "scale"):
        return [{
            "type": "INIT",
            "translation": tuple(float(v) for v in mesh_input.location),
            "rotation": tuple(float(v) for v in mesh_input.rotation_euler),
            "scale": tuple(float(v) for v in mesh_input.scale),
        }]
    if isinstance(mesh_input, list):
        for item in mesh_input:
            trans = extract_mesh_transforms(item)
            if trans:
                return trans
    return []


def extract_geometry_transforms(
    geom: Any
) -> tuple[tuple[float, float, float], tuple[float, float, float], tuple[float, float, float]] | list[Any]:
    """
    Extracts consolidated (translation, rotation, scale) from a geometry item
    (SlicePaths, SlicingObject, SlicingCollection, SlicingGroup, or Blender Object).
    Rotation is returned as Euler angles in radians (rx, ry, rz).
    If geom is a list/tuple of geometry items (not a single SlicePaths),
    returns a list of consolidated (translation, rotation, scale) tuples.
    """
    if geom is None:
        return (0.0, 0.0, 0.0), (0.0, 0.0, 0.0), (1.0, 1.0, 1.0)

    # Check if geom is a list or tuple of geometry items (excluding SlicePaths)
    if isinstance(geom, (list, tuple)) and type(geom).__name__ != "SlicePaths" and not getattr(geom, "is_slice_paths", False):
        return [extract_geometry_transforms(item) for item in geom]
    elif hasattr(geom, "__iter__") and not hasattr(geom, "shape") and not isinstance(geom, (str, bytes, dict, set)) and type(geom).__name__ != "SlicePaths" and not getattr(geom, "is_slice_paths", False) and not hasattr(geom, "mesh") and not hasattr(geom, "objects") and not hasattr(geom, "collections") and not hasattr(geom, "location"):
        return [extract_geometry_transforms(item) for item in geom]

    transforms = extract_mesh_transforms(geom)
    if not transforms and hasattr(geom, "transforms") and isinstance(geom.transforms, list):
        transforms = [dict(t) for t in geom.transforms if isinstance(t, dict)]

    if not transforms and hasattr(geom, "location") and hasattr(geom, "rotation_euler") and hasattr(geom, "scale"):
        t = tuple(float(v) for v in geom.location)
        r = tuple(float(v) for v in geom.rotation_euler)
        s = tuple(float(v) for v in geom.scale)
        return (t[0], t[1], t[2]), (r[0], r[1], r[2]), (s[0], s[1], s[2])

    cur_t = [0.0, 0.0, 0.0]
    cur_r = [0.0, 0.0, 0.0]
    cur_s = [1.0, 1.0, 1.0]

    if transforms:
        for t in transforms:
            if not isinstance(t, dict):
                continue
            tr = t.get("translation", t.get("offset", (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))
            cur_t[0] += float(tr[0])
            cur_t[1] += float(tr[1])
            cur_t[2] += float(tr[2])
            cur_r[0] += float(rot[0])
            cur_r[1] += float(rot[1])
            cur_r[2] += float(rot[2])
            cur_s[0] *= float(sc[0])
            cur_s[1] *= float(sc[1])
            cur_s[2] *= float(sc[2])
            if "angle_deg" in t and not any(float(v) != 0.0 for v in rot):
                cur_r[2] += math.radians(float(t["angle_deg"]))

    return (
        (float(cur_t[0]), float(cur_t[1]), float(cur_t[2])),
        (float(cur_r[0]), float(cur_r[1]), float(cur_r[2])),
        (float(cur_s[0]), float(cur_s[1]), float(cur_s[2])),
    )


def compute_delta_transforms(
    path_transforms: list[dict[str, Any]] | None,
    mesh_transforms: list[dict[str, Any]] | None
) -> list[dict[str, Any]]:
    """
    Computes the difference between path_transforms and mesh_transforms,
    so that only missing transforms are applied to the cutting mesh.
    Prevents double-transformation when the mesh was already transformed upstream.
    """
    if not path_transforms:
        return []
    if not mesh_transforms:
        return path_transforms

    # Consolidate path transforms
    p_t = [0.0, 0.0, 0.0]
    p_r = [0.0, 0.0, 0.0]
    p_s = [1.0, 1.0, 1.0]
    p_angle_deg = 0.0
    p_center = [0.0, 0.0, 0.0]

    for t in path_transforms:
        if not isinstance(t, dict):
            continue
        tr = t.get("translation", t.get("offset", (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))
        p_t[0] += float(tr[0])
        p_t[1] += float(tr[1])
        p_t[2] += float(tr[2])
        p_r[0] += float(rot[0])
        p_r[1] += float(rot[1])
        p_r[2] += float(rot[2])
        p_s[0] *= float(sc[0])
        p_s[1] *= float(sc[1])
        p_s[2] *= float(sc[2])
        if "angle_deg" in t:
            p_angle_deg += float(t["angle_deg"])
            if "center" in t:
                p_center = [float(v) for v in t["center"]]

    # Consolidate mesh transforms
    m_t = [0.0, 0.0, 0.0]
    m_r = [0.0, 0.0, 0.0]
    m_s = [1.0, 1.0, 1.0]
    m_angle_deg = 0.0
    m_center = [0.0, 0.0, 0.0]

    for t in mesh_transforms:
        if not isinstance(t, dict):
            continue
        tr = t.get("translation", t.get("offset", (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))
        m_t[0] += float(tr[0])
        m_t[1] += float(tr[1])
        m_t[2] += float(tr[2])
        m_r[0] += float(rot[0])
        m_r[1] += float(rot[1])
        m_r[2] += float(rot[2])
        m_s[0] *= float(sc[0])
        m_s[1] *= float(sc[1])
        m_s[2] *= float(sc[2])
        if "angle_deg" in t:
            m_angle_deg += float(t["angle_deg"])
            if "center" in t:
                m_center = [float(v) for v in t["center"]]

    # Compute deltas
    d_t = (p_t[0] - m_t[0], p_t[1] - m_t[1], p_t[2] - m_t[2])
    d_r = (p_r[0] - m_r[0], p_r[1] - m_r[1], p_r[2] - m_r[2])
    d_s = (
        p_s[0] / m_s[0] if abs(m_s[0]) > 1e-12 else 1.0,
        p_s[1] / m_s[1] if abs(m_s[1]) > 1e-12 else 1.0,
        p_s[2] / m_s[2] if abs(m_s[2]) > 1e-12 else 1.0,
    )
    d_angle = p_angle_deg - m_angle_deg

    delta_transforms: list[dict[str, Any]] = []

    has_scale = any(abs(v - 1.0) > 1e-6 for v in d_s)
    has_rot = any(abs(v) > 1e-6 for v in d_r)
    has_trans = any(abs(v) > 1e-6 for v in d_t)
    has_angle = abs(d_angle) > 1e-6

    if has_scale or has_rot or has_trans or has_angle:
        dt: dict[str, Any] = {"type": "transform"}
        if has_scale:
            dt["scale"] = d_s
        if has_rot:
            dt["rotation"] = d_r
        if has_trans:
            dt["translation"] = d_t
        if has_angle:
            dt["angle_deg"] = d_angle
            dt["center"] = p_center
        delta_transforms.append(dt)

    return delta_transforms


def apply_transforms_to_triangles(
    tris: np.ndarray,
    transforms: list[dict[str, Any]] | None
) -> np.ndarray:
    """
    Applies a list of transform dicts (translation, rotation, scale, angle_deg)
    to a 3D triangle array of shape (F, 3, 3) or (F, 4, 3).
    """
    if tris.size == 0 or not transforms:
        return tris

    out_tris = tris.copy()
    for t in transforms:
        if not isinstance(t, dict):
            continue

        # 1. Scale
        if "scale" in t:
            sc = t["scale"]
            if any(float(v) != 1.0 for v in sc):
                s_arr = np.array(sc, dtype=np.float64)
                out_tris[:, :3, :] *= s_arr

        # 2. Rotation
        if "angle_deg" in t:
            ang_deg = float(t["angle_deg"])
            if ang_deg != 0.0:
                rad = math.radians(ang_deg)
                cos_a = math.cos(rad)
                sin_a = math.sin(rad)
                center = t.get("center", (0.0, 0.0, 0.0))
                c_arr = np.array(center, dtype=np.float64)[:2]

                pts_2d = out_tris[:, :3, :2] - c_arr
                rot_mat = np.array([[cos_a, -sin_a], [sin_a, cos_a]], dtype=np.float64)
                out_tris[:, :3, :2] = (pts_2d @ rot_mat.T) + c_arr
        elif "rotation" in t:
            rot = t["rotation"]
            if any(float(v) != 0.0 for v in rot):
                from .geometry import euler_rotation_matrix
                R = euler_rotation_matrix(float(rot[0]), float(rot[1]), float(rot[2]))
                center = t.get("center", None)
                if center is not None:
                    c_arr = np.array(center, dtype=np.float64)
                    out_tris[:, :3, :] = (out_tris[:, :3, :] - c_arr) @ R.T + c_arr
                else:
                    out_tris[:, :3, :] = out_tris[:, :3, :] @ R.T

        # 3. Translation / Offset
        tr = t.get("translation", t.get("offset", None))
        if tr is not None and any(float(v) != 0.0 for v in tr):
            t_arr = np.array(tr, dtype=np.float64)
            out_tris[:, :3, :] += t_arr

    return out_tris


def apply_transforms_to_slice_paths(
    paths_obj: Any,
    transforms: list[dict[str, Any]] | None
) -> Any:
    """
    Applies a list of transform dicts (translation, rotation, scale, angle_deg)
    to a SlicePaths object, transforming all coordinates on a copy.
    """
    if paths_obj is None or len(paths_obj) == 0 or not transforms:
        return paths_obj

    new_paths: list[np.ndarray] = []
    for p in paths_obj:
        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 dim == 0:
            new_paths.append(new_p)
            continue

        for t in transforms:
            if not isinstance(t, dict):
                continue

            # 1. Scale
            if "scale" in t:
                sc = t["scale"]
                if any(float(v) != 1.0 for v in sc):
                    s_arr = np.array(sc, dtype=np.float64)
                    if dim >= 3:
                        new_p[:, :3] *= s_arr[:3]
                    else:
                        new_p[:, :dim] *= s_arr[:dim]

            # 2. Rotation
            if "angle_deg" in t:
                ang_deg = float(t["angle_deg"])
                if ang_deg != 0.0:
                    rad = math.radians(ang_deg)
                    cos_a = math.cos(rad)
                    sin_a = math.sin(rad)
                    center = t.get("center", (0.0, 0.0, 0.0))
                    c_arr = np.array(center, dtype=np.float64)[:2]

                    pts_2d = new_p[:, :2] - c_arr
                    rot_mat = np.array([[cos_a, -sin_a], [sin_a, cos_a]], dtype=np.float64)
                    new_p[:, :2] = (pts_2d @ rot_mat.T) + c_arr
            elif "rotation" in t:
                rot = t["rotation"]
                if any(float(v) != 0.0 for v in rot):
                    from .geometry import euler_rotation_matrix
                    R = euler_rotation_matrix(float(rot[0]), float(rot[1]), float(rot[2]))
                    center = t.get("center", None)
                    if center is not None:
                        c_arr = np.array(center, dtype=np.float64)
                        if dim >= 3:
                            new_p[:, :3] = (new_p[:, :3] - c_arr[:3]) @ R.T + c_arr[:3]
                        elif dim == 2:
                            new_p[:, :2] = (new_p[:, :2] - c_arr[:2]) @ R[:2, :2].T + c_arr[:2]
                    else:
                        if dim >= 3:
                            new_p[:, :3] = new_p[:, :3] @ R.T
                        elif dim == 2:
                            new_p[:, :2] = new_p[:, :2] @ R[:2, :2].T

            # 3. Translation / Offset
            tr = t.get("translation", t.get("offset", None))
            if tr is not None and any(float(v) != 0.0 for v in tr):
                t_arr = np.array(tr, dtype=np.float64)
                if dim >= 3:
                    new_p[:, :3] += t_arr[:3]
                else:
                    new_p[:, :dim] += t_arr[:dim]

        new_paths.append(new_p)

    out_sp = paths_obj.copy_with(new_paths)
    res_transforms = [dict(t) for t in getattr(paths_obj, "transforms", [])]
    for dt in transforms:
        if isinstance(dt, dict):
            res_transforms = update_transform(
                res_transforms,
                t_arr=dt.get("translation", dt.get("offset")),
                r_arr=dt.get("rotation"),
                s_arr=dt.get("scale")
            )
            if "angle_deg" in dt and res_transforms:
                res_transforms[0]["angle_deg"] = float(dt["angle_deg"])
                if "center" in dt:
                    res_transforms[0]["center"] = list(dt["center"])
    out_sp.transforms = res_transforms
    return out_sp


def _intersect_segments_2d_py(a: tuple[float, float], b: tuple[float, float], c: tuple[float, float], d: tuple[float, float]) -> float | None:
    denom = (b[0] - a[0]) * (d[1] - c[1]) - (b[1] - a[1]) * (d[0] - c[0])
    if abs(denom) < 1e-12:
        return None
    t = ((c[0] - a[0]) * (d[1] - c[1]) - (c[1] - a[1]) * (d[0] - c[0])) / denom
    u = ((c[0] - a[0]) * (b[1] - a[1]) - (c[1] - a[1]) * (b[0] - a[0])) / denom
    if -1e-6 <= t <= 1.0 + 1e-6 and -1e-6 <= u <= 1.0 + 1e-6:
        return max(0.0, min(1.0, t))
    return None


def _intersect_segment_triangle_py(p0: np.ndarray, p1: np.ndarray, v0: np.ndarray, v1: np.ndarray, v2: np.ndarray) -> tuple[float, np.ndarray] | None:
    e1 = v1 - v0
    e2 = v2 - v0
    n = np.cross(e1, e2)
    n_len_sq = float(np.dot(n, n))
    if n_len_sq < 1e-14:
        return None

    d0 = float(np.dot(p0 - v0, n))
    d1 = float(np.dot(p1 - v0, n))

    if (d0 > 1e-6 and d1 > 1e-6) or (d0 < -1e-6 and d1 < -1e-6):
        return None

    denom = d0 - d1
    if abs(denom) > 1e-12:
        t = d0 / denom
        if t < -1e-6 or t > 1.0 + 1e-6:
            return None
        t = max(0.0, min(1.0, t))
        ipt = p0 + t * (p1 - p0)

        c0 = np.cross(v1 - v0, ipt - v0)
        c1 = np.cross(v2 - v1, ipt - v1)
        c2 = np.cross(v0 - v2, ipt - v2)

        tol = -1e-6 * n_len_sq
        if np.dot(c0, n) >= tol and np.dot(c1, n) >= tol and np.dot(c2, n) >= tol:
            return t, ipt
    elif abs(d0) < 1e-6:
        # Coplanar segment - project along dominant normal axis
        ax0, ax1 = 0, 1
        if abs(n[0]) >= abs(n[1]) and abs(n[0]) >= abs(n[2]):
            ax0, ax1 = 1, 2
        elif abs(n[1]) >= abs(n[0]) and abs(n[1]) >= abs(n[2]):
            ax0, ax1 = 0, 2
        verts = [v0, v1, v2]
        best_t = None
        for e in range(3):
            va = verts[e]
            vb = verts[(e + 1) % 3]
            t_cand = _intersect_segments_2d_py(
                (float(p0[ax0]), float(p0[ax1])), (float(p1[ax0]), float(p1[ax1])),
                (float(va[ax0]), float(va[ax1])), (float(vb[ax0]), float(vb[ax1]))
            )
            if t_cand is not None:
                if best_t is None or t_cand < best_t:
                    best_t = t_cand
        if best_t is not None:
            return best_t, p0 + best_t * (p1 - p0)
    return None


def _slice_triangles_at_z_np(active_tris: np.ndarray, pz: float) -> np.ndarray:
    if len(active_tris) == 0:
        return np.empty((0, 4), dtype=np.float64)

    v0 = active_tris[:, 0, :]
    v1 = active_tris[:, 1, :]
    v2 = active_tris[:, 2, :]

    edges = [(v0, v1), (v1, v2), (v2, v0)]
    pts_list = []
    for va, vb in edges:
        dza = va[:, 2] - pz
        dzb = vb[:, 2] - pz
        cross = (dza * dzb <= 0) & (np.abs(vb[:, 2] - va[:, 2]) > 1e-12)
        denom = np.where(np.abs(vb[:, 2] - va[:, 2]) > 1e-12, vb[:, 2] - va[:, 2], 1.0)
        u = np.clip((pz - va[:, 2]) / denom, 0.0, 1.0)
        qx = va[:, 0] + u * (vb[:, 0] - va[:, 0])
        qy = va[:, 1] + u * (vb[:, 1] - va[:, 1])
        q = np.stack([qx, qy], axis=-1)
        pts_list.append((cross, q))

    c0, q0 = pts_list[0]
    c1, q1 = pts_list[1]
    c2, q2 = pts_list[2]

    segs = []
    for idx in range(len(active_tris)):
        p_cand = []
        if c0[idx]:
            p_cand.append(q0[idx])
        if c1[idx] and (len(p_cand) == 0 or np.linalg.norm(q1[idx] - p_cand[0]) > 1e-6):
            p_cand.append(q1[idx])
        if c2[idx] and (len(p_cand) == 0 or (np.linalg.norm(q2[idx] - p_cand[0]) > 1e-6 and (len(p_cand) < 2 or np.linalg.norm(q2[idx] - p_cand[1]) > 1e-6))):
            p_cand.append(q2[idx])
        if len(p_cand) >= 2:
            segs.append([p_cand[0][0], p_cand[0][1], p_cand[1][0], p_cand[1][1]])

    if not segs:
        return np.empty((0, 4), dtype=np.float64)
    return np.array(segs, dtype=np.float64)


def _intersect_path_2d_segs_np(
    path: np.ndarray,
    mesh_segs: np.ndarray,
    pz: float,
    mode_int: int,
    ref_pt: np.ndarray
) -> tuple[int, float, np.ndarray] | None:
    if len(path) < 4 or len(mesh_segs) == 0:
        return None
    pts = path[:, :3]
    P_start = pts[:-1, :2]  # (M, 2)
    P_end = pts[1:, :2]    # (M, 2)

    Q_start = mesh_segs[:, :2]  # (K, 2)
    Q_end = mesh_segs[:, 2:]    # (K, 2)

    A = P_start[:, np.newaxis, :]  # (M, 1, 2)
    B = P_end[:, np.newaxis, :]    # (M, 1, 2)
    C = Q_start[np.newaxis, :, :]  # (1, K, 2)
    D = Q_end[np.newaxis, :, :]    # (1, K, 2)

    dx_AB = B[..., 0] - A[..., 0]
    dy_AB = B[..., 1] - A[..., 1]
    dx_CD = D[..., 0] - C[..., 0]
    dy_CD = D[..., 1] - C[..., 1]
    dx_AC = C[..., 0] - A[..., 0]
    dy_AC = C[..., 1] - A[..., 1]

    denom = dx_AB * dy_CD - dy_AB * dx_CD
    valid = np.abs(denom) > 1e-12

    denom_safe = np.where(valid, denom, 1.0)
    t_mat = (dx_AC * dy_CD - dy_AC * dx_CD) / denom_safe
    u_mat = (dx_AC * dy_AB - dy_AC * dx_AB) / denom_safe

    hit_mask = valid & (t_mat >= -1e-6) & (t_mat <= 1.0 + 1e-6) & (u_mat >= -1e-6) & (u_mat <= 1.0 + 1e-6)
    hit_m, hit_k = np.where(hit_mask)

    if len(hit_m) > 0:
        t_vals = np.clip(t_mat[hit_m, hit_k], 0.0, 1.0)
        p0_hits = pts[hit_m]
        p1_hits = pts[hit_m + 1]
        ipts = p0_hits + t_vals[:, np.newaxis] * (p1_hits - p0_hits)

        if mode_int == 1:  # FIRST
            return int(hit_m[0]), float(t_vals[0]), ipts[0]

        dists_sq = np.sum((ipts - ref_pt) ** 2, axis=1)
        best_idx = int(np.argmin(dists_sq))
        return int(hit_m[best_idx]), float(t_vals[best_idx]), ipts[best_idx]

    return None


def slice_triangles_at_z(tris: np.ndarray, pz: float) -> np.ndarray:
    """
    Slices a 3D triangle array (F, 3, 3) or (F, 4, 3) at a specific Z height.
    Returns 2D line segments of shape (K, 4) where each segment is [x0, y0, x1, y1].
    """
    if len(tris) == 0:
        return np.empty((0, 4), dtype=np.float64)

    if _c_seam_lib is not None and hasattr(_c_seam_lib, "slice_triangles_at_z_c"):
        try:
            tris_c = np.ascontiguousarray(tris[:, :3, :], dtype=np.float64)
            num_tris = len(tris_c)
            out_segs = np.empty((num_tris, 4), dtype=np.float64)
            num_segs = _c_seam_lib.slice_triangles_at_z_c(
                tris_c.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                num_tris,
                ctypes.c_double(pz),
                out_segs.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
            )
            return out_segs[:num_segs]
        except Exception:
            pass

    return _slice_triangles_at_z_np(tris, pz)


def intersect_path_with_2d_segments(
    path: np.ndarray,
    mesh_segs: np.ndarray,
    pz: float,
    mode: str = "closest",
    mesh_center: np.ndarray | None = None
) -> tuple[int, float, np.ndarray] | None:
    """
    Intersects a closed planar path with pre-computed 2D line segments at height pz.
    Returns (segment_index, t_parameter, intersection_point) or None.
    """
    if len(path) < 4 or len(mesh_segs) == 0:
        return None

    pts = np.ascontiguousarray(path[:, :3], dtype=np.float64)
    segs = np.ascontiguousarray(mesh_segs, dtype=np.float64)
    num_pts = len(pts)
    num_segs = len(segs)

    mode_int = 0
    if mode.lower() == "first":
        mode_int = 1
    elif mode.lower() == "centroid":
        mode_int = 2

    ref_pt = pts[0]
    if mode_int == 2:
        ref_pt = mesh_center if mesh_center is not None else np.array([0.0, 0.0, pz], dtype=np.float64)
    ref_pt_c = np.ascontiguousarray(ref_pt, dtype=np.float64)

    if _c_seam_lib is not None and hasattr(_c_seam_lib, "intersect_path_with_2d_segments_c"):
        try:
            out_seg = ctypes.c_int(-1)
            out_t = ctypes.c_double(0.0)
            out_pt = (ctypes.c_double * 3)()

            found = _c_seam_lib.intersect_path_with_2d_segments_c(
                pts.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                num_pts,
                3,
                segs.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                num_segs,
                ctypes.c_double(pz),
                mode_int,
                ref_pt_c.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                ctypes.byref(out_seg),
                ctypes.byref(out_t),
                out_pt
            )
            if found:
                return int(out_seg.value), float(out_t.value), np.array([out_pt[0], out_pt[1], out_pt[2]], dtype=np.float64)
            return None
        except Exception:
            pass

    return _intersect_path_2d_segs_np(pts, segs, pz, mode_int, ref_pt_c)


def find_path_mesh_intersection(
    path: np.ndarray,
    tri_verts: np.ndarray,
    mode: str = "closest",
    mesh_center: np.ndarray | None = None
) -> tuple[int, float, np.ndarray] | None:
    """
    Finds intersection between a closed polygon path and a mesh of triangles.
    Returns (segment_index, t_parameter, intersection_point) or None.
    """
    if len(path) < 4 or not is_closed_path(path) or len(tri_verts) == 0:
        return None

    pts = np.ascontiguousarray(path[:, :3], dtype=np.float64)
    tris = np.ascontiguousarray(tri_verts[:, :3, :], dtype=np.float64)
    num_pts = len(pts)
    num_tris = len(tris)
    stride = 3

    mode_int = 0
    if mode.lower() == "first":
        mode_int = 1
    elif mode.lower() == "centroid":
        mode_int = 2

    ref_pt = pts[0]
    if mode_int == 2:
        ref_pt = mesh_center if mesh_center is not None else np.mean(tris.reshape(-1, 3), axis=0)
    ref_pt_c = np.ascontiguousarray(ref_pt, dtype=np.float64)

    if _c_seam_lib is not None and hasattr(_c_seam_lib, "find_path_mesh_intersection_c"):
        try:
            out_seg = ctypes.c_int(-1)
            out_t = ctypes.c_double(0.0)
            out_pt = (ctypes.c_double * 3)()

            found = _c_seam_lib.find_path_mesh_intersection_c(
                pts.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                num_pts,
                stride,
                tris.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                num_tris,
                mode_int,
                ref_pt_c.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
                ctypes.byref(out_seg),
                ctypes.byref(out_t),
                out_pt
            )
            if found:
                return int(out_seg.value), float(out_t.value), np.array([out_pt[0], out_pt[1], out_pt[2]], dtype=np.float64)
            return None
        except Exception:
            pass

    # 2. Vectorized / Planar Fallback
    p_min = np.min(pts, axis=0) - 1e-4
    p_max = np.max(pts, axis=0) + 1e-4

    t_min = np.min(tris, axis=1)
    t_max = np.max(tris, axis=1)

    active_mask = np.all(t_max >= p_min, axis=1) & np.all(t_min <= p_max, axis=1)
    if not np.any(active_mask):
        return None

    active_tris = tris[active_mask]
    pz = pts[0, 2]
    is_planar = (p_max[2] - p_min[2]) < 1e-3

    if is_planar:
        segs = slice_triangles_at_z(active_tris, pz)
        if len(segs) > 0:
            return intersect_path_with_2d_segments(path, segs, pz, mode=mode, mesh_center=mesh_center)
        return None

    # 3D fallback for non-planar path
    best_hit = None
    best_dist = 1e30
    M = len(pts) - 1
    for i in range(M):
        p0 = pts[i]
        p1 = pts[i + 1]
        for tr in active_tris:
            hit = _intersect_segment_triangle_py(p0, p1, tr[0], tr[1], tr[2])
            if hit is not None:
                t_cand, ipt = hit
                if mode_int == 1:
                    return i, t_cand, ipt
                d_sq = float(np.sum((ipt - ref_pt) ** 2))
                if d_sq < best_dist:
                    best_dist = d_sq
                    best_hit = (i, t_cand, ipt)
    return best_hit


def reposition_path_seam(
    path: np.ndarray,
    seg_idx: int,
    t: float,
    intersect_pt: np.ndarray,
    extrusion_rate: np.ndarray | None = None,
    feedrate: np.ndarray | None = None,
    extra_attributes: dict[str, np.ndarray] | None = None
) -> tuple[np.ndarray, np.ndarray | None, np.ndarray | None, dict[str, np.ndarray] | None]:
    """
    Reposition the seam of a closed path at the given intersection (seg_idx, t, intersect_pt).
    Properly rotates/interpolates per-vertex attributes (extrusion_rate, feedrate, extra_attributes).
    """
    if len(path) < 4:
        return path, extrusion_rate, feedrate, extra_attributes

    M = len(path) - 1
    k = seg_idx

    # Check snapping
    if t <= 1e-6:
        if k == 0:
            return path, extrusion_rate, feedrate, extra_attributes
        new_p = np.concatenate([path[k:-1], path[:k+1]], axis=0)
        new_ext = np.concatenate([extrusion_rate[k:-1], extrusion_rate[:k+1]]) if extrusion_rate is not None and len(extrusion_rate) == len(path) else extrusion_rate
        new_feed = np.concatenate([feedrate[k:-1], feedrate[:k+1]]) if feedrate is not None and len(feedrate) == len(path) else feedrate
        new_extra = {}
        if extra_attributes:
            for attr_name, attr_arr in extra_attributes.items():
                if len(attr_arr) == len(path):
                    new_extra[attr_name] = np.concatenate([attr_arr[k:-1], attr_arr[:k+1]])
                else:
                    new_extra[attr_name] = attr_arr
        return new_p, new_ext, new_feed, new_extra if extra_attributes else extra_attributes

    if 1.0 - t <= 1e-6:
        next_k = (k + 1) % M
        if next_k == 0:
            return path, extrusion_rate, feedrate, extra_attributes
        new_p = np.concatenate([path[next_k:-1], path[:next_k+1]], axis=0)
        new_ext = np.concatenate([extrusion_rate[next_k:-1], extrusion_rate[:next_k+1]]) if extrusion_rate is not None and len(extrusion_rate) == len(path) else extrusion_rate
        new_feed = np.concatenate([feedrate[next_k:-1], feedrate[:next_k+1]]) if feedrate is not None and len(feedrate) == len(path) else feedrate
        new_extra = {}
        if extra_attributes:
            for attr_name, attr_arr in extra_attributes.items():
                if len(attr_arr) == len(path):
                    new_extra[attr_name] = np.concatenate([attr_arr[next_k:-1], attr_arr[:next_k+1]])
                else:
                    new_extra[attr_name] = attr_arr
        return new_p, new_ext, new_feed, new_extra if extra_attributes else extra_attributes

    # Insert intersection point I
    I_pt = np.array(intersect_pt, dtype=path.dtype)
    if I_pt.shape != path[0].shape:
        I_full = np.zeros_like(path[0])
        I_full[:min(len(I_full), len(I_pt))] = I_pt[:min(len(I_full), len(I_pt))]
        I_pt = I_full

    new_p = np.concatenate([[I_pt], path[k+1:-1], path[:k+1], [I_pt]], axis=0)

    new_ext = None
    if extrusion_rate is not None and len(extrusion_rate) == len(path):
        ext_I = (1.0 - t) * extrusion_rate[k] + t * extrusion_rate[k + 1]
        new_ext = np.concatenate([[ext_I], extrusion_rate[k+1:-1], extrusion_rate[:k+1], [ext_I]])
    else:
        new_ext = extrusion_rate

    new_feed = None
    if feedrate is not None and len(feedrate) == len(path):
        feed_I = (1.0 - t) * feedrate[k] + t * feedrate[k + 1]
        new_feed = np.concatenate([[feed_I], feedrate[k+1:-1], feedrate[:k+1], [feed_I]])
    else:
        new_feed = feedrate

    new_extra = {}
    if extra_attributes:
        for attr_name, attr_arr in extra_attributes.items():
            if len(attr_arr) == len(path):
                attr_I = (1.0 - t) * attr_arr[k] + t * attr_arr[k + 1]
                new_extra[attr_name] = np.concatenate([[attr_I], attr_arr[k+1:-1], attr_arr[:k+1], [attr_I]])
            else:
                new_extra[attr_name] = attr_arr

    return new_p, new_ext, new_feed, new_extra if extra_attributes else extra_attributes


def shift_path_seam_by_distance(
    path: np.ndarray,
    dist_mm: float,
    extrusion_rate: np.ndarray | None = None,
    feedrate: np.ndarray | None = None,
    extra_attributes: dict[str, np.ndarray] | None = None
) -> tuple[np.ndarray, np.ndarray | None, np.ndarray | None, dict[str, np.ndarray] | None]:
    """
    Shifts the starting point (seam) of a closed path along its perimeter by dist_mm.
    Positive dist_mm moves forward along the path direction, negative moves backward.
    Preserves and interpolates per-vertex attributes.
    """
    if len(path) < 4 or not is_closed_path(path) or abs(dist_mm) < 1e-6:
        return path, extrusion_rate, feedrate, extra_attributes

    pts = path[:, :3]
    diffs = np.diff(pts, axis=0)
    seg_lens = np.sqrt(np.sum(diffs ** 2, axis=1))
    total_len = float(np.sum(seg_lens))

    if total_len <= 1e-6:
        return path, extrusion_rate, feedrate, extra_attributes

    target_d = dist_mm % total_len
    if target_d < 0:
        target_d += total_len

    if target_d <= 1e-6 or abs(target_d - total_len) <= 1e-6:
        return path, extrusion_rate, feedrate, extra_attributes

    cum = np.concatenate(([0.0], np.cumsum(seg_lens)))
    seg_idx = int(np.searchsorted(cum, target_d, side='right') - 1)
    seg_idx = max(0, min(len(seg_lens) - 1, seg_idx))

    l_seg = seg_lens[seg_idx]
    if l_seg > 1e-7:
        t = float((target_d - cum[seg_idx]) / l_seg)
    else:
        t = 0.0
    t = max(0.0, min(1.0, t))

    p0 = path[seg_idx]
    p1 = path[seg_idx + 1]
    ipt = (1.0 - t) * p0 + t * p1

    return reposition_path_seam(
        path, seg_idx, t, ipt,
        extrusion_rate=extrusion_rate,
        feedrate=feedrate,
        extra_attributes=extra_attributes
    )


class SlicePaths(list[np.ndarray]):
    is_slice_paths: bool = True
    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."""
        res = self.copy_with([p + offset_vec for p in self])
        res.transforms = update_transform(getattr(self, "transforms", []), t_arr=offset_vec)
        return res

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

    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([])

        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_multiple_slice(p, base_distance, count, width_step, 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_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 accelerated Ramer-Douglas-Peucker (RDP) algorithm,
        preserving closed loop topology, valid polygon constraints, and per-path metadata/attributes.
        """
        if not self or tolerance <= 0:
            return self

        # Use accelerated C implementation if available
        if _c_seam_lib is not None and hasattr(_c_seam_lib, "rdp_simplify_batch_c"):
            try:
                paths_c = [np.ascontiguousarray(p, dtype=np.float64) for p in self]
                num_paths = len(paths_c)

                in_ptrs = (ctypes.POINTER(ctypes.c_double) * num_paths)()
                out_ptrs = (ctypes.POINTER(ctypes.c_double) * num_paths)()
                num_pts = (ctypes.c_int * num_paths)()
                strides = (ctypes.c_int * num_paths)()
                out_counts = (ctypes.c_int * num_paths)()

                out_bufs = [np.empty_like(p) for p in paths_c]
                for i in range(num_paths):
                    in_ptrs[i] = paths_c[i].ctypes.data_as(ctypes.POINTER(ctypes.c_double))
                    out_ptrs[i] = out_bufs[i].ctypes.data_as(ctypes.POINTER(ctypes.c_double))
                    num_pts[i] = len(paths_c[i])
                    strides[i] = paths_c[i].shape[1] if paths_c[i].ndim > 1 else 1

                _c_seam_lib.rdp_simplify_batch_c(
                    in_ptrs,
                    num_pts,
                    strides,
                    num_paths,
                    float(tolerance),
                    1 if preserve_topology else 0,
                    out_ptrs,
                    out_counts
                )

                simplified_paths = []
                for i in range(num_paths):
                    c = out_counts[i]
                    orig_p = self[i]
                    if c >= 2 and (not is_closed_path(orig_p) or c >= (4 if preserve_topology else 3)):
                        simplified_paths.append(out_bufs[i][:c].copy().astype(orig_p.dtype))
                    else:
                        simplified_paths.append(orig_p)

                return self.copy_with(simplified_paths)
            except Exception:
                pass

        # Shapely fallback
        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 next((sp for sp in paths_list if len(sp) > 0), paths_list[0])
        merged = base.copy_with([])

        for sp in paths_list:
            if not sp:
                continue
            merged.extend_with_groups(sp)
        merged.transforms = [dict(t) for t in getattr(base, "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()

        non_empty = [sp for sp in paths_list if len(sp) > 0]
        if not non_empty:
            return SlicePaths()

        base = primary_source if primary_source is not None else non_empty[0]
        
        all_groups = set()
        grouped_inputs = []
        group_z: dict[int, int] = {}
        for p_list in non_empty:
            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()}))
                if g not in group_z and len(p_list[i]) > 0:
                    group_z[g] = int(round(p_list[i][0, 2] * SCALE))
            grouped_inputs.append(d)

        sorted_groups = sorted(list(all_groups), key=lambda g: (group_z.get(g, 0), g))

        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_names = []
        merged_objs = []
        for sp in non_empty:
            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": []})

        merged.transforms = [dict(t) for t in getattr(base, "transforms", [])]
        if merged_objs:
            merged.objects_info = {"objects": merged_objs}
            merged.name = merged_names[0] if merged_names else ""

        return merged

    def get_infill_mask(self) -> SlicePaths:
        """Computes the 2D solid infill mask for each layer, resolving exterior boundaries and interior holes."""
        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.0],
            [s,  c, 0.0],
            [0.0, 0.0, 1.0]
        ], dtype=np.float64)

        new_paths = []
        for p in self:
            # Shift to center, rotate, shift back, preserving exact original Z
            new_p = (p - center) @ rot_mat.T + center
            if len(p) > 0 and p.shape[1] >= 3:
                new_p[:, 2] = p[:, 2]
            new_paths.append(new_p)
            
        res = self.copy_with(new_paths)
        r_arr = np.array([0.0, 0.0, float(angle_rad)], dtype=np.float64)
        res.transforms = update_transform(getattr(self, "transforms", []), r_arr=r_arr)
        if res.transforms:
            res.transforms[0]["angle_deg"] = float(angle_deg)
            if center is not None:
                res.transforms[0]["center"] = center.tolist() if isinstance(center, np.ndarray) else list(center)
        return res

    def get_first_layer_midpoint(self) -> np.ndarray:
        if not self:
            return np.zeros(3)
        
        # Find true minimum Z layer in integer microns across non-empty paths
        min_z_um = min((int(round(p[0, 2] * SCALE)) for p in self if len(p) > 0), default=None)
        if min_z_um is None:
            return np.zeros(3)
        
        first_layer_pts = [p for p in self if len(p) > 0 and int(round(p[0, 2] * SCALE)) == min_z_um]
        if not first_layer_pts:
            return np.zeros(3)
            
        pts = np.concatenate(first_layer_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)

    @classmethod
    def path_operations(cls, paths_list: list[SlicePaths], operation: str = 'ZIP', connect_along_boundary: bool = False, acquire_transforms: bool = False) -> SlicePaths:
        """Executes Path Operations across multiple SlicePaths inputs."""
        return path_operations_slice(paths_list, operation=operation, connect_along_boundary=connect_along_boundary, acquire_transforms=acquire_transforms)

    @classmethod
    def minimum_coverage(cls, paths_list: list[SlicePaths], connect_along_boundary: bool = False) -> SlicePaths:
        """Computes Minimum Coverage across multiple SlicePaths inputs."""
        return minimum_coverage_slice(paths_list, connect_along_boundary=connect_along_boundary)

    @classmethod
    def maximum_coverage(cls, paths_list: list[SlicePaths], connect_along_boundary: bool = False) -> SlicePaths:
        """Computes Maximum Coverage across multiple SlicePaths inputs."""
        return maximum_coverage_slice(paths_list, connect_along_boundary=connect_along_boundary)

    @classmethod
    def intersection(cls, paths_list: list[SlicePaths], connect_along_boundary: bool = False) -> SlicePaths:
        """Computes Intersection across multiple SlicePaths inputs."""
        return intersection_slice(paths_list, connect_along_boundary=connect_along_boundary)

    @classmethod
    def subtract(cls, paths_list: list[SlicePaths], connect_along_boundary: bool = False) -> SlicePaths:
        """Computes Difference (Subtract) across multiple SlicePaths inputs."""
        return subtract_slice(paths_list, connect_along_boundary=connect_along_boundary)

    @classmethod
    def xor(cls, paths_list: list[SlicePaths], connect_along_boundary: bool = False) -> SlicePaths:
        """Computes Exclusive Or (XOR) across multiple SlicePaths inputs."""
        return xor_slice(paths_list, 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 to the back (+Y).
        Preserves all per-path and per-vertex attributes.
        """
        if not self:
            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 is_closed_path(p):
                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

            M = len(p) - 1
            best_idx = find_path_seam(p)
            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])
                continue

            else:
                new_p = np.empty_like(p)
                split = M - best_idx
                new_p[:split] = p[best_idx:-1]
                new_p[split:] = p[:best_idx + 1]
                new_paths.append(new_p)

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

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

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

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

    def stagger_seams(self, stagger_mm: float) -> SlicePaths:
        """
        Shifts the seams of closed paths layer-by-layer, alternating between +stagger_mm (even layers)
        and -stagger_mm (odd layers) along the perimeter to avoid vertical seam alignment weak spots.
        """
        return self.place_seams(mesh_input=None, stagger_mm=stagger_mm)

    def place_seams(self, mesh_input: Any = None, mode: str = "closest", acquire_transforms: bool = True, stagger_mm: float = 0.0) -> SlicePaths:
        """
        Repositions the starting point (seam) of closed slice paths (including outer contours
        and interior holes) to the intersection points between the slice paths and the mesh.
        If acquire_transforms is True, applies upstream transforms from slicepaths metadata to the cutting mesh.
        If stagger_mm is non-zero, offsets the seam along the path perimeter (+stagger_mm on even layers, -stagger_mm on odd layers).
        If there is no intersection for a given path, its seam remains unchanged (aside from any staggering).
        """
        if not self:
            return self

        tris = extract_mesh_triangles(mesh_input) if mesh_input is not None else np.empty((0, 3, 3), dtype=np.float64)
        if len(tris) == 0 and abs(stagger_mm) < 1e-6:
            return self

        if len(tris) > 0:
            if acquire_transforms and getattr(self, "transforms", None):
                mesh_transforms = extract_mesh_transforms(mesh_input)
                delta_transforms = compute_delta_transforms(self.transforms, mesh_transforms)
                if delta_transforms:
                    tris = apply_transforms_to_triangles(tris, delta_transforms)

            mesh_center = np.mean(tris[:, :3, :].reshape(-1, 3), axis=0)

            # Precompute triangle Z bounds
            tz_min = np.min(tris[:, :3, 2], axis=1) - 1e-4
            tz_max = np.max(tris[:, :3, 2], axis=1) + 1e-4
            mesh_z_min = float(np.min(tz_min)) if len(tz_min) > 0 else 0.0
            mesh_z_max = float(np.max(tz_max)) if len(tz_max) > 0 else 0.0
        else:
            mesh_center = np.zeros(3)
            tz_min = np.empty(0, dtype=np.float64)
            tz_max = np.empty(0, dtype=np.float64)
            mesh_z_min = 0.0
            mesh_z_max = 0.0

        n_paths = len(self)
        res_paths: list[np.ndarray | None] = [None] * n_paths
        res_ext: list[np.ndarray | None] = [None] * n_paths
        res_feed: list[np.ndarray | None] = [None] * n_paths
        res_extra: dict[str, list[np.ndarray | None]] = {k: [None] * n_paths for k in self.extra_attributes}

        # Group planar paths by layer integer micron key
        z_to_indices: dict[int, list[int]] = {}
        non_planar_indices: list[int] = []

        for i, p in enumerate(self):
            if len(p) < 4 or not is_closed_path(p):
                res_paths[i] = p
                res_ext[i] = self.extrusion_rate[i]
                res_feed[i] = self.feedrate[i]
                for k in self.extra_attributes:
                    res_extra[k][i] = self.extra_attributes[k][i]
                continue

            pz_min = float(np.min(p[:, 2]))
            pz_max = float(np.max(p[:, 2]))

            if (pz_max - pz_min) < 1e-3:
                z_um = int(round(float(p[0, 2]) * SCALE))
                z_to_indices.setdefault(z_um, []).append(i)
            else:
                non_planar_indices.append(i)

        # Process planar paths layer-by-layer in sorted Z order
        sorted_z_keys = sorted(z_to_indices.keys())
        for layer_idx, z_um in enumerate(sorted_z_keys):
            layer_stagger = float(stagger_mm) if (layer_idx % 2 == 0) else -float(stagger_mm)
            indices = z_to_indices[z_um]
            first_p = self[indices[0]]
            pz = float(first_p[0, 2])

            mesh_segs = np.empty((0, 2, 2), dtype=np.float64)
            if len(tris) > 0 and (mesh_z_min <= pz <= mesh_z_max):
                t_mask = (tz_min <= pz) & (tz_max >= pz)
                active_tris = tris[t_mask]
                if len(active_tris) > 0:
                    mesh_segs = slice_triangles_at_z(active_tris, pz)

            for idx in indices:
                p = self[idx]
                curr_ext = self.extrusion_rate[idx]
                curr_feed = self.feedrate[idx]
                curr_extra = {k: self.extra_attributes[k][idx] for k in self.extra_attributes} if self.extra_attributes else None

                # 1. Reposition to cutting mesh intersection if available
                if len(mesh_segs) > 0:
                    inter = intersect_path_with_2d_segments(p, mesh_segs, pz, mode=mode, mesh_center=mesh_center)
                    if inter is not None:
                        seg_idx, t, ipt = inter
                        p, curr_ext, curr_feed, curr_extra = reposition_path_seam(
                            p, seg_idx, t, ipt,
                            extrusion_rate=curr_ext,
                            feedrate=curr_feed,
                            extra_attributes=curr_extra
                        )

                # 2. Shift seam along path if seam staggering is requested
                if abs(layer_stagger) > 1e-6 and is_closed_path(p):
                    p, curr_ext, curr_feed, curr_extra = shift_path_seam_by_distance(
                        p, layer_stagger,
                        extrusion_rate=curr_ext,
                        feedrate=curr_feed,
                        extra_attributes=curr_extra
                    )

                res_paths[idx] = p
                res_ext[idx] = curr_ext if curr_ext is not None else self.extrusion_rate[idx]
                res_feed[idx] = curr_feed if curr_feed is not None else self.feedrate[idx]
                if curr_extra:
                    for k in self.extra_attributes:
                        res_extra[k][idx] = curr_extra.get(k, self.extra_attributes[k][idx])
                else:
                    for k in self.extra_attributes:
                        res_extra[k][idx] = self.extra_attributes[k][idx]

        # Process any non-planar paths
        for np_idx, idx in enumerate(non_planar_indices):
            layer_idx = len(sorted_z_keys) + np_idx
            layer_stagger = float(stagger_mm) if (layer_idx % 2 == 0) else -float(stagger_mm)
            p = self[idx]
            curr_ext = self.extrusion_rate[idx]
            curr_feed = self.feedrate[idx]
            curr_extra = {k: self.extra_attributes[k][idx] for k in self.extra_attributes} if self.extra_attributes else None

            if len(tris) > 0:
                inter = find_path_mesh_intersection(p, tris, mode=mode, mesh_center=mesh_center)
                if inter is not None:
                    seg_idx, t, ipt = inter
                    p, curr_ext, curr_feed, curr_extra = reposition_path_seam(
                        p, seg_idx, t, ipt,
                        extrusion_rate=curr_ext,
                        feedrate=curr_feed,
                        extra_attributes=curr_extra
                    )

            if abs(layer_stagger) > 1e-6 and is_closed_path(p):
                p, curr_ext, curr_feed, curr_extra = shift_path_seam_by_distance(
                    p, layer_stagger,
                    extrusion_rate=curr_ext,
                    feedrate=curr_feed,
                    extra_attributes=curr_extra
                )

            res_paths[idx] = p
            res_ext[idx] = curr_ext if curr_ext is not None else self.extrusion_rate[idx]
            res_feed[idx] = curr_feed if curr_feed is not None else self.feedrate[idx]
            if curr_extra:
                for k in self.extra_attributes:
                    res_extra[k][idx] = curr_extra.get(k, self.extra_attributes[k][idx])
            else:
                for k in self.extra_attributes:
                    res_extra[k][idx] = self.extra_attributes[k][idx]

        return self.copy_with(
            res_paths,  # type: ignore[arg-type]
            extrusion_rate=res_ext,  # type: ignore[arg-type]
            feedrate=res_feed,  # type: ignore[arg-type]
            extra_attributes=res_extra  # type: ignore[arg-type]
        )


def _slice_single_mesh(mesh_array: np.ndarray, z_values: np.ndarray | list[float], ignore_holes: bool = False) -> list[np.ndarray]:
    """
    Slices a single 3D triangle mesh (shape (F, 4, 3) or (F, 3, 3)) across multiple Z heights.
    Uses manifold3d for robust solid slicing with exact coplanar/topological resolution,
    and falls back to trimesh for open/non-manifold surface meshes.
    """
    if len(mesh_array) == 0:
        return []
    
    z_arr = np.asarray(z_values, dtype=np.float64)
    if len(z_arr) == 0:
        return []

    unique_verts, faces = weld_vertices(mesh_array)
    if len(unique_verts) == 0:
        return []

    all_paths: list[np.ndarray] = []

    # 1. Primary: manifold3d for solid meshes (fastest, guaranteed watertight topological output)
    try:
        mesh_m = manifold3d.Mesh(
            vert_properties=np.ascontiguousarray(unique_verts),
            tri_verts=np.ascontiguousarray(faces)
        )
        m = manifold3d.Manifold(mesh_m)
        # If mesh has inverted normals (negative volume), flip triangle winding
        if m.status() == manifold3d.Error.NoError and m.volume() < 0:
            faces_flipped = faces[:, [0, 2, 1]]
            mesh_m_flipped = manifold3d.Mesh(
                vert_properties=np.ascontiguousarray(unique_verts),
                tri_verts=np.ascontiguousarray(faces_flipped)
            )
            m = manifold3d.Manifold(mesh_m_flipped)

        if m.status() == manifold3d.Error.NoError and not m.is_empty():
            min_z = float(np.min(unique_verts[:, 2]))
            max_z = float(np.max(unique_verts[:, 2]))
            
            for z in z_arr:
                zf = float(z)
                if zf < min_z - 1e-4 or zf > max_z + 1e-4:
                    continue
                # If slicing at boundary face, apply micro-epsilon so boundary is cleanly captured
                slice_z = zf - 1e-5 if abs(zf - max_z) < 1e-5 else (zf + 1e-5 if abs(zf - min_z) < 1e-5 else zf)
                cs = m.slice(slice_z)
                polys = cs.to_polygons()
                for poly in polys:
                    N = len(poly)
                    if N < 3:
                        continue
                    if ignore_holes:
                        if _c_seam_lib is not None and hasattr(_c_seam_lib, "polygon_signed_area_2d_c"):
                            poly_c = np.ascontiguousarray(poly, dtype=np.float64)
                            area = _c_seam_lib.polygon_signed_area_2d_c(poly_c.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), N, 2)
                            if area <= 0:
                                continue
                        else:
                            p_x = poly[:, 0]
                            p_y = poly[:, 1]
                            signed_area = np.sum(p_x[:-1] * p_y[1:] - p_x[1:] * p_y[:-1]) + (p_x[-1] * p_y[0] - p_x[0] * p_y[-1])
                            if signed_area <= 0:
                                continue
                    p3d = np.empty((N + 1, 3), dtype=np.float64)
                    p3d[:N, :2] = poly
                    p3d[:N, 2] = zf
                    p3d[N] = p3d[0]
                    p3d = align_path_seam_sharpest(p3d)
                    all_paths.append(p3d)
            return all_paths
    except Exception:
        pass

    # 2. Fallback: trimesh multiplane slicing for open sheets / non-manifold meshes
    try:
        min_z = float(np.min(unique_verts[:, 2]))
        max_z = float(np.max(unique_verts[:, 2]))
        slice_heights = np.where(
            np.abs(z_arr - max_z) < 1e-5,
            z_arr - 1e-5,
            np.where(np.abs(z_arr - min_z) < 1e-5, z_arr + 1e-5, z_arr)
        )
        tm = trimesh.Trimesh(vertices=unique_verts, faces=faces, process=False)
        lines, _, _ = trimesh.intersections.mesh_multiplane(
            mesh=tm,
            plane_origin=[0.0, 0.0, 0.0],
            plane_normal=[0.0, 0.0, 1.0],
            heights=slice_heights
        )
        for z, line_segs in zip(z_arr, lines):
            if line_segs is None or len(line_segs) == 0:
                continue
            zf = float(z)
            seg2d = np.asarray(line_segs, dtype=np.float64)
            z_col = np.full((len(seg2d), 2, 1), zf, dtype=np.float64)
            seg3d = np.concatenate([seg2d, z_col], axis=-1)
            slice_paths = segments_to_slice(seg3d, ignore_holes=ignore_holes)
            all_paths.extend(slice_paths)
    except Exception:
        pass

    return all_paths


def slice_at_range(models: list[np.ndarray], z_values: np.ndarray | list[float], ignore_holes: bool = False) -> list[list[np.ndarray]]:
    """
    Top-level slicing function that slices multiple model meshes across given Z heights.
    Returns a list of path lists (one per model).
    """
    if not models:
        return []

    mmps: list[list[np.ndarray]] = []
    for model in models:
        model_paths = _slice_single_mesh(model, z_values, ignore_holes=ignore_holes)
        mmps.append(model_paths)
    return mmps

def link_segments(segments: np.ndarray) -> list[np.ndarray]:
    """Links segments (N, 2, 3) into continuous paths using 1-micron integer grid."""
    if len(segments) == 0: return []
    z = segments[0, 0, 2]
    pts = segments[:, :, :2]
    
    # Quantize to integer micron: (x_um, y_um) packed into 64-bit scalar key
    pts_um = np.round(pts * SCALE).astype(np.int64)
    k1_arr = (pts_um[:, 0, 0] << 32) | (pts_um[:, 0, 1] & 0xFFFFFFFF)
    k2_arr = (pts_um[:, 1, 0] << 32) | (pts_um[:, 1, 1] & 0xFFFFFFFF)

    adj: dict[int, list[tuple[int, int]]] = {}
    for i in range(len(segments)):
        k1 = int(k1_arr[i])
        k2 = int(k2_arr[i])
        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
        k1 = int(k1_arr[i])
        k2 = int(k2_arr[i])
        if k1 == k2:
            continue
        
        fwd_pts = [pts[i, 0], pts[i, 1]]
        bwd_pts = []
        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 = int(k1_arr[idx])
                    other_pt = pts[idx, 1] if pk1 == curr_k else pts[idx, 0]
                    fwd_pts.append(other_pt)
                    used.add(idx)
                    curr_k = next_k
                    found = True
                    break
            if not found:
                break
            
        # Backward (O(1) append to bwd_pts)
        curr_k = k1
        while True:
            found = False
            for next_k, idx in adj.get(curr_k, []):
                if idx not in used:
                    pk1 = int(k1_arr[idx])
                    other_pt = pts[idx, 1] if pk1 == curr_k else pts[idx, 0]
                    bwd_pts.append(other_pt)
                    used.add(idx)
                    curr_k = next_k
                    found = True
                    break
            if not found:
                break
            
        all_pts = bwd_pts[::-1] + fwd_pts if bwd_pts else fwd_pts
        p_arr = np.empty((len(all_pts), 3), dtype=np.float64)
        p_arr[:, :2] = all_pts
        p_arr[:, 2] = z

        if len(p_arr) > 2 and is_closed_path(p_arr):
            p_arr[-1] = p_arr[0]
            
        paths.append(p_arr)
    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 is_closed_path(p):
            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).round().astype(np.int64).tolist()
            if abs(pyclipper.Area(scaled)) < 1:
                continue
            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, dtype=np.float64) / 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:
        # Ensure it's closed for logic.
        working_loop = loop
        if not is_closed_path(loop):
            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

def _group_closed_paths_into_shapes(closed_paths: list[np.ndarray]) -> list[list[np.ndarray]]:
    """
    Groups a list of closed paths on a layer into distinct polygon shapes,
    where each shape consists of one outer boundary and zero or more inner holes.
    Uses topological containment with representative points:
    - Paths contained in an even number of enclosing paths (0, 2, ...) are Outer boundaries.
    - Paths contained in an odd number of enclosing paths (1, 3, ...) are Holes belonging
      to the immediately enclosing Outer boundary.
    """
    if not closed_paths:
        return []
    if len(closed_paths) == 1:
        return [[closed_paths[0]]]

    polys: list[tuple[int, np.ndarray, Polygon, float]] = []
    for i, p in enumerate(closed_paths):
        pts = p[:, :2]
        if len(pts) < 3:
            continue
        try:
            poly = Polygon(pts)
            if not poly.is_valid:
                poly = poly.buffer(0)
            if not poly.is_empty:
                polys.append((i, p, poly, float(poly.area)))
        except Exception:
            continue

    if not polys:
        return [[p] for p in closed_paths]
    if len(polys) == 1:
        return [[polys[0][1]]]

    # Sort descending by area so enclosing paths come before enclosed ones
    polys.sort(key=lambda x: x[3], reverse=True)

    # For each polygon, find all strictly larger polygons that contain its representative point
    shapes_map: dict[int, list[np.ndarray]] = {}
    outer_entries: list[tuple[int, Polygon, float]] = []

    for idx, p, poly, area in polys:
        rep_pt = poly.representative_point()
        # Find all larger polygons containing this representative point
        containers = [
            (c_idx, c_poly, c_area)
            for (c_idx, c_p, c_poly, c_area) in polys
            if c_area > area and (c_poly.contains(rep_pt) or c_poly.buffer(1e-5).contains(rep_pt))
        ]

        if len(containers) % 2 == 0:
            # Even enclosing count -> Outer boundary
            shapes_map[idx] = [p]
            outer_entries.append((idx, poly, area))
        else:
            # Odd enclosing count -> Hole belonging to the smallest enclosing outer
            # Find the enclosing outer with the smallest area
            enclosing_outers = [
                (o_idx, o_poly, o_area)
                for (o_idx, o_poly, o_area) in outer_entries
                if o_area > area and (o_poly.contains(rep_pt) or o_poly.buffer(1e-5).contains(rep_pt))
            ]
            if enclosing_outers:
                # Smallest area is the last one because outer_entries is in descending order
                parent_idx = enclosing_outers[-1][0]
                shapes_map[parent_idx].append(p)
            elif outer_entries:
                shapes_map[outer_entries[-1][0]].append(p)
            else:
                shapes_map[idx] = [p]

    return list(shapes_map.values())

@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:
        if len(s) < 2:
            continue
        z = int(round(s[0, 2] * SCALE))
        all_z_keys.add(z)
        subject_dict.setdefault(z, []).append(s)

    for s in object:
        if len(s) < 2:
            continue
        z = int(round(s[0, 2] * SCALE))
        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))
                group += 1
            continue

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

        z = float(k) / SCALE

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

        layer_results = []

        sub_closed = [s for s in subject_slices if is_closed_path(s)]
        sub_open = [s for s in subject_slices if not is_closed_path(s)]
        obj_closed = [s for s in object_slices if is_closed_path(s)]
        obj_open = [s for s in object_slices if not is_closed_path(s)]

        sub_closed_scaled = [(s[:, :2][:-1] * SCALE).round().astype(np.int64).tolist() for s in sub_closed]
        sub_closed_pts = [s[:, :2] for s in sub_closed]
        obj_closed_scaled = [(s[:, :2][:-1] * SCALE).round().astype(np.int64).tolist() for s in obj_closed]
        obj_closed_pts = [s[:, :2] for s in obj_closed]

        def process_open_op(open_paths, clips_scaled, clips_pts, op_type):
            for s in open_paths:
                path_pts = s[:, :2]
                pc.Clear()
                if not connect_along_boundary:
                    pc.AddPath((path_pts * SCALE).round().astype(np.int64).tolist(), pyclipper.PT_SUBJECT, False)
                    for c in clips_scaled:
                        pc.AddPath(c, pyclipper.PT_CLIP, True)
                    solution = pc.Execute2(op_type, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD)
                    for path in pyclipper.OpenPathsFromPolyTree(solution):
                        fpath = np.array(path, dtype=np.float64) / SCALE
                        layer_results.append(np.column_stack([fpath, np.full(len(fpath), z)]))
                else:
                    pc.AddPath((path_pts * SCALE).round().astype(np.int64).tolist(), pyclipper.PT_SUBJECT, False)
                    for c in clips_scaled:
                        pc.AddPath(c, pyclipper.PT_CLIP, True)
                    solution = pc.Execute2(op_type, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD)
                    open_results = pyclipper.OpenPathsFromPolyTree(solution)
                    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, dtype=np.float64) / 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_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':
            if sub_closed and obj_closed:
                pc.Clear()
                for c in sub_closed_scaled:
                    pc.AddPath(c, pyclipper.PT_SUBJECT, True)
                for c in obj_closed_scaled:
                    pc.AddPath(c, pyclipper.PT_CLIP, True)
                solution = pc.Execute2(pyclipper.CT_INTERSECTION, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD)
                for path in pyclipper.ClosedPathsFromPolyTree(solution):
                    fpath = np.array(path, dtype=np.float64) / SCALE
                    fpath_3d = np.column_stack([fpath, np.full(len(fpath), z)])
                    fpath_3d = np.vstack([fpath_3d, fpath_3d[0]])
                    layer_results.append(align_path_seam_sharpest(fpath_3d))
                for path in pyclipper.OpenPathsFromPolyTree(solution):
                    fpath = np.array(path, dtype=np.float64) / SCALE
                    layer_results.append(np.column_stack([fpath, np.full(len(fpath), z)]))

            if sub_open and obj_closed_scaled:
                process_open_op(sub_open, obj_closed_scaled, obj_closed_pts, pyclipper.CT_INTERSECTION)
            if obj_open and sub_closed_scaled:
                process_open_op(obj_open, sub_closed_scaled, sub_closed_pts, pyclipper.CT_INTERSECTION)

        elif operation == 'DIFFERENCE':
            if sub_closed and obj_closed:
                pc.Clear()
                for c in sub_closed_scaled:
                    pc.AddPath(c, pyclipper.PT_SUBJECT, True)
                for c in obj_closed_scaled:
                    pc.AddPath(c, pyclipper.PT_CLIP, True)
                solution = pc.Execute2(pyclipper.CT_DIFFERENCE, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD)
                for path in pyclipper.ClosedPathsFromPolyTree(solution):
                    fpath = np.array(path, dtype=np.float64) / SCALE
                    fpath_3d = np.column_stack([fpath, np.full(len(fpath), z)])
                    fpath_3d = np.vstack([fpath_3d, fpath_3d[0]])
                    layer_results.append(align_path_seam_sharpest(fpath_3d))
                for path in pyclipper.OpenPathsFromPolyTree(solution):
                    fpath = np.array(path, dtype=np.float64) / SCALE
                    layer_results.append(np.column_stack([fpath, np.full(len(fpath), z)]))
            elif sub_closed and not obj_closed:
                layer_results.extend(sub_closed)

            if sub_open and obj_closed_scaled:
                process_open_op(sub_open, obj_closed_scaled, obj_closed_pts, pyclipper.CT_DIFFERENCE)
            elif sub_open and not obj_closed:
                layer_results.extend(sub_open)

        else:
            # UNION, XOR
            pc.Clear()
            for c in sub_closed_scaled:
                pc.AddPath(c, pyclipper.PT_SUBJECT, True)
            for s in sub_open:
                path_pts = s[:, :2]
                scaled = (path_pts * SCALE).round().astype(np.int64).tolist()
                pc.AddPath(scaled, pyclipper.PT_SUBJECT, False)

            for c in obj_closed_scaled:
                pc.AddPath(c, pyclipper.PT_CLIP, True)
            for s in obj_open:
                path_pts = s[:, :2]
                scaled = (path_pts * SCALE).round().astype(np.int64).tolist()
                pc.AddPath(scaled, pyclipper.PT_SUBJECT, False)

            op_type = pyclipper.CT_UNION if operation == 'UNION' else pyclipper.CT_XOR
            solution = pc.Execute2(op_type, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD)
            for path in pyclipper.ClosedPathsFromPolyTree(solution):
                fpath = np.array(path, dtype=np.float64) / SCALE
                fpath_3d = np.column_stack([fpath, np.full(len(fpath), z)])
                fpath_3d = np.vstack([fpath_3d, fpath_3d[0]])
                layer_results.append(align_path_seam_sharpest(fpath_3d))
            for path in pyclipper.OpenPathsFromPolyTree(solution):
                fpath = np.array(path, dtype=np.float64) / 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", "transforms", "objects_info", "name"]:
            setattr(results, attr, copy.deepcopy(val) if isinstance(val, (list, dict)) else val)
    for attr, val in object.__dict__.items():
        if attr not in ["groups", "extrusion_rate", "feedrate", "extra_attributes", "transforms", "objects_info", "name"]:
            if not hasattr(results, attr) or getattr(results, attr) is None or getattr(results, attr) == "" or getattr(results, attr) == []:
                setattr(results, attr, copy.deepcopy(val) if isinstance(val, (list, dict)) else val)

    merged_transforms = []
    merged_names = []
    merged_objs = []
    for sp in [subject, object]:
        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": []})

        subj_trans = getattr(subject, "transforms", None)
        obj_trans = getattr(object, "transforms", None)
        results.transforms = [dict(t) for t in subj_trans] if subj_trans else ([dict(t) for t in obj_trans] if obj_trans else [])
    if merged_objs:
        results.objects_info = {"objects": merged_objs}
        results.name = merged_names[0] if merged_names else ""

    return results

def minimum_coverage_slice(
    paths_list: list[SlicePaths],
    connect_along_boundary: bool = False
) -> SlicePaths:
    """
    Computes the Minimum Coverage across multiple SlicePaths inputs (or within a single SlicePaths input).
    For each layer (grouped by integer micron Z), it calculates the intersection of all objects having paths on that layer.
    In layers where some objects have no paths, the empty paths are ignored (so the result doesn't end up being zero for that layer).
    """
    valid_paths = [sp for sp in paths_list if isinstance(sp, SlicePaths) and len(sp) > 0]
    if not valid_paths:
        return SlicePaths()

    primary = valid_paths[0]

    all_z_keys: set[int] = set()
    input_layer_maps: list[dict[int, list[np.ndarray]]] = []

    for sp in valid_paths:
        layer_map: dict[int, list[np.ndarray]] = {}
        for p in sp:
            if len(p) < 2:
                continue
            z_key = int(round(p[0, 2] * SCALE))
            all_z_keys.add(z_key)
            layer_map.setdefault(z_key, []).append(p)
        input_layer_maps.append(layer_map)

    results = SlicePaths()
    group = 0

    for z_key in sorted(all_z_keys):
        active_indices = [
            i for i, lmap in enumerate(input_layer_maps) if z_key in lmap and len(lmap[z_key]) > 0
        ]

        if not active_indices:
            continue

        if len(active_indices) == 1:
            # Only one object has paths on this layer -> empty paths of other objects ignored
            idx = active_indices[0]
            layer_paths = input_layer_maps[idx][z_key]

            # If this single object has multiple closed paths on this layer (e.g. concentric perimeters):
            closed = [p for p in layer_paths if is_closed_path(p)]
            open_p = [p for p in layer_paths if not is_closed_path(p)]

            if len(closed) > 1 and len(valid_paths) == 1:
                # Intersect concentric perimeters to extract the innermost perimeter
                curr_sp = valid_paths[idx].copy_with([closed[0]])
                for c in closed[1:]:
                    inter_res = boolean_slice(
                        curr_sp,
                        valid_paths[idx].copy_with([c]),
                        operation='INTERSECTION',
                        connect_along_boundary=connect_along_boundary
                    )
                    if inter_res:
                        curr_sp = inter_res
                    else:
                        # If disjoint islands, preserve both!
                        curr_sp.extend([c])
                if open_p:
                    curr_sp.extend(open_p)
                if curr_sp:
                    results.extend(curr_sp)
                    results.groups.extend([group] * len(curr_sp))
                    group += 1
            else:
                results.extend(layer_paths)
                results.groups.extend([group] * len(layer_paths))
                group += 1
            continue

        # Multiple active inputs on this layer -> intersect across active inputs
        idx0 = active_indices[0]
        curr_sp = valid_paths[idx0].copy_with(input_layer_maps[idx0][z_key])

        for next_idx in active_indices[1:]:
            next_sp = valid_paths[next_idx].copy_with(input_layer_maps[next_idx][z_key])
            curr_sp = boolean_slice(
                curr_sp,
                next_sp,
                operation='INTERSECTION',
                connect_along_boundary=connect_along_boundary
            )
            if not curr_sp:
                break

        if curr_sp:
            results.extend(curr_sp)
            results.groups.extend([group] * len(curr_sp))
            group += 1

    calc_ext, calc_feed, calc_extra = _transfer_group_attributes(primary, results, results.groups)
    results.extrusion_rate = calc_ext
    results.feedrate = calc_feed
    results.extra_attributes = calc_extra

    for attr, val in primary.__dict__.items():
        if attr not in ["groups", "extrusion_rate", "feedrate", "extra_attributes", "transforms", "objects_info", "name"]:
            setattr(results, attr, copy.deepcopy(val) if isinstance(val, (list, dict)) else val)

    merged_names = []
    merged_objs = []
    for sp in valid_paths:
        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": []})

    results.transforms = [dict(t) for t in getattr(primary, "transforms", [])]
    if merged_objs:
        results.objects_info = {"objects": merged_objs}
        results.name = merged_names[0] if merged_names else ""

    return results

def maximum_coverage_slice(
    paths_list: list[SlicePaths],
    connect_along_boundary: bool = False
) -> SlicePaths:
    """
    Computes the Maximum Coverage across multiple SlicePaths inputs (or within a single SlicePaths input).
    For each layer (grouped by integer micron Z), it calculates the boolean union of all polygons.
    """
    valid_paths = [sp for sp in paths_list if isinstance(sp, SlicePaths) and len(sp) > 0]
    if not valid_paths:
        return SlicePaths()

    if len(valid_paths) == 1:
        # Single input: if it has multiple perimeters/shapes on a layer, union them
        sp = valid_paths[0]
        layer_map: dict[int, list[np.ndarray]] = {}
        for p in sp:
            if len(p) < 2:
                continue
            z_key = int(round(p[0, 2] * SCALE))
            layer_map.setdefault(z_key, []).append(p)

        results = SlicePaths()
        group = 0
        for z_key in sorted(layer_map.keys()):
            l_paths = layer_map[z_key]
            if len(l_paths) <= 1:
                results.extend(l_paths)
                results.groups.extend([group] * len(l_paths))
                group += 1
            else:
                closed = [p for p in l_paths if is_closed_path(p)]
                open_p = [p for p in l_paths if not is_closed_path(p)]
                if len(closed) <= 1:
                    results.extend(l_paths)
                    results.groups.extend([group] * len(l_paths))
                    group += 1
                else:
                    # Union all closed paths on this layer
                    curr_sp = sp.copy_with([closed[0]])
                    for c in closed[1:]:
                        curr_sp = boolean_slice(
                            curr_sp,
                            sp.copy_with([c]),
                            operation='UNION',
                            connect_along_boundary=connect_along_boundary
                        )
                    if open_p:
                        curr_sp.extend(open_p)
                        curr_sp.groups.extend([group] * len(open_p))
                    if curr_sp:
                        results.extend(curr_sp)
                        results.groups.extend([group] * len(curr_sp))
                        group += 1

        calc_ext, calc_feed, calc_extra = _transfer_group_attributes(sp, results, results.groups)
        results.extrusion_rate = calc_ext
        results.feedrate = calc_feed
        results.extra_attributes = calc_extra
        for attr, val in sp.__dict__.items():
            if attr not in ["groups", "extrusion_rate", "feedrate", "extra_attributes", "transforms", "objects_info", "name"]:
                setattr(results, attr, copy.deepcopy(val) if isinstance(val, (list, dict)) else val)

        results.transforms = [dict(t) for t in getattr(sp, "transforms", [])]
        results.name = getattr(sp, "name", "")
        results.objects_info = copy.deepcopy(getattr(sp, "objects_info", None))
        return results

    # Multiple inputs: sequentially union all inputs with boolean_slice
    accumulated = valid_paths[0]
    for sp in valid_paths[1:]:
        accumulated = boolean_slice(
            accumulated,
            sp,
            operation='UNION',
            connect_along_boundary=connect_along_boundary
        )
    return accumulated

def intersection_slice(
    paths_list: list[SlicePaths],
    connect_along_boundary: bool = False
) -> SlicePaths:
    """Computes the Intersection across multiple SlicePaths inputs (same behavior as Boolean Paths)."""
    if not paths_list:
        return SlicePaths()
    if len(paths_list) == 1:
        return paths_list[0]

    accumulated = paths_list[0]
    for sp in paths_list[1:]:
        if not sp or not accumulated:
            return SlicePaths()
        accumulated = boolean_slice(
            accumulated,
            sp,
            operation='INTERSECTION',
            connect_along_boundary=connect_along_boundary
        )
        if not accumulated:
            return SlicePaths()
    return accumulated

def subtract_slice(
    paths_list: list[SlicePaths],
    connect_along_boundary: bool = False
) -> SlicePaths:
    """Computes the Difference (Subtract) of subsequent SlicePaths inputs from the first input."""
    if not paths_list:
        return SlicePaths()
    if len(paths_list) == 1:
        return paths_list[0]

    accumulated = paths_list[0]
    for sp in paths_list[1:]:
        if not accumulated:
            return SlicePaths()
        if not sp:
            continue
        accumulated = boolean_slice(
            accumulated,
            sp,
            operation='DIFFERENCE',
            connect_along_boundary=connect_along_boundary
        )
    return accumulated

def xor_slice(
    paths_list: list[SlicePaths],
    connect_along_boundary: bool = False
) -> SlicePaths:
    """Computes the Symmetric Difference (XOR) across multiple SlicePaths inputs."""
    if not paths_list:
        return SlicePaths()
    if len(paths_list) == 1:
        return paths_list[0]

    accumulated = paths_list[0]
    for sp in paths_list[1:]:
        if not sp:
            continue
        if not accumulated:
            accumulated = sp
            continue
        accumulated = boolean_slice(
            accumulated,
            sp,
            operation='XOR',
            connect_along_boundary=connect_along_boundary
        )
    return accumulated

def path_operations_slice(
    paths_list: list[SlicePaths],
    operation: str = 'ZIP',
    connect_along_boundary: bool = False,
    acquire_transforms: bool = False
) -> SlicePaths:
    """Dispatches Path Operations across multiple SlicePaths inputs."""
    if not paths_list:
        return SlicePaths()

    if acquire_transforms and len(paths_list) > 1:
        primary = paths_list[0]
        primary_transforms = getattr(primary, "transforms", None)
        if primary_transforms:
            aligned_list = [primary]
            for sec in paths_list[1:]:
                sec_transforms = getattr(sec, "transforms", None)
                delta = compute_delta_transforms(primary_transforms, sec_transforms)
                if delta:
                    aligned_list.append(apply_transforms_to_slice_paths(sec, delta))
                else:
                    aligned_list.append(sec)
            paths_list = aligned_list

    op = str(operation).upper()
    if op in ('ZIP', 'ZIP_PATHS'):
        return SlicePaths.zip_paths(paths_list)
    elif op in ('EXTEND', 'EXTEND_PATHS'):
        return SlicePaths.extend_paths(paths_list)
    elif op in ('MAX_COVERAGE', 'MAXIMUM_COVERAGE', 'UNION'):
        return maximum_coverage_slice(paths_list, connect_along_boundary=connect_along_boundary)
    elif op in ('MIN_COVERAGE', 'MINIMUM_COVERAGE'):
        return minimum_coverage_slice(paths_list, connect_along_boundary=connect_along_boundary)
    elif op in ('INTERSECTION', 'INTERSECT'):
        return intersection_slice(paths_list, connect_along_boundary=connect_along_boundary)
    elif op in ('SUBTRACT', 'DIFFERENCE'):
        return subtract_slice(paths_list, connect_along_boundary=connect_along_boundary)
    elif op == 'XOR':
        return xor_slice(paths_list, connect_along_boundary=connect_along_boundary)
    else:
        return SlicePaths.zip_paths(paths_list)

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.asarray(geom.exterior.coords, dtype=np.float64)
        N = len(ext_coords)
        ext_coords_3d = np.empty((N, 3), dtype=np.float64)
        ext_coords_3d[:, :2] = ext_coords
        ext_coords_3d[:, 2] = z
        paths.append(ext_coords_3d)
        for interior in geom.interiors:
            int_coords = np.asarray(interior.coords, dtype=np.float64)
            M = len(int_coords)
            int_coords_3d = np.empty((M, 3), dtype=np.float64)
            int_coords_3d[:, :2] = int_coords
            int_coords_3d[:, 2] = 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 in integer microns
    groups: dict[int, list[np.ndarray]] = {}
    for p in slice_paths:
        z = p[0, 2]
        z_key = int(round(z * SCALE))
        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 = float(z_key) / SCALE

        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.asarray(geom.coords, dtype=np.float64)
        N = len(coords)
        coords_3d = np.empty((N, 3), dtype=np.float64)
        coords_3d[:, :2] = coords
        coords_3d[:, 2] = z
        paths.append(coords_3d)
    elif isinstance(geom, MultiLineString):
        for line in geom.geoms:
            coords = np.asarray(line.coords, dtype=np.float64)
            N = len(coords)
            coords_3d = np.empty((N, 3), dtype=np.float64)
            coords_3d[:, :2] = coords
            coords_3d[:, 2] = 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 in integer microns
    groups: dict[int, list[np.ndarray]] = {}
    for p in slice_paths:
        z = p[0, 2]
        z_key = int(round(z * SCALE))
        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 = float(z_key) / SCALE

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

            pts_2d = p[:, :2]
            is_closed = is_closed_path(p)
            
            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

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 in 1-micron space."""
    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 is_closed_path(p):
            pts_2d = pts_2d[:-1]
        scaled = (pts_2d * SCALE).round().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:
            if len(node.Contour) >= 3:
                exterior = np.array(node.Contour, dtype=np.float64) / SCALE
                interiors = []
                child_islands = []
                for child in node.Childs:
                    if child.IsHole:
                        if child.Contour and len(child.Contour) >= 3:
                            interiors.append(np.array(child.Contour, dtype=np.float64) / SCALE)
                        for grandchild in child.Childs:
                            child_islands.append(grandchild)
                    else:
                        child_islands.append(child)
                poly = Polygon(exterior, interiors)
                if not poly.is_valid:
                    poly = poly.buffer(0)
                if not poly.is_empty:
                    polygons.append(poly)
                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=64)
def get_infill_mask_slice(slice_paths: list[np.ndarray]) -> tuple[list[np.ndarray], list[int]]:
    if not slice_paths:
        return [], []

    # Group paths by their Z height in integer microns
    groups: dict[int, list[np.ndarray]] = {}
    for p in slice_paths:
        if len(p) == 0:
            continue
        z = p[0, 2]
        z_key = int(round(z * SCALE))
        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 = float(z_key) / SCALE
        
        # Convert layer paths to a single Shapely geometry (handling holes and multiple islands)
        geom = paths_to_shapely_polygon(layer_paths)
        if geom.is_empty:
            continue
            
        # Convert resulting polygon back to 3D paths
        result_paths = geometry_to_paths(geom, z)
        
        for path in result_paths:
            paths_list.append(path)
            groups_list.append(group_id)
            
        group_id += 1
            
    return paths_list, groups_list

@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(),
        key=lambda g: int(round(group_map[g][0][0, 2] * SCALE)) if group_map[g] and len(group_map[g][0]) > 0 else g
    )
    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 = float(int(round(group_map[g][0][0, 2] * SCALE))) / SCALE 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 in integer microns
    groups: dict[int, list[np.ndarray]] = {}
    for p in slice_paths:
        z = p[0, 2]
        z_key = int(round(z * SCALE))
        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=int(round(tolerance * SCALE))) # pyright: ignore[reportAttributeAccessIssue]
    scale_inv = 1.0 / SCALE
    delta_scaled = int(round(distance * SCALE))

    # Buffer each layer independently
    group_id = 0
    for z_key in sorted(groups.keys()):
        layer_paths = groups[z_key]
        z = float(z_key) * scale_inv
        
        pco.Clear()
        scaled_paths = []
        for p in layer_paths:
            if len(p) < 3:
                continue
            path_pts = p[:-1, :2] if is_closed_path(p) else p[:, :2]
            if len(path_pts) < 3:
                continue
            scaled_paths.append(np.rint(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(delta_scaled)
        
        for path in solution:
            if not path:
                continue
            N = len(path)
            fpath_3d = np.empty((N + 1, 3), dtype=np.float64)
            arr = np.asarray(path, dtype=np.float64)
            fpath_3d[:N, :2] = arr * scale_inv
            fpath_3d[:N, 2] = z
            fpath_3d[N] = 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=64)
def buffer_multiple_slice(
    slice_paths: list[np.ndarray],
    base_distance: float,
    count: int,
    width_step: float,
    miter: str = 'round',
    tolerance: float = 0.01
) -> tuple[list[np.ndarray], list[int]]:
    """Generates multiple concentric buffers per layer in a single pass."""
    if not slice_paths or count <= 0:
        return [], []

    groups: dict[int, list[np.ndarray]] = {}
    for p in slice_paths:
        z = p[0, 2]
        z_key = int(round(z * SCALE))
        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=int(round(tolerance * SCALE))) # pyright: ignore[reportAttributeAccessIssue]

    scaled_dists = [int(round((base_distance + i * width_step) * SCALE)) for i in range(count)]
    scale_inv = 1.0 / SCALE

    group_id = 0
    for z_key in sorted(groups.keys()):
        layer_paths = groups[z_key]
        z = float(z_key) * scale_inv

        pco.Clear()
        scaled_paths = []
        for p in layer_paths:
            if len(p) < 3:
                continue
            path_pts = p[:-1, :2] if is_closed_path(p) else p[:, :2]
            if len(path_pts) < 3:
                continue
            scaled_paths.append(np.rint(path_pts * SCALE).astype(np.int64).tolist())

        if not scaled_paths:
            continue

        pco.AddPaths(scaled_paths, jt, pyclipper.ET_CLOSEDPOLYGON) # pyright: ignore[reportAttributeAccessIssue]

        for d_scaled in scaled_dists:
            solution = pco.Execute(d_scaled)
            for path in solution:
                if not path:
                    continue
                N = len(path)
                fpath_3d = np.empty((N + 1, 3), dtype=np.float64)
                arr = np.asarray(path, dtype=np.float64)
                fpath_3d[:N, :2] = arr * scale_inv
                fpath_3d[:N, 2] = z
                fpath_3d[N] = 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 using 1-micron integer grid."""
    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()

    density_um = int(round(density * SCALE))
    if density_um <= 0:
        return [], []

    # Generate the global static grid lines in de-rotated space using integer microns
    rg_min_x_um = int(np.floor(rg_min_x * SCALE)) - density_um
    rg_max_x_um = int(np.ceil(rg_max_x * SCALE)) + density_um
    rg_min_y_um = int(np.floor(rg_min_y * SCALE)) - density_um
    rg_max_y_um = int(np.ceil(rg_max_y * SCALE)) + density_um

    vertical_lines = []
    x_um = rg_min_x_um
    while x_um <= rg_max_x_um:
        vertical_lines.append([
            (x_um, rg_min_y_um),
            (x_um, rg_max_y_um)
        ])
        x_um += density_um

    horizontal_lines = []
    y_um = rg_min_y_um
    while y_um <= rg_max_y_um:
        horizontal_lines.append([
            (rg_min_x_um, y_um),
            (rg_max_x_um, y_um)
        ])
        y_um += density_um

    # Group by Z in integer microns
    groups: dict[int, list[np.ndarray]] = {}
    if hasattr(slice_paths, 'groups') and 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] * SCALE))
            groups.setdefault(z_key, []).append(p)
        
    output_paths: list[np.ndarray] = []
    output_groups: list[int] = []

    # Sort groups by Z in integer microns to determine layer index for alternating rectilinear infill
    unique_groups = sorted(
        groups.keys(),
        key=lambda g: int(round(groups[g][0][0, 2] * SCALE)) if groups[g] and len(groups[g][0]) > 0 else g
    )
    group_to_layer_index = {g: i for i, g in enumerate(unique_groups)}

    for g in unique_groups:
        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 = float(int(round(layer_paths[0][0, 2] * SCALE))) / SCALE if layer_paths and len(layer_paths[0]) > 0 else 0.0
        
        # 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 in integer micron space
        pc = pyclipper.Pyclipper() # pyright: ignore[reportAttributeAccessIssue]
        
        for p in active_paths:
            path_pts = p[:, :2]
            if is_closed_path(p):
                path_pts = path_pts[:-1]
            if len(path_pts) < 3: continue
            scaled = (path_pts * SCALE).round().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, dtype=np.float64) / 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]:
    if len(segments) == 0:
        return {}
    z_um = np.round(segments[:, 0, 2] * SCALE).astype(np.int64)
    unique_z_um = np.unique(z_um)
    return {int(z): segments[z_um == z] for z in unique_z_um}

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

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

        for k, v in kwargs.items():
            setattr(self, k, v)

    @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 t in self.transforms:
                    if isinstance(t, dict):
                        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))
                        if "offset" in t:
                            tr = t["offset"]

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

                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_um: int | None = None
        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]
            z0_um = int(round(z0 * SCALE))

            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 last_z_um is None or z0_um != last_z_um:
                    gcode.append(f";BEFORE_LAYER_CHANGE")
                    gcode.append(f"G0 Z{z0:.3f} F720")
                    gcode.append(f";AFTER_LAYER_CHANGE")
                    last_z_um = z0_um
                    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 last_z_um is None or z0_um != last_z_um:
                    gcode.append(f";BEFORE_LAYER_CHANGE")
                    gcode.append(f"G0 Z{z0:.3f} F720")
                    gcode.append(f";AFTER_LAYER_CHANGE")
                    last_z_um = z0_um
                    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)