from __future__ import annotations
import numpy as np
import pyclipper
from typing import Any

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]

    def __init__(self, paths: list[np.ndarray] = [], groups: list[int] = [], extrusion_rate: 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]
    
    def copy_with(self, paths: list[np.ndarray], groups: list[int] | None = None, extrusion_rate: 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
        
        new_obj = SlicePaths(paths, groups if groups is not None else self.groups, extrusion_rate)
        # Copy any other attributes that might have been added
        for k, v in self.__dict__.items():
            if k not in ["groups", "extrusion_rate"]:
                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 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 boolean(self, other: SlicePaths, operation: str = 'UNION') -> SlicePaths:
        """Performs a 2D boolean operation with another set of paths."""
        res = boolean_slice(self, other, operation)
        # Boolean operations merge paths across layers, so we don't preserve original grouping as precisely.
        return self.copy_with(res)

    def extend_with_groups(self, other: SlicePaths | list[np.ndarray]):
        if isinstance(other, SlicePaths):
            self.extend(other)
            self.groups.extend(other.groups)
        else:
            self.extend(other)
            self.groups.extend(list(range(0,len(other))))

    def append_with_group(self, path: np.ndarray, group: int = 0):
        self.append(path)
        self.groups.append(group)

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 []
    
    pc = pyclipper.Pyclipper() # pyright: ignore[reportAttributeAccessIssue]
    for p in linked:
        # Pyclipper expects open paths (last point != first point) for ET_CLOSEDPOLYGON
        path_pts = p[:, :2]
        if len(path_pts) > 2 and np.allclose(path_pts[0], path_pts[-1]):
            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]
    
    z = segments[0, 0, 2]
    res = []
    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)
    return res

@ndarray_lru_cache(maxsize=128)
def boolean_slice(subject: list[np.ndarray], object: list[np.ndarray], operation: str = 'UNION') -> list[np.ndarray]:
    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: list[np.ndarray] = []

    for k in all_z_keys:
        subject_slices = subject_dict.get(k)
        object_slices = object_dict.get(k)

        if subject_slices is None:
            if object_slices and operation == 'UNION':
                for s in object_slices: results.append(s)                
            continue

        if object_slices is None:
            if operation != 'INTERSECTION':
                for s in subject_slices: results.append(s)
            continue

        z = subject_slices[0][0, 2]

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

        # Subject
        for s in subject_slices:
            path_pts = s[:, :2]
            scaled = (path_pts * SCALE).astype(np.int64).tolist()
            pc.AddPath(scaled, pyclipper.PT_SUBJECT, True) # pyright: ignore[reportAttributeAccessIssue]

        # Object
        for s in object_slices:
            path_pts = s[:, :2]
            scaled = (path_pts * SCALE).astype(np.int64).tolist()
            pc.AddPath(scaled, pyclipper.PT_CLIP, True) # pyright: ignore[reportAttributeAccessIssue]
    
        # Determine operation type
        op_type = pyclipper.CT_UNION # pyright: ignore[reportAttributeAccessIssue]
        if operation == 'UNION': op_type = pyclipper.CT_UNION # pyright: ignore[reportAttributeAccessIssue]
        elif operation == 'INTERSECTION': op_type = pyclipper.CT_INTERSECTION # pyright: ignore[reportAttributeAccessIssue]
        elif operation == 'DIFFERENCE': op_type = pyclipper.CT_DIFFERENCE # pyright: ignore[reportAttributeAccessIssue]
        elif operation == 'XOR': op_type = pyclipper.CT_XOR # pyright: ignore[reportAttributeAccessIssue]

        solution = pc.Execute(op_type, 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)])
            fpath_3d = np.vstack([fpath_3d, fpath_3d[0]])
            results.append(fpath_3d)
            
    return results

@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
        self.extrusion_rates = getattr(paths, "extrusion_rate", [np.ones(len(p)) for p in paths])
        self.bed_center_x: float = 250.
        self.bed_center_y: float = 250.
        self.feedrate: float = 1200.
        self.extrusion_rate: float = 0.017

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

        last_z = -1.0
        max_z = 0

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

            e_rates = self.extrusion_rates[i]

            x0, y0, z0 = coords[0]

            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

            gcode.append(f"G0 X{x0:.3f} Y{y0:.3f} Z{z0:.3f}")
            old_x, old_y, old_z = x0, y0, z0

            fr = self.feedrate
            old_fr = 0
            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_rates[j]

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

                gcode += [cur_gcode]

                max_z = max(max_z, z)
                old_fr = fr

        return "\n".join(gcode)