simplify faster

core/c_seam.c

@@ -804,3 +804,207 @@ EXPORT int shift_path_seam_by_distance_c(
     memcpy(out_pts, in_pts, num_pts * stride * sizeof(double));
     return num_pts;
 }
+
+static inline double point_to_segment_dist_sq(
+    const double* p, const double* a, const double* b, int dim
+) {
+    double ab_sq = 0.0;
+    double ap_dot_ab = 0.0;
+    for (int i = 0; i < dim; ++i) {
+        double d_ab = b[i] - a[i];
+        double d_ap = p[i] - a[i];
+        ab_sq += d_ab * d_ab;
+        ap_dot_ab += d_ap * d_ab;
+    }
+    if (ab_sq <= 1e-14) {
+        double d_sq = 0.0;
+        for (int i = 0; i < dim; ++i) {
+            double d = p[i] - a[i];
+            d_sq += d * d;
+        }
+        return d_sq;
+    }
+    double t = ap_dot_ab / ab_sq;
+    if (t < 0.0) t = 0.0;
+    else if (t > 1.0) t = 1.0;
+    
+    double dist_sq = 0.0;
+    for (int i = 0; i < dim; ++i) {
+        double proj = a[i] + t * (b[i] - a[i]);
+        double d = p[i] - proj;
+        dist_sq += d * d;
+    }
+    return dist_sq;
+}
+
+EXPORT int rdp_simplify_polyline_c(
+    const double* in_pts,
+    int num_pts,
+    int stride,
+    double tol,
+    int preserve_topology,
+    double* out_pts
+) {
+    if (num_pts <= 2 || tol <= 0.0) {
+        memcpy(out_pts, in_pts, num_pts * stride * sizeof(double));
+        return num_pts;
+    }
+
+    double tol_sq = tol * tol;
+    char* keep = (char*)calloc(num_pts, sizeof(char));
+    if (!keep) {
+        memcpy(out_pts, in_pts, num_pts * stride * sizeof(double));
+        return num_pts;
+    }
+
+    keep[0] = 1;
+    keep[num_pts - 1] = 1;
+
+    // Check if closed loop
+    int is_closed = 0;
+    double d_close_sq = 0.0;
+    for (int s = 0; s < 3 && s < stride; ++s) {
+        double d = in_pts[s] - in_pts[(num_pts - 1) * stride + s];
+        d_close_sq += d * d;
+    }
+    if (d_close_sq < 1e-12) {
+        is_closed = 1;
+    }
+
+    // Stack for iterative RDP: pairs of (start_idx, end_idx)
+    int stack_cap = 256;
+    int stack_size = 0;
+    int stack_buf[512];
+    int* stack = stack_buf;
+    int heap_alloc = 0;
+
+    if (is_closed) {
+        double max_d_sq = -1.0;
+        int max_idx = -1;
+        const double* p0 = in_pts;
+        for (int i = 1; i < num_pts - 1; ++i) {
+            double d_sq = 0.0;
+            for (int s = 0; s < 3 && s < stride; ++s) {
+                double diff = in_pts[i * stride + s] - p0[s];
+                d_sq += diff * diff;
+            }
+            if (d_sq > max_d_sq) {
+                max_d_sq = d_sq;
+                max_idx = i;
+            }
+        }
+        if (max_idx > 0 && max_d_sq > tol_sq) {
+            keep[max_idx] = 1;
+            stack[stack_size * 2 + 0] = 0;
+            stack[stack_size * 2 + 1] = max_idx;
+            stack_size++;
+            stack[stack_size * 2 + 0] = max_idx;
+            stack[stack_size * 2 + 1] = num_pts - 1;
+            stack_size++;
+        }
+    } else {
+        stack[stack_size * 2 + 0] = 0;
+        stack[stack_size * 2 + 1] = num_pts - 1;
+        stack_size++;
+    }
+
+    while (stack_size > 0) {
+        stack_size--;
+        int start = stack[stack_size * 2 + 0];
+        int end = stack[stack_size * 2 + 1];
+
+        if (end - start <= 1) continue;
+
+        const double* a = in_pts + (start * stride);
+        const double* b = in_pts + (end * stride);
+
+        double max_d_sq = 0.0;
+        int max_idx = -1;
+
+        for (int i = start + 1; i < end; ++i) {
+            const double* p = in_pts + (i * stride);
+            double d_sq = point_to_segment_dist_sq(p, a, b, (stride >= 3 ? 3 : stride));
+            if (d_sq > max_d_sq) {
+                max_d_sq = d_sq;
+                max_idx = i;
+            }
+        }
+
+        if (max_d_sq > tol_sq && max_idx > 0) {
+            keep[max_idx] = 1;
+            if (stack_size + 2 >= stack_cap) {
+                stack_cap *= 2;
+                if (!heap_alloc) {
+                    stack = (int*)malloc(stack_cap * 2 * sizeof(int));
+                    if (stack) {
+                        memcpy(stack, stack_buf, stack_size * 2 * sizeof(int));
+                        heap_alloc = 1;
+                    } else {
+                        stack = stack_buf;
+                    }
+                } else {
+                    int* new_st = (int*)realloc(stack, stack_cap * 2 * sizeof(int));
+                    if (new_st) stack = new_st;
+                }
+            }
+            if (stack_size + 2 < stack_cap) {
+                stack[stack_size * 2 + 0] = start;
+                stack[stack_size * 2 + 1] = max_idx;
+                stack_size++;
+                stack[stack_size * 2 + 0] = max_idx;
+                stack[stack_size * 2 + 1] = end;
+                stack_size++;
+            }
+        }
+    }
+
+    if (heap_alloc) {
+        free(stack);
+    }
+
+    if (preserve_topology && is_closed) {
+        int kept_count = 0;
+        for (int i = 0; i < num_pts; ++i) if (keep[i]) kept_count++;
+        if (kept_count < 4 && num_pts >= 4) {
+            keep[num_pts / 3] = 1;
+            keep[(2 * num_pts) / 3] = 1;
+        }
+    }
+
+    int out_count = 0;
+    for (int i = 0; i < num_pts; ++i) {
+        if (keep[i]) {
+            memcpy(out_pts + (out_count * stride), in_pts + (i * stride), stride * sizeof(double));
+            out_count++;
+        }
+    }
+
+    free(keep);
+    return out_count;
+}
+
+EXPORT void rdp_simplify_batch_c(
+    const double** in_pts_arr,
+    const int* num_pts_arr,
+    const int* stride_arr,
+    int num_paths,
+    double tol,
+    int preserve_topology,
+    double** out_pts_arr,
+    int* out_counts_arr
+) {
+    #ifdef _OPENMP
+    #pragma omp parallel for schedule(static)
+    #endif
+    for (int i = 0; i < num_paths; ++i) {
+        out_counts_arr[i] = rdp_simplify_polyline_c(
+            in_pts_arr[i],
+            num_pts_arr[i],
+            stride_arr[i],
+            tol,
+            preserve_topology,
+            out_pts_arr[i]
+        );
+    }
+}
+

