diff --git a/core/c_seam.c b/core/c_seam.c index eaa765a..ce9054b 100644 --- a/core/c_seam.c +++ b/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] + ); + } +} + diff --git a/core/slice.py b/core/slice.py index c5dd693..fa88808 100644 --- a/core/slice.py +++ b/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 diff --git a/scratch/benchmark_simplify.py b/scratch/benchmark_simplify.py new file mode 100644 index 0000000..7811d7d --- /dev/null +++ b/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)}") diff --git a/scratch/test_batch_c_rdp.py b/scratch/test_batch_c_rdp.py new file mode 100644 index 0000000..11103bd --- /dev/null +++ b/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 +#include +#include + +#ifdef _OPENMP +#include +#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)}") diff --git a/scratch/test_c_rdp.py b/scratch/test_c_rdp.py new file mode 100644 index 0000000..5c0f2d7 --- /dev/null +++ b/scratch/test_c_rdp.py @@ -0,0 +1,227 @@ +import numpy as np +import time +import ctypes + +c_code = """ +#include +#include +#include + +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") diff --git a/tests/test_simplify.py b/tests/test_simplify.py new file mode 100644 index 0000000..86686ef --- /dev/null +++ b/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()