nodeslicer/core/slice.py

from __future__ import annotations
from numpy import float64
import numpy as np
import pyclipper
from typing import Any, Literal
from shapely.geometry.base import BaseGeometry

from .geometry import ndarray_lru_cache

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

class SlicePaths(list[np.ndarray]):
    groups: list[int]
    extrusion_rate: list[np.ndarray]
    feedrate: list[np.ndarray]
    extra_attributes: dict[str, list[np.ndarray]]

    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):
        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 {}
    
    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 but preserving all other metadata."""
        if extrusion_rate is None:
            # Only reuse extrusion_rate if path lengths match exactly
            if len(paths) == len(self) and all(len(p) == len(op) for p, op in zip(paths, self)):
                extrusion_rate = self.extrusion_rate
        
        if feedrate is None:
            if len(paths) == len(self) and all(len(p) == len(op) for p, op in zip(paths, self)):
                feedrate = self.feedrate
        
        if extra_attributes is None:
            if len(paths) == len(self) and all(len(p) == len(op) for p, op in zip(paths, self)):
                extra_attributes = self.extra_attributes

        new_obj = SlicePaths(paths, groups if groups is not None else self.groups, extrusion_rate, feedrate, extra_attributes)
        # Copy any other attributes that might have been added
        for k, v in self.__dict__.items():
            if k not in ["groups", "extrusion_rate", "feedrate", "extra_attributes"]:
                setattr(new_obj, k, v)
        return new_obj

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

            zf = float(z)

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

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

        mmps.append(all_model_paths)
    return mmps

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

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

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

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

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

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

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

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

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

    return all_segments

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

@ndarray_lru_cache(maxsize=128)
def segments_to_slice(segments: np.ndarray) -> 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) > 2 and np.allclose(path_pts[0], path_pts[-1]):
            closed_paths.append(p)
        else:
            open_paths.append(p)
            
    res = []
    z = segments[0, 0, 2]
    
    if closed_paths:
        pc = pyclipper.Pyclipper() # pyright: ignore[reportAttributeAccessIssue]
        for p in closed_paths:
            path_pts = p[:, :2]
            path_pts = path_pts[:-1]
            scaled = (path_pts * SCALE).astype(np.int64).tolist()
            pc.AddPath(scaled, pyclipper.PT_SUBJECT, closed = True) # pyright: ignore[reportAttributeAccessIssue]
        
        # Use EVENODD for robust hole detection.
        solution = pc.Execute(pyclipper.CT_UNION, pyclipper.PFT_EVENODD, pyclipper.PFT_EVENODD) # pyright: ignore[reportAttributeAccessIssue]
        
        for path in solution:
            if not path: continue
            fpath = np.array(path) / SCALE
            fpath_3d = np.column_stack([fpath, np.full(len(fpath), z)])
            # Close loop by appending the first point for internal consistency
            fpath_3d = np.vstack([fpath_3d, fpath_3d[0]])
            res.append(fpath_3d)
            
    res.extend(open_paths)
    return res

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

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

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

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

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

    results: SlicePaths = SlicePaths()

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

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

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

        z = subject_slices[0][0, 2]

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

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

        layer_results = []

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

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

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

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

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

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

        results.extend(layer_results)
        results.groups.extend([group] * len(layer_results))
        group += 1
            
    return results

def geometry_to_paths(geom: BaseGeometry, 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 geom.geom_type == 'Polygon':
        ext_coords = np.array(geom.exterior.coords)
        ext_coords_3d = np.column_stack([ext_coords, np.full(len(ext_coords), z)])
        paths.append(ext_coords_3d)
        for interior in geom.interiors:
            int_coords = np.array(interior.coords)
            int_coords_3d = np.column_stack([int_coords, np.full(len(int_coords), z)])
            paths.append(int_coords_3d)
    elif geom.geom_type == 'MultiPolygon':
        for poly in geom.geoms:
            paths.extend(geometry_to_paths(poly, z))
    elif geom.geom_type == 'GeometryCollection':
        for sub_geom in geom.geoms:
            paths.extend(geometry_to_paths(sub_geom, z))
    
    return paths

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

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

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

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

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

            pts_2d = p[:, :2]
            line = LineString(pts_2d)
            
            poly = line.buffer(
                distance,
                cap_style=cap_style,
                join_style=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: BaseGeometry, 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 geom.geom_type in ('LineString', 'LinearRing'):
        coords = np.array(geom.coords)
        coords_3d = np.column_stack([coords, np.full(len(coords), z)])
        paths.append(coords_3d)
    elif geom.geom_type == 'MultiLineString':
        for line in geom.geoms:
            coords = np.array(line.coords)
            coords_3d = np.column_stack([coords, np.full(len(coords), z)])
            paths.append(coords_3d)
    elif geom.geom_type == 'GeometryCollection':
        for sub_geom in geom.geoms:
            paths.extend(offset_geometry_to_paths(sub_geom, z))
    elif geom.geom_type in ('Polygon', 'MultiPolygon'):
        boundary = geom.boundary
        paths.extend(offset_geometry_to_paths(boundary, z))
    return paths

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

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

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

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

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

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

        group_id += 1

    return paths_list, groups_list

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

class Pather(object):

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

    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)

        def segment_length_3d(p1, p2) -> float:
            dx = p2[0] - p1[0]
            dy = p2[1] - p1[1]
            dz = p2[2] - p1[2]
            return (dx * dx + dy * dy + dz * dz) ** 0.5

        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_EXTRUSION_RATE={self.extrusion_rate}",
                f";OFFSET_X={offset_x:.6f}",
                f";OFFSET_Y={offset_y:.6f}",
                "G21 ; set units to mm",
                "G90 ; absolute XYZ positioning",
                "M83 ; relative extrusion mode",
                "G92 E0",
            ]
            gcode += header

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

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

        last_z = -1.0
        max_z = 0
        last_nozzle_z = 0.0

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

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

            x0, y0, z0 = coords[0]

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

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

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

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

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

            old_fr = active_feedrate
            old_x, old_y, old_z = x0, y0, z0
            old_extras = {prefix: None for prefix in extras}
            for j, (p1, p2) in enumerate(zip(coords, coords[1:])):
                seg_len: float = segment_length_3d(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)