core/slice.py

@@ -120,12 +120,15 @@ try:
     if rebuild:
         import subprocess
         for compiler in ("gcc", "clang"):
-            try:
-                subprocess.run([compiler, "-O3", "-fPIC", "-shared", _c_src, "-o", _so_path, "-lm"], check=True, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
-                if os.path.exists(_so_path):
-                    break
-            except Exception:
-                pass
+            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 = [
@@ -214,6 +217,30 @@ try:
                 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
 
@@ -1178,10 +1205,57 @@ class SlicePaths(list[np.ndarray]):
         return new_paths
 
     def simplify(self, tolerance: float = 0.1, preserve_topology: bool = True) -> SlicePaths:
-        """Simplifies paths using Shapely while conserving all per-path metadata and attributes."""
+        """
+        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 = []
@@ -1504,8 +1578,6 @@ class SlicePaths(list[np.ndarray]):
         """
         return self.place_seams(mesh_input=None, stagger_mm=stagger_mm)
 
-    from ..infra.debug import debug
-    @debug
     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

scratch/benchmark_simplify.py

@@ -0,0 +1,71 @@
+import sys
+from pathlib import Path
+from unittest.mock import MagicMock
+
+PROJECT_ROOT = Path(__file__).resolve().parent.parent
+sys.path.append(str(PROJECT_ROOT.parent))
+
+class DummyNodeSocket: pass
+class DummyNode: pass
+class DummyObject: pass
+class DummyMesh: pass
+class DummyNodeSocketShader: pass
+class DummyNodeCustomGroup: pass
+class DummyNodeTree: pass
+class DummyOperator: pass
+class DummyMenu: pass
+class DummyPanel: pass
+class DummyPropertyGroup: pass
+
+mock_types = MagicMock()
+mock_types.NodeSocket = DummyNodeSocket
+mock_types.Node = DummyNode
+mock_types.Object = DummyObject
+mock_types.Mesh = DummyMesh
+mock_types.NodeSocketShader = DummyNodeSocketShader
+mock_types.NodeCustomGroup = DummyNodeCustomGroup
+mock_types.NodeTree = DummyNodeTree
+mock_types.Operator = DummyOperator
+mock_types.Menu = DummyMenu
+mock_types.Panel = DummyPanel
+mock_types.PropertyGroup = DummyPropertyGroup
+
+mock_bpy = MagicMock()
+mock_bpy.context.scene.unit_settings.scale_length = 0.001
+mock_bpy.types = mock_types
+
+sys.modules['bpy'] = mock_bpy
+sys.modules['bpy.types'] = mock_types
+sys.modules['bpy.props'] = MagicMock()
+sys.modules['bpy.utils'] = MagicMock()
+sys.modules['bpy.utils.previews'] = MagicMock()
+sys.modules['bpy.app'] = MagicMock()
+sys.modules['bpy.app.handlers'] = MagicMock()
+
+import numpy as np
+import time
+from NodeSlicer.core.slice import SlicePaths
+
+n_layers = 100
+paths_per_layer = 5
+pts_per_path = 300
+
+paths = []
+np.random.seed(42)
+for z in np.linspace(0, 20, n_layers):
+    for r in [5.0, 10.0, 15.0, 20.0, 25.0]:
+        theta = np.linspace(0, 2*np.pi, pts_per_path)
+        x = r * np.cos(theta) + np.random.normal(0, 0.005, pts_per_path)
+        y = r * np.sin(theta) + np.random.normal(0, 0.005, pts_per_path)
+        x[-1] = x[0]
+        y[-1] = y[0]
+        p = np.column_stack([x, y, np.full(pts_per_path, z)])
+        paths.append(p)
+
+sp = SlicePaths(paths)
+print(f"Total paths: {len(sp)}, Total points: {sum(len(p) for p in sp)}")
+
+t0 = time.perf_counter()
+res_shapely = sp.simplify(0.05, preserve_topology=True)
+t1 = time.perf_counter()
+print(f"Current Shapely simplify: {t1 - t0:.4f} s, output points: {sum(len(p) for p in res_shapely)}")

scratch/test_batch_c_rdp.py

@@ -0,0 +1,256 @@
+import numpy as np
+import time
+import ctypes
+import subprocess, tempfile, os, sys
+
+c_code = """
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+static inline double point_to_segment_dist_sq(
+    const double* p, const double* a, const double* b, int dim
+) {
+    double ab_sq = 0.0;
+    double ap_dot_ab = 0.0;
+    for (int i = 0; i < dim; ++i) {
+        double d_ab = b[i] - a[i];
+        double d_ap = p[i] - a[i];
+        ab_sq += d_ab * d_ab;
+        ap_dot_ab += d_ap * d_ab;
+    }
+    if (ab_sq <= 1e-14) {
+        double d_sq = 0.0;
+        for (int i = 0; i < dim; ++i) {
+            double d = p[i] - a[i];
+            d_sq += d * d;
+        }
+        return d_sq;
+    }
+    double t = ap_dot_ab / ab_sq;
+    if (t < 0.0) t = 0.0;
+    else if (t > 1.0) t = 1.0;
+    
+    double dist_sq = 0.0;
+    for (int i = 0; i < dim; ++i) {
+        double proj = a[i] + t * (b[i] - a[i]);
+        double d = p[i] - proj;
+        dist_sq += d * d;
+    }
+    return dist_sq;
+}
+
+int rdp_simplify_polyline_c(
+    const double* in_pts,
+    int num_pts,
+    int stride,
+    double tol,
+    int preserve_topology,
+    double* out_pts
+) {
+    if (num_pts <= 2 || tol <= 0.0) {
+        memcpy(out_pts, in_pts, num_pts * stride * sizeof(double));
+        return num_pts;
+    }
+
+    double tol_sq = tol * tol;
+    char* keep = (char*)calloc(num_pts, sizeof(char));
+    if (!keep) {
+        memcpy(out_pts, in_pts, num_pts * stride * sizeof(double));
+        return num_pts;
+    }
+
+    keep[0] = 1;
+    keep[num_pts - 1] = 1;
+
+    // Check if closed loop
+    int is_closed = 0;
+    double d_close_sq = 0.0;
+    for (int s = 0; s < 3 && s < stride; ++s) {
+        double d = in_pts[s] - in_pts[(num_pts - 1) * stride + s];
+        d_close_sq += d * d;
+    }
+    if (d_close_sq < 1e-12) {
+        is_closed = 1;
+    }
+
+    // Stack for iterative RDP: pairs of (start_idx, end_idx)
+    int stack_cap = 256;
+    int stack_size = 0;
+    int stack_buf[512];
+    int* stack = stack_buf;
+    int heap_alloc = 0;
+
+    if (is_closed) {
+        double max_d_sq = -1.0;
+        int max_idx = -1;
+        const double* p0 = in_pts;
+        for (int i = 1; i < num_pts - 1; ++i) {
+            double d_sq = 0.0;
+            for (int s = 0; s < 3 && s < stride; ++s) {
+                double diff = in_pts[i * stride + s] - p0[s];
+                d_sq += diff * diff;
+            }
+            if (d_sq > max_d_sq) {
+                max_d_sq = d_sq;
+                max_idx = i;
+            }
+        }
+        if (max_idx > 0 && max_d_sq > tol_sq) {
+            keep[max_idx] = 1;
+            stack[stack_size * 2 + 0] = 0;
+            stack[stack_size * 2 + 1] = max_idx;
+            stack_size++;
+            stack[stack_size * 2 + 0] = max_idx;
+            stack[stack_size * 2 + 1] = num_pts - 1;
+            stack_size++;
+        }
+    } else {
+        stack[stack_size * 2 + 0] = 0;
+        stack[stack_size * 2 + 1] = num_pts - 1;
+        stack_size++;
+    }
+
+    while (stack_size > 0) {
+        stack_size--;
+        int start = stack[stack_size * 2 + 0];
+        int end = stack[stack_size * 2 + 1];
+
+        if (end - start <= 1) continue;
+
+        const double* a = in_pts + (start * stride);
+        const double* b = in_pts + (end * stride);
+
+        double max_d_sq = 0.0;
+        int max_idx = -1;
+
+        for (int i = start + 1; i < end; ++i) {
+            const double* p = in_pts + (i * stride);
+            double d_sq = point_to_segment_dist_sq(p, a, b, (stride >= 3 ? 3 : stride));
+            if (d_sq > max_d_sq) {
+                max_d_sq = d_sq;
+                max_idx = i;
+            }
+        }
+
+        if (max_d_sq > tol_sq && max_idx > 0) {
+            keep[max_idx] = 1;
+            if (stack_size + 2 >= stack_cap) {
+                stack_cap *= 2;
+                if (!heap_alloc) {
+                    stack = (int*)malloc(stack_cap * 2 * sizeof(int));
+                    memcpy(stack, stack_buf, stack_size * 2 * sizeof(int));
+                    heap_alloc = 1;
+                } else {
+                    stack = (int*)realloc(stack, stack_cap * 2 * sizeof(int));
+                }
+            }
+            stack[stack_size * 2 + 0] = start;
+            stack[stack_size * 2 + 1] = max_idx;
+            stack_size++;
+            stack[stack_size * 2 + 0] = max_idx;
+            stack[stack_size * 2 + 1] = end;
+            stack_size++;
+        }
+    }
+
+    if (heap_alloc) {
+        free(stack);
+    }
+
+    if (preserve_topology && is_closed) {
+        int kept_count = 0;
+        for (int i = 0; i < num_pts; ++i) if (keep[i]) kept_count++;
+        if (kept_count < 4 && num_pts >= 4) {
+            keep[num_pts / 3] = 1;
+            keep[(2 * num_pts) / 3] = 1;
+        }
+    }
+
+    int out_count = 0;
+    for (int i = 0; i < num_pts; ++i) {
+        if (keep[i]) {
+            memcpy(out_pts + (out_count * stride), in_pts + (i * stride), stride * sizeof(double));
+            out_count++;
+        }
+    }
+
+    free(keep);
+    return out_count;
+}
+
+// Batch simplify with OpenMP parallelization
+void rdp_simplify_batch_c(
+    const double** in_pts_arr,
+    const int* num_pts_arr,
+    const int* stride_arr,
+    int num_paths,
+    double tol,
+    int preserve_topology,
+    double** out_pts_arr,
+    int* out_counts_arr
+) {
+    #pragma omp parallel for schedule(static)
+    for (int i = 0; i < num_paths; ++i) {
+        out_counts_arr[i] = rdp_simplify_polyline_c(
+            in_pts_arr[i],
+            num_pts_arr[i],
+            stride_arr[i],
+            tol,
+            preserve_topology,
+            out_pts_arr[i]
+        );
+    }
+}
+"""
+
+with tempfile.NamedTemporaryFile("w", suffix=".c", delete=False) as f:
+    f.write(c_code)
+    c_path = f.name
+
+so_path = c_path.replace(".c", ".so")
+subprocess.run(["gcc", "-O3", "-fopenmp", "-fPIC", "-shared", c_path, "-o", so_path], check=True)
+lib = ctypes.CDLL(so_path)
+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)
+]
+
+sys.path.append(os.path.dirname(__file__))
+from benchmark_simplify import sp
+
+paths_c = [np.ascontiguousarray(p, dtype=np.float64) for p in sp]
+num_paths = len(paths_c)
+
+# Setup pointers
+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]
+
+t0 = time.perf_counter()
+lib.rdp_simplify_batch_c(
+    in_ptrs, num_pts, strides, num_paths, 0.05, 1, out_ptrs, out_counts
+)
+res = [out_bufs[i][:out_counts[i]].copy() for i in range(num_paths)]
+t1 = time.perf_counter()
+
+print(f"OpenMP Batch C RDP simplify: {t1 - t0:.5f} s, output points: {sum(len(p) for p in res)}")

scratch/test_c_rdp.py

@@ -0,0 +1,227 @@
+import numpy as np
+import time
+import ctypes
+
+c_code = """
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+static inline double point_to_segment_dist_sq(
+    const double* p, const double* a, const double* b, int dim
+) {
+    double ab_sq = 0.0;
+    double ap_dot_ab = 0.0;
+    for (int i = 0; i < dim; ++i) {
+        double d_ab = b[i] - a[i];
+        double d_ap = p[i] - a[i];
+        ab_sq += d_ab * d_ab;
+        ap_dot_ab += d_ap * d_ab;
+    }
+    if (ab_sq <= 1e-14) {
+        double d_sq = 0.0;
+        for (int i = 0; i < dim; ++i) {
+            double d = p[i] - a[i];
+            d_sq += d * d;
+        }
+        return d_sq;
+    }
+    double t = ap_dot_ab / ab_sq;
+    if (t < 0.0) t = 0.0;
+    else if (t > 1.0) t = 1.0;
+    
+    double dist_sq = 0.0;
+    for (int i = 0; i < dim; ++i) {
+        double proj = a[i] + t * (b[i] - a[i]);
+        double d = p[i] - proj;
+        dist_sq += d * d;
+    }
+    return dist_sq;
+}
+
+int rdp_simplify_polyline_c(
+    const double* in_pts,
+    int num_pts,
+    int stride,
+    double tol,
+    int preserve_topology,
+    double* out_pts
+) {
+    if (num_pts <= 2 || tol <= 0.0) {
+        memcpy(out_pts, in_pts, num_pts * stride * sizeof(double));
+        return num_pts;
+    }
+
+    double tol_sq = tol * tol;
+    char* keep = (char*)calloc(num_pts, sizeof(char));
+    if (!keep) {
+        memcpy(out_pts, in_pts, num_pts * stride * sizeof(double));
+        return num_pts;
+    }
+
+    keep[0] = 1;
+    keep[num_pts - 1] = 1;
+
+    // Check if closed loop
+    int is_closed = 0;
+    double d_close_sq = 0.0;
+    for (int s = 0; s < 3 && s < stride; ++s) {
+        double d = in_pts[s] - in_pts[(num_pts - 1) * stride + s];
+        d_close_sq += d * d;
+    }
+    if (d_close_sq < 1e-12) {
+        is_closed = 1;
+    }
+
+    // Stack for iterative RDP: pairs of (start_idx, end_idx)
+    int stack_cap = 256;
+    int stack_size = 0;
+    int* stack = (int*)malloc(stack_cap * 2 * sizeof(int));
+
+    if (is_closed) {
+        // For closed loop, find point furthest from start point
+        double max_d_sq = -1.0;
+        int max_idx = -1;
+        const double* p0 = in_pts;
+        for (int i = 1; i < num_pts - 1; ++i) {
+            double d_sq = 0.0;
+            for (int s = 0; s < 3 && s < stride; ++s) {
+                double diff = in_pts[i * stride + s] - p0[s];
+                d_sq += diff * diff;
+            }
+            if (d_sq > max_d_sq) {
+                max_d_sq = d_sq;
+                max_idx = i;
+            }
+        }
+        if (max_idx > 0 && max_d_sq > tol_sq) {
+            keep[max_idx] = 1;
+            stack[stack_size * 2 + 0] = 0;
+            stack[stack_size * 2 + 1] = max_idx;
+            stack_size++;
+            stack[stack_size * 2 + 0] = max_idx;
+            stack[stack_size * 2 + 1] = num_pts - 1;
+            stack_size++;
+        }
+    } else {
+        stack[stack_size * 2 + 0] = 0;
+        stack[stack_size * 2 + 1] = num_pts - 1;
+        stack_size++;
+    }
+
+    while (stack_size > 0) {
+        stack_size--;
+        int start = stack[stack_size * 2 + 0];
+        int end = stack[stack_size * 2 + 1];
+
+        if (end - start <= 1) continue;
+
+        const double* a = in_pts + (start * stride);
+        const double* b = in_pts + (end * stride);
+
+        double max_d_sq = 0.0;
+        int max_idx = -1;
+
+        for (int i = start + 1; i < end; ++i) {
+            const double* p = in_pts + (i * stride);
+            double d_sq = point_to_segment_dist_sq(p, a, b, (stride >= 3 ? 3 : stride));
+            if (d_sq > max_d_sq) {
+                max_d_sq = d_sq;
+                max_idx = i;
+            }
+        }
+
+        if (max_d_sq > tol_sq && max_idx > 0) {
+            keep[max_idx] = 1;
+            if (stack_size + 2 >= stack_cap) {
+                stack_cap *= 2;
+                stack = (int*)realloc(stack, stack_cap * 2 * sizeof(int));
+            }
+            stack[stack_size * 2 + 0] = start;
+            stack[stack_size * 2 + 1] = max_idx;
+            stack_size++;
+            stack[stack_size * 2 + 0] = max_idx;
+            stack[stack_size * 2 + 1] = end;
+            stack_size++;
+        }
+    }
+
+    free(stack);
+
+    // If preserve_topology is requested on closed loop and fewer than 4 points kept:
+    if (preserve_topology && is_closed) {
+        int kept_count = 0;
+        for (int i = 0; i < num_pts; ++i) if (keep[i]) kept_count++;
+        if (kept_count < 4 && num_pts >= 4) {
+            // Keep at least 3 distinct vertices (e.g. 0, num_pts/3, 2*num_pts/3, num_pts-1)
+            keep[num_pts / 3] = 1;
+            keep[(2 * num_pts) / 3] = 1;
+        }
+    }
+
+    // Collect output points
+    int out_count = 0;
+    for (int i = 0; i < num_pts; ++i) {
+        if (keep[i]) {
+            memcpy(out_pts + (out_count * stride), in_pts + (i * stride), stride * sizeof(double));
+            out_count++;
+        }
+    }
+
+    free(keep);
+    return out_count;
+}
+"""
+
+import subprocess, tempfile, os
+
+with tempfile.NamedTemporaryFile("w", suffix=".c", delete=False) as f:
+    f.write(c_code)
+    c_path = f.name
+
+so_path = c_path.replace(".c", ".so")
+subprocess.run(["gcc", "-O3", "-fPIC", "-shared", c_path, "-o", so_path], check=True)
+lib = ctypes.CDLL(so_path)
+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)
+]
+lib.rdp_simplify_polyline_c.restype = ctypes.c_int
+
+def simplify_c(p: np.ndarray, tol: float, preserve_topology: bool = True) -> np.ndarray:
+    if len(p) <= 2 or tol <= 0:
+        return p
+    p_c = np.ascontiguousarray(p, dtype=np.float64)
+    out_buf = np.empty_like(p_c)
+    n_out = lib.rdp_simplify_polyline_c(
+        p_c.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+        len(p_c),
+        p_c.shape[1],
+        float(tol),
+        1 if preserve_topology else 0,
+        out_buf.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
+    )
+    return out_buf[:n_out].copy()
+
+# Test benchmark on dataset
+import sys
+sys.path.append(os.path.dirname(__file__))
+from benchmark_simplify import sp
+
+print(f"Testing on {len(sp)} paths...")
+t0 = time.perf_counter()
+res_c = [simplify_c(p, 0.05, preserve_topology=True) for p in sp]
+t1 = time.perf_counter()
+print(f"C RDP simplify: {t1 - t0:.4f} s, output points: {sum(len(p) for p in res_c)}")
+
+# Check accuracy/geometry similarity
+from shapely.geometry import LineString
+p_sample = sp[0]
+p_shapely = np.array(LineString(p_sample).simplify(0.05, preserve_topology=True).coords)
+p_c_sample = simplify_c(p_sample, 0.05, preserve_topology=True)
+
+print(f"Sample path: original {len(p_sample)} pts -> Shapely: {len(p_shapely)} pts, C: {len(p_c_sample)} pts")

tests/test_simplify.py

@@ -0,0 +1,144 @@
+import unittest
+import numpy as np
+import sys
+from pathlib import Path
+from unittest.mock import MagicMock
+
+PROJECT_ROOT = Path(__file__).resolve().parent.parent
+sys.path.append(str(PROJECT_ROOT.parent))
+
+# Setup dummy bpy
+class DummyNodeSocket: pass
+class DummyNode: pass
+class DummyObject: pass
+class DummyMesh: pass
+class DummyNodeSocketShader: pass
+class DummyNodeCustomGroup: pass
+class DummyNodeTree: pass
+class DummyOperator: pass
+class DummyMenu: pass
+class DummyPanel: pass
+class DummyPropertyGroup: pass
+
+mock_types = MagicMock()
+mock_types.NodeSocket = DummyNodeSocket
+mock_types.Node = DummyNode
+mock_types.Object = DummyObject
+mock_types.Mesh = DummyMesh
+mock_types.NodeSocketShader = DummyNodeSocketShader
+mock_types.NodeCustomGroup = DummyNodeCustomGroup
+mock_types.NodeTree = DummyNodeTree
+mock_types.Operator = DummyOperator
+mock_types.Menu = DummyMenu
+mock_types.Panel = DummyPanel
+mock_types.PropertyGroup = DummyPropertyGroup
+
+mock_bpy = MagicMock()
+mock_bpy.context.scene.unit_settings.scale_length = 0.001
+mock_bpy.types = mock_types
+
+sys.modules['bpy'] = mock_bpy
+sys.modules['bpy.types'] = mock_types
+sys.modules['bpy.props'] = MagicMock()
+sys.modules['bpy.utils'] = MagicMock()
+sys.modules['bpy.utils.previews'] = MagicMock()
+sys.modules['bpy.app'] = MagicMock()
+sys.modules['bpy.app.handlers'] = MagicMock()
+
+from NodeSlicer.core.slice import SlicePaths, is_closed_path
+
+
+class TestSimplifyPaths(unittest.TestCase):
+
+    def test_colinear_points_elimination(self):
+        # A straight line segment with 50 redundant colinear points from (0,0,0) to (100,0,0)
+        xs = np.linspace(0, 100, 50)
+        ys = np.zeros_like(xs)
+        zs = np.zeros_like(xs)
+        line = np.column_stack([xs, ys, zs])
+
+        sp = SlicePaths([line])
+        res = sp.simplify(0.01)
+
+        self.assertEqual(len(res), 1)
+        self.assertEqual(len(res[0]), 2)
+        np.testing.assert_allclose(res[0][0], [0.0, 0.0, 0.0])
+        np.testing.assert_allclose(res[0][-1], [100.0, 0.0, 0.0])
+
+    def test_closed_square_with_subdivisions(self):
+        # 100 points on each edge of a 20x20 square (400 points total)
+        edge1 = np.column_stack([np.linspace(-10, 10, 100), np.full(100, -10.0), np.zeros(100)])
+        edge2 = np.column_stack([np.full(100, 10.0), np.linspace(-10, 10, 100), np.zeros(100)])
+        edge3 = np.column_stack([np.linspace(10, -10, 100), np.full(100, 10.0), np.zeros(100)])
+        edge4 = np.column_stack([np.full(100, -10.0), np.linspace(10, -10, 100), np.zeros(100)])
+
+        square_dense = np.concatenate([edge1[:-1], edge2[:-1], edge3[:-1], edge4], axis=0)
+        self.assertTrue(is_closed_path(square_dense))
+
+        sp = SlicePaths([square_dense])
+        res = sp.simplify(0.01, preserve_topology=True)
+
+        self.assertEqual(len(res), 1)
+        self.assertTrue(is_closed_path(res[0]))
+        # A clean 20x20 square should reduce to 5 vertices (4 corners + loop closure)
+        self.assertEqual(len(res[0]), 5)
+        np.testing.assert_allclose(res[0][0], res[0][-1])
+
+    def test_closed_loop_topology_preservation(self):
+        # Small circle with diameter 0.02mm, tolerance 0.05mm
+        # If preserve_topology=True, it should not collapse into < 4 vertices (stays a valid triangle loop)
+        theta = np.linspace(0, 2 * np.pi, 20)
+        r = 0.01
+        circle = np.column_stack([r * np.cos(theta), r * np.sin(theta), np.zeros(20)])
+        circle[-1] = circle[0]
+
+        sp = SlicePaths([circle])
+        res = sp.simplify(0.05, preserve_topology=True)
+
+        self.assertEqual(len(res), 1)
+        self.assertTrue(is_closed_path(res[0]))
+        self.assertGreaterEqual(len(res[0]), 4)
+
+    def test_attribute_transfer(self):
+        # Path with extrusion_rate and feedrate attributes
+        xs = np.linspace(0, 100, 50)
+        ys = np.zeros_like(xs)
+        zs = np.zeros_like(xs)
+        line = np.column_stack([xs, ys, zs])
+        ext = [np.linspace(1.0, 2.0, 50)]
+        feed = [np.linspace(1000.0, 2000.0, 50)]
+
+        sp = SlicePaths([line], extrusion_rate=ext, feedrate=feed)
+        res = sp.simplify(0.01)
+
+        self.assertEqual(len(res), 1)
+        self.assertEqual(len(res[0]), 2)
+        self.assertEqual(len(res.extrusion_rate[0]), 2)
+        self.assertEqual(len(res.feedrate[0]), 2)
+        np.testing.assert_allclose(res.extrusion_rate[0], [1.0, 2.0], atol=1e-3)
+        np.testing.assert_allclose(res.feedrate[0], [1000.0, 2000.0], atol=1e-3)
+
+    def test_batch_performance(self):
+        # 200 layers, 5 paths per layer = 1,000 paths (300,000 vertices)
+        paths = []
+        for z in np.linspace(0, 40, 200):
+            for r in [5.0, 10.0, 15.0, 20.0, 25.0]:
+                theta = np.linspace(0, 2 * np.pi, 300)
+                x = r * np.cos(theta) + np.random.normal(0, 0.005, 300)
+                y = r * np.sin(theta) + np.random.normal(0, 0.005, 300)
+                x[-1] = x[0]
+                y[-1] = y[0]
+                paths.append(np.column_stack([x, y, np.full(300, z)]))
+
+        sp = SlicePaths(paths)
+        import time
+        t0 = time.perf_counter()
+        res = sp.simplify(0.05)
+        duration = time.perf_counter() - t0
+
+        self.assertEqual(len(res), 1000)
+        self.assertLess(duration, 0.1, f"Batch simplify of 300,000 points took {duration:.4f}s (must be < 0.1s)")
+
+
+if __name__ == '__main__':
+    unittest.main()