diff --git a/CMakeLists.txt b/CMakeLists.txt index c62b8b2..f83b0ae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,6 +31,7 @@ include_directories( # ----------------------------------------------------------------------------- add_library(cudaslicer_offset SHARED src/gpu_offset.cu + src/gpu_sdf.cu src/cpu_offset.cpp ) set_target_properties(cudaslicer_offset PROPERTIES @@ -45,6 +46,7 @@ target_link_libraries(cudaslicer_offset PUBLIC cudart) # ----------------------------------------------------------------------------- add_library(cudaslicer_offset_static STATIC src/gpu_offset.cu + src/gpu_sdf.cu src/cpu_offset.cpp ) set_target_properties(cudaslicer_offset_static PROPERTIES diff --git a/examples/gpu_sdf_csg_example.py b/examples/gpu_sdf_csg_example.py new file mode 100644 index 0000000..40cb4ed --- /dev/null +++ b/examples/gpu_sdf_csg_example.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +""" +GPU-Resident Signed Distance Field (SDF) CSG & Infill Example +============================================================= +Demonstrates chaining 3D mesh slicing, Boolean CSG operations (Union, Intersection, +Difference, Smooth Blend), constant-distance offsetting, volume shelling, procedural +3D gyroid infill generation, and contour vectorization entirely in GPU VRAM without +host-to-device memory roundtrips. +""" + +import os +import sys +import time +import numpy as np + +# Ensure local build is in path +sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "python")) +import cudaslicer +import trimesh + + +def main(): + print("=" * 70) + print(" CUDASlicer: GPU Signed Distance Field (SDF) & CSG Engine") + print("=" * 70) + + if not cudaslicer.has_cuda(): + print("[Error] CUDA-capable GPU is required for this demo.") + return + + # 1. Load or generate 3D models + stl_path = os.path.join(os.path.dirname(__file__), "..", "benchy.stl") + if os.path.exists(stl_path): + print(f"\n[1] Loading main model: {stl_path}") + mesh_a = trimesh.load(stl_path) + else: + print("\n[1] Generating synthetic main cylinder...") + mesh_a = trimesh.creation.cylinder(radius=15.0, height=30.0) + + # Generate a modifier box volume to subtract from the model + print("[1] Generating modifier cutting box...") + mesh_b = trimesh.creation.box(extents=(12.0, 60.0, 25.0)) + mesh_b.apply_translation((0.0, 0.0, 15.0)) + + # 2. Direct GPU Slicing to SDF Stacks + layer_height_mm = 0.20 + pitch_mm = 0.10 + + t0 = time.perf_counter() + sdf_a = cudaslicer.slice_mesh_to_sdf(mesh_a, layer_height=layer_height_mm, pitch_mm=pitch_mm) + sdf_b = cudaslicer.slice_mesh_to_sdf(mesh_b, layer_height=layer_height_mm, pitch_mm=pitch_mm) + t_slice = (time.perf_counter() - t0) * 1000.0 + + print(f"\n[2] Direct GPU Slicing & Signed Distance Generation ({t_slice:.2f} ms):") + print(f" - Main Model Stack: {sdf_a}") + print(f" Layer Z range (µm): {sdf_a.layer_z_um[0]} µm to {sdf_a.layer_z_um[-1]} µm ({len(sdf_a.layer_z_um)} layers)") + print(f" - Modifier Box Stack: {sdf_b}") + print(f" Layer Z range (µm): {sdf_b.layer_z_um[0]} µm to {sdf_b.layer_z_um[-1]} µm ({len(sdf_b.layer_z_um)} layers)") + + # 3. GPU-Resident Boolean Operations (Zero CPU transfers) + t0 = time.perf_counter() + # Chained expression: (Model A \ (Modifier B + 0.5mm offset)) + modifier_expanded = sdf_b.offset(delta_mm=0.5) + sdf_carved = sdf_a.difference(modifier_expanded) + t_csg = (time.perf_counter() - t0) * 1000.0 + + print(f"\n[3] Chained GPU CSG (Carve / Subtraction) ({t_csg:.2f} ms):") + print(f" - Result Stack: {sdf_carved}") + print(f" - Automatic spatial bounding grid expansion: {sdf_carved.shape}") + + # 4. Procedural 3D Gyroid Infill Generation on GPU + t0 = time.perf_counter() + sdf_gyroid = sdf_carved.generate_infill( + pattern="gyroid", + cell_size_mm=4.0, + density=0.25, + wall_offset_mm=1.20 + ) + t_infill = (time.perf_counter() - t0) * 1000.0 + + print(f"\n[4] GPU Procedural Gyroid Infill Generation ({t_infill:.2f} ms):") + print(f" - Infill Stack: {sdf_gyroid}") + + # 5. Vectorize directly to polygon perimeters and infill contours + t0 = time.perf_counter() + layer_results = sdf_gyroid.vectorize( + perimeters_count=2, + layer_width=0.45, + export_infill_mask=True + ) + t_vec = (time.perf_counter() - t0) * 1000.0 + + total_contours = sum( + sum(len(p["contours"]) for p in l.perimeters) + len(l.infill_contours) + for l in layer_results + ) + + print(f"\n[5] Vectorization to Slicer Contours ({t_vec:.2f} ms):") + print(f" - Total Layers Output: {len(layer_results)}") + print(f" - Total Vector Contours Extracted: {total_contours}") + + # Sample mid layer + mid_idx = len(layer_results) // 2 + sample_layer = layer_results[mid_idx] + print(f"\n[6] Sample Mid Layer (Layer {mid_idx}):") + print(f" - Z Height: {sample_layer.z:.3f} mm ({sample_layer.z_um} µm)") + print(f" - Perimeters: {len(sample_layer.perimeters)}") + for i, p in enumerate(sample_layer.perimeters): + print(f" * Perimeter {i+1} (offset {p['offset_mm']:.2f} mm): {len(p['contours'])} contour loop(s)") + print(f" - Infill Contours: {len(sample_layer.infill_contours)} contour loop(s)") + + print("\n" + "=" * 70) + print(" Pipeline completed successfully with 100% GPU-resident computation!") + print("=" * 70) + + +if __name__ == "__main__": + main() diff --git a/include/gpu_offset.h b/include/gpu_offset.h index 1bf67c6..4bab76d 100644 --- a/include/gpu_offset.h +++ b/include/gpu_offset.h @@ -105,6 +105,8 @@ typedef struct { int use_gpu_rasterizer; // 1 = GPU-accelerated rasterization (when CUDA enabled), 0 = CPU multi-threaded float first_layer_height; // First layer height in mm (default: 0.0f -> uses layer_height) int export_infill_mask; // 1 = export binary 2D infill mask into GpuLayerResult (default: 1) + float simplify_tolerance_mm; // RDP simplification tolerance in mm (default: 0.02 mm, 0 to disable) + int subpixel_lerp; // 1 = Continuous SDF sub-pixel LERP vectorization (default: 1), 0 = discrete bitpack int backend; // GpuOffsetBackend: 0 = AUTO, 1 = CUDA, 2 = CPU (default: 0) } GpuOffsetConfig; diff --git a/include/gpu_sdf.h b/include/gpu_sdf.h new file mode 100644 index 0000000..a6b6717 --- /dev/null +++ b/include/gpu_sdf.h @@ -0,0 +1,209 @@ +#ifndef GPU_SDF_H +#define GPU_SDF_H + +#ifdef _WIN32 + #ifdef GPU_OFFSET_EXPORTS + #define GPU_SDF_API __declspec(dllexport) + #else + #define GPU_SDF_API __declspec(dllimport) + #endif +#else + #define GPU_SDF_API __attribute__((visibility("default"))) +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +#include +#include +#include "gpu_offset.h" + +/** + * @brief CSG Operation Types for GPU SDF combinations + */ +typedef enum { + GPU_SDF_OP_UNION = 0, // A U B: max(d_A, d_B) + GPU_SDF_OP_INTERSECTION = 1, // A ^ B: min(d_A, d_B) + GPU_SDF_OP_DIFFERENCE = 2, // A \ B: min(d_A, -d_B) + GPU_SDF_OP_SMOOTH_UNION = 3, // Smooth union with blend radius k + GPU_SDF_OP_SMOOTH_INTERSECT= 4, // Smooth intersection with blend radius k + GPU_SDF_OP_SMOOTH_DIFF = 5 // Smooth difference with blend radius k +} GpuSdfOp; + +/** + * @brief Layer matching policy when operating on two stacks + */ +typedef enum { + GPU_SDF_LAYER_MATCH_A_DOMINANT = 0, // Keep layers matching stack A (resample B at A's Z) + GPU_SDF_LAYER_MATCH_UNION = 1, // Output union of all unique layer Z heights + GPU_SDF_LAYER_MATCH_INTERSECTION = 2 // Output only exact matching layer Z heights +} GpuSdfLayerMatchMode; + +/** + * @brief Procedural 3D Infill Pattern Types + */ +typedef enum { + GPU_SDF_INFILL_SOLID = 0, + GPU_SDF_INFILL_GRID = 1, + GPU_SDF_INFILL_GYROID = 2, + GPU_SDF_INFILL_SCHWARZ_P = 3, + GPU_SDF_INFILL_SCHWARZ_D = 4 +} GpuSdfInfillPattern; + +/** + * @brief Spatial metadata for a 2D/3D GPU Signed Distance Field stack + * All spatial coordinates and pitches are stored in integer micrometers (µm) + */ +typedef struct { + int32_t min_x_um; // Minimum X bound in micrometers (1 µm = 0.001 mm) + int32_t min_y_um; // Minimum Y bound in micrometers + int32_t pitch_um; // Pixel raster pitch in micrometers (e.g. 100 µm = 0.1 mm) + int width; // Grid columns (X) + int height; // Grid rows (Y) + size_t num_layers; // Total number of layers (Z) +} GpuSdfMetadata; + +/** + * @brief Opaque handle representing a GPU-resident Signed Distance Field stack. + * Distance values are in millimeters (float): + * d > 0 : Inside the object (positive Euclidean distance to boundary) + * d < 0 : Outside the object (negative Euclidean distance to boundary) + * d = 0 : Exact object boundary + */ +typedef struct GpuSdfStack_T* GpuSdfStackHandle; + +/** + * @brief Create an empty GPU SDF stack with given metadata and layer Z heights (in µm). + */ +GPU_SDF_API GpuSdfStackHandle gpu_sdf_create( + const GpuSdfMetadata* meta, + const int32_t* layer_z_heights_um +); + +/** + * @brief Create a GPU SDF stack from 3D mesh triangles directly on GPU. + */ +GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_mesh( + const GpuTriangle3D* triangles, + size_t num_triangles, + float min_x, float min_y, float min_z, + float max_x, float max_y, float max_z, + float layer_height_mm, + float first_layer_height_mm, + float pitch_mm, + float pad_mm +); + +/** + * @brief Create a GPU SDF stack from raw host binary occupancy masks (0 = background, >0 = object). + */ +GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_masks( + const uint8_t* h_masks, + int width, + int height, + size_t num_layers, + const int32_t* layer_z_heights_um, + int32_t min_x_um, + int32_t min_y_um, + int32_t pitch_um +); + +/** + * @brief Clone an existing GPU SDF stack (allocates new GPU memory). + */ +GPU_SDF_API GpuSdfStackHandle gpu_sdf_clone(GpuSdfStackHandle stack); + +/** + * @brief Free GPU memory allocated for a GPU SDF stack. + */ +GPU_SDF_API void gpu_sdf_destroy(GpuSdfStackHandle stack); + +/** + * @brief Query metadata of a GPU SDF stack. + */ +GPU_SDF_API int gpu_sdf_get_metadata(GpuSdfStackHandle stack, GpuSdfMetadata* out_meta); + +/** + * @brief Retrieve layer Z heights in micrometers (copies to host buffer of size num_layers). + */ +GPU_SDF_API int gpu_sdf_get_layer_z_um(GpuSdfStackHandle stack, int32_t* out_z_um); + +/** + * @brief Retrieve raw SDF distance data (float mm) copied from GPU to Host. + * Buffer size must be: width * height * num_layers * sizeof(float). + */ +GPU_SDF_API int gpu_sdf_to_host(GpuSdfStackHandle stack, float* h_out_data); + +/** + * @brief Get raw GPU device pointer for zero-copy interop (CuPy / PyTorch). + */ +GPU_SDF_API float* gpu_sdf_get_device_ptr(GpuSdfStackHandle stack); + +/** + * @brief Get raw GPU device pointer for layer Z heights (int32_t* on device). + */ +GPU_SDF_API int32_t* gpu_sdf_get_device_z_ptr(GpuSdfStackHandle stack); + +/** + * @brief Execute a CSG operation between two GPU SDF stacks without host roundtrips. + * Handles spatial bounding box differences and layer Z matching automatically. + * + * @param a First operand stack (GPU) + * @param b Second operand stack (GPU) + * @param op CSG operation (UNION, INTERSECTION, DIFFERENCE, SMOOTH_...) + * @param match_mode Layer Z matching mode + * @param blend_radius_mm Blend radius parameter for smooth operations + * @return GpuSdfStackHandle Resulting stack on GPU (caller must free with gpu_sdf_destroy) + */ +GPU_SDF_API GpuSdfStackHandle gpu_sdf_csg_op( + GpuSdfStackHandle a, + GpuSdfStackHandle b, + GpuSdfOp op, + GpuSdfLayerMatchMode match_mode, + float blend_radius_mm +); + +/** + * @brief Apply a constant distance offset (dilation/expansion if delta > 0, inset/erosion if delta < 0). + */ +GPU_SDF_API GpuSdfStackHandle gpu_sdf_offset( + GpuSdfStackHandle stack, + float delta_mm +); + +/** + * @brief Apply hollowing / shelling with a given wall thickness. + */ +GPU_SDF_API GpuSdfStackHandle gpu_sdf_shell( + GpuSdfStackHandle stack, + float thickness_mm +); + +/** + * @brief Generate a procedural 3D infill pattern intersected with the interior region. + */ +GPU_SDF_API GpuSdfStackHandle gpu_sdf_generate_infill( + GpuSdfStackHandle stack, + GpuSdfInfillPattern pattern, + float cell_size_mm, + float density_ratio, + float wall_offset_mm +); + +/** + * @brief Vectorize GPU SDF stack directly to polygon perimeter contours and infill mask. + */ +GPU_SDF_API int gpu_sdf_vectorize( + GpuSdfStackHandle stack, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + size_t* out_num_layers, + GpuOffsetStats* stats +); + +#ifdef __cplusplus +} +#endif + +#endif // GPU_SDF_H diff --git a/python/cudaslicer/__init__.py b/python/cudaslicer/__init__.py index cd18bda..2e66c68 100644 --- a/python/cudaslicer/__init__.py +++ b/python/cudaslicer/__init__.py @@ -347,11 +347,225 @@ def slice_mesh_and_offset( # Alias for convenience slice_mesh = slice_mesh_and_offset + +class SdfStack: + """ + GPU-resident Signed Distance Field (SDF) stack representing a 3D model or modifier volume. + Enables chained GPU-accelerated Boolean CSG operations (union, intersection, difference, + smooth blend, offset, shelling, procedural 3D infill) without host-to-device data transfers. + """ + def __init__(self, backend_stack: Any): + self._stack = backend_stack + + @property + def metadata(self) -> Dict[str, Any]: + """Dictionary containing spatial bounds, grid dimensions, and layer count.""" + return self._stack.get_metadata() + + @property + def layer_z_um(self) -> np.ndarray: + """1D array of layer heights in micrometers (µm).""" + return self._stack.get_layer_z_um() + + @property + def layer_z_mm(self) -> np.ndarray: + """1D array of layer heights in millimeters (mm).""" + return self._stack.get_layer_z_mm() + + @property + def shape(self) -> Tuple[int, int, int]: + """Grid shape (num_layers, height, width).""" + m = self.metadata + return (m["num_layers"], m["height"], m["width"]) + + def numpy(self) -> np.ndarray: + """Download raw 3D distance field array (num_layers, H, W) to CPU NumPy float32.""" + return self._stack.to_numpy() + + def clone(self) -> "SdfStack": + """Duplicate this SDF stack in GPU VRAM.""" + return SdfStack(self._stack.clone()) + + def union(self, other: "SdfStack", blend_radius_mm: float = 0.0) -> "SdfStack": + """Compute Boolean Union (A ∪ B) on GPU. Optional blend_radius_mm for smooth fillet union.""" + if not isinstance(other, SdfStack): + raise TypeError(f"Operand must be an SdfStack, got {type(other)}") + return SdfStack(self._stack.op_union(other._stack, float(blend_radius_mm))) + + def intersection(self, other: "SdfStack", blend_radius_mm: float = 0.0) -> "SdfStack": + """Compute Boolean Intersection (A ∩ B) on GPU. Optional blend_radius_mm for smooth intersection.""" + if not isinstance(other, SdfStack): + raise TypeError(f"Operand must be an SdfStack, got {type(other)}") + return SdfStack(self._stack.op_intersection(other._stack, float(blend_radius_mm))) + + def difference(self, other: "SdfStack", blend_radius_mm: float = 0.0) -> "SdfStack": + """Compute Boolean Difference / Subtraction (A \\ B) on GPU. Optional blend_radius_mm for smooth difference.""" + if not isinstance(other, SdfStack): + raise TypeError(f"Operand must be an SdfStack, got {type(other)}") + return SdfStack(self._stack.op_difference(other._stack, float(blend_radius_mm))) + + def offset(self, delta_mm: float) -> "SdfStack": + """Constant-distance offset/dilation (delta > 0) or inset/erosion (delta < 0) on GPU.""" + return SdfStack(self._stack.offset(float(delta_mm))) + + def shell(self, thickness_mm: float) -> "SdfStack": + """Hollow / shell the volume with specified wall thickness on GPU.""" + return SdfStack(self._stack.shell(float(thickness_mm))) + + def generate_infill( + self, + pattern: str = "gyroid", + cell_size_mm: float = 5.0, + density: float = 0.20, + wall_offset_mm: float = 1.0 + ) -> "SdfStack": + """Generate continuous 3D procedural lattice (gyroid, schwarz_p, schwarz_d, grid) on GPU.""" + return SdfStack(self._stack.generate_infill(str(pattern), float(cell_size_mm), float(density), float(wall_offset_mm))) + + def vectorize( + self, + perimeters_count: int = 2, + layer_width: float = 0.45, + export_infill_mask: bool = True, + simplify_tolerance_mm: float = 0.02, + subpixel_lerp: bool = True, + offset_distances: Optional[Sequence[float]] = None + ) -> List[LayerResult]: + """ + Vectorize GPU SDF stack directly to smooth polygon perimeter contours and infill masks. + + Parameters + ---------- + perimeters_count : int, default=2 + Number of perimeter passes. + layer_width : float, default=0.45 + Extrusion line width in mm. + export_infill_mask : bool, default=True + Whether to export binary 2D infill mask. + simplify_tolerance_mm : float, default=0.02 + Ramer-Douglas-Peucker (RDP) polygon simplification tolerance in mm (eliminates micro-steps). + subpixel_lerp : bool, default=True + Continuous SDF linear interpolation for sub-pixel boundary placement at arbitrary angles. + offset_distances : sequence of float, optional + Custom offset distances. + """ + offsets_float = [float(d) for d in offset_distances] if offset_distances is not None else [] + raw_res = self._stack.vectorize( + int(perimeters_count), + float(layer_width), + bool(export_infill_mask), + float(simplify_tolerance_mm), + bool(subpixel_lerp), + offsets_float + ) + return [LayerResult(r) for r in raw_res] + + # Operator Overloads for fluent Python syntax + def __or__(self, other: "SdfStack") -> "SdfStack": + return self.union(other) + + def __and__(self, other: "SdfStack") -> "SdfStack": + return self.intersection(other) + + def __sub__(self, other: "SdfStack") -> "SdfStack": + return self.difference(other) + + def __add__(self, delta: float) -> "SdfStack": + if isinstance(delta, (int, float)): + return self.offset(float(delta)) + raise TypeError(f"Unsupported operand type for +: 'SdfStack' and '{type(delta)}'") + + def __repr__(self) -> str: + m = self.metadata + return f"" + + +def slice_mesh_to_sdf( + mesh_or_triangles: Any, + layer_height: float = 0.20, + first_layer_height: Optional[float] = None, + pitch_mm: float = 0.10, + pad_mm: float = 1.5 +) -> SdfStack: + """ + Slice a 3D mesh directly to a GPU-resident Signed Distance Field (SDF) stack. + + Parameters + ---------- + mesh_or_triangles : trimesh.Trimesh OR np.ndarray + Mesh or triangle array of shape (N, 3, 3) or (N, 9). + layer_height : float, default=0.20 + Z slicing layer height in millimeters. + first_layer_height : float, optional + First layer height in millimeters. Defaults to `layer_height` if None. + pitch_mm : float, default=0.10 + Pixel raster pitch in millimeters. + pad_mm : float, default=1.5 + Bounding box padding in millimeters. + + Returns + ------- + SdfStack + GPU-resident SDF stack object ready for CSG operations or vectorization. + """ + if not _HAS_BACKEND or _backend is None: + raise RuntimeError("CUDASlicer backend extension (_cudaslicer_backend) is not available.") + + if hasattr(mesh_or_triangles, "triangles"): + tris = np.ascontiguousarray(mesh_or_triangles.triangles, dtype=np.float32) + bounds = mesh_or_triangles.bounds + min_x, min_y, min_z = float(bounds[0][0]), float(bounds[0][1]), float(bounds[0][2]) + max_x, max_y, max_z = float(bounds[1][0]), float(bounds[1][1]), float(bounds[1][2]) + else: + tris = np.ascontiguousarray(mesh_or_triangles, dtype=np.float32) + if tris.ndim == 2 and tris.shape[1] == 9: + tris_3d = tris.reshape(-1, 3, 3) + else: + tris_3d = tris + min_x = float(tris_3d[:, :, 0].min()) + max_x = float(tris_3d[:, :, 0].max()) + min_y = float(tris_3d[:, :, 1].min()) + max_y = float(tris_3d[:, :, 1].max()) + min_z = float(tris_3d[:, :, 2].min()) + max_z = float(tris_3d[:, :, 2].max()) + + h_first = float(first_layer_height) if first_layer_height is not None and first_layer_height > 0 else float(layer_height) + + backend_sdf = _backend.slice_mesh_to_sdf( + tris, min_x, min_y, min_z, max_x, max_y, max_z, + float(layer_height), float(h_first), float(pitch_mm), float(pad_mm) + ) + return SdfStack(backend_sdf) + + +def masks_to_sdf( + masks: np.ndarray, + layer_z_um: Optional[Sequence[int]] = None, + min_x_mm: float = 0.0, + min_y_mm: float = 0.0, + pitch_mm: float = 0.10 +) -> SdfStack: + """ + Compute a GPU-resident Signed Distance Field from 2D (H, W) or 3D (N, H, W) binary masks. + """ + if not _HAS_BACKEND or _backend is None: + raise RuntimeError("CUDASlicer backend extension (_cudaslicer_backend) is not available.") + + masks_u8 = np.ascontiguousarray(masks, dtype=np.uint8) + z_arr = np.ascontiguousarray(layer_z_um, dtype=np.int32) if layer_z_um is not None else np.array([], dtype=np.int32) + + backend_sdf = _backend.masks_to_sdf(masks_u8, z_arr, float(min_x_mm), float(min_y_mm), float(pitch_mm)) + return SdfStack(backend_sdf) + + __all__ = [ "offset_masks", "offset_polygons", "slice_mesh_and_offset", "slice_mesh", + "slice_mesh_to_sdf", + "masks_to_sdf", + "SdfStack", "LayerResult", "is_available", "has_cuda", @@ -361,3 +575,4 @@ __all__ = [ "get_cuda_device_count", "__version__", ] + diff --git a/python/cudaslicer/bindings.cpp b/python/cudaslicer/bindings.cpp index 9c9e60b..f315f1e 100644 --- a/python/cudaslicer/bindings.cpp +++ b/python/cudaslicer/bindings.cpp @@ -362,8 +362,352 @@ py::list py_slice_mesh_and_offset( return py_layer_results; } +#include "gpu_sdf.h" +#include + +class PyGpuSdfStack { +private: + GpuSdfStackHandle handle_; + +public: + PyGpuSdfStack(GpuSdfStackHandle handle) : handle_(handle) { + if (!handle_) throw std::runtime_error("Failed to allocate GPU SDF Stack."); + } + + ~PyGpuSdfStack() { + if (handle_) { + gpu_sdf_destroy(handle_); + handle_ = nullptr; + } + } + + PyGpuSdfStack(const PyGpuSdfStack&) = delete; + PyGpuSdfStack& operator=(const PyGpuSdfStack&) = delete; + + PyGpuSdfStack(PyGpuSdfStack&& other) noexcept : handle_(other.handle_) { + other.handle_ = nullptr; + } + + PyGpuSdfStack& operator=(PyGpuSdfStack&& other) noexcept { + if (this != &other) { + if (handle_) gpu_sdf_destroy(handle_); + handle_ = other.handle_; + other.handle_ = nullptr; + } + return *this; + } + + GpuSdfStackHandle get_handle() const { return handle_; } + + std::shared_ptr clone() const { + return std::make_shared(gpu_sdf_clone(handle_)); + } + + py::dict get_metadata() const { + GpuSdfMetadata meta; + if (gpu_sdf_get_metadata(handle_, &meta) != GPU_OFFSET_SUCCESS) { + throw std::runtime_error("Failed to get SDF metadata."); + } + py::dict d; + d["width"] = meta.width; + d["height"] = meta.height; + d["num_layers"] = meta.num_layers; + d["min_x_um"] = meta.min_x_um; + d["min_y_um"] = meta.min_y_um; + d["pitch_um"] = meta.pitch_um; + d["min_x_mm"] = (float)meta.min_x_um * 0.001f; + d["min_y_mm"] = (float)meta.min_y_um * 0.001f; + d["pitch_mm"] = (float)meta.pitch_um * 0.001f; + return d; + } + + py::array_t get_layer_z_um() const { + GpuSdfMetadata meta; + gpu_sdf_get_metadata(handle_, &meta); + py::array_t arr((ssize_t)meta.num_layers); + auto buf = arr.request(); + gpu_sdf_get_layer_z_um(handle_, static_cast(buf.ptr)); + return arr; + } + + py::array_t get_layer_z_mm() const { + GpuSdfMetadata meta; + gpu_sdf_get_metadata(handle_, &meta); + std::vector z_um(meta.num_layers); + gpu_sdf_get_layer_z_um(handle_, z_um.data()); + py::array_t arr((ssize_t)meta.num_layers); + auto r = arr.mutable_unchecked<1>(); + for (size_t i = 0; i < meta.num_layers; ++i) { + r(i) = (float)z_um[i] * 0.001f; + } + return arr; + } + + py::array_t to_numpy() const { + GpuSdfMetadata meta; + gpu_sdf_get_metadata(handle_, &meta); + std::vector shape = { (ssize_t)meta.num_layers, (ssize_t)meta.height, (ssize_t)meta.width }; + py::array_t out_arr(shape); + auto buf = out_arr.request(); + if (gpu_sdf_to_host(handle_, static_cast(buf.ptr)) != GPU_OFFSET_SUCCESS) { + throw std::runtime_error("Failed to copy SDF to NumPy host buffer."); + } + return out_arr; + } + + std::shared_ptr csg_op( + const PyGpuSdfStack& other, + const std::string& op_name, + const std::string& match_mode_name = "dominant", + float blend_radius_mm = 0.0f + ) const { + GpuSdfOp op = GPU_SDF_OP_UNION; + if (op_name == "union" || op_name == "or" || op_name == "+") op = GPU_SDF_OP_UNION; + else if (op_name == "intersection" || op_name == "and" || op_name == "&") op = GPU_SDF_OP_INTERSECTION; + else if (op_name == "difference" || op_name == "subtract" || op_name == "-") op = GPU_SDF_OP_DIFFERENCE; + else if (op_name == "smooth_union") op = GPU_SDF_OP_SMOOTH_UNION; + else if (op_name == "smooth_intersection" || op_name == "smooth_intersect") op = GPU_SDF_OP_SMOOTH_INTERSECT; + else if (op_name == "smooth_difference" || op_name == "smooth_diff") op = GPU_SDF_OP_SMOOTH_DIFF; + else throw std::invalid_argument("Unknown CSG op: " + op_name); + + GpuSdfLayerMatchMode match = GPU_SDF_LAYER_MATCH_A_DOMINANT; + if (match_mode_name == "dominant" || match_mode_name == "a") match = GPU_SDF_LAYER_MATCH_A_DOMINANT; + else if (match_mode_name == "union") match = GPU_SDF_LAYER_MATCH_UNION; + else if (match_mode_name == "intersection" || match_mode_name == "intersect") match = GPU_SDF_LAYER_MATCH_INTERSECTION; + else throw std::invalid_argument("Unknown layer match mode: " + match_mode_name); + + GpuSdfStackHandle res = gpu_sdf_csg_op(handle_, other.handle_, op, match, blend_radius_mm); + if (!res) throw std::runtime_error("CSG operation failed."); + return std::make_shared(res); + } + + std::shared_ptr op_union(const PyGpuSdfStack& other, float blend_k = 0.0f) const { + return csg_op(other, blend_k > 0.0f ? "smooth_union" : "union", "dominant", blend_k); + } + + std::shared_ptr op_intersection(const PyGpuSdfStack& other, float blend_k = 0.0f) const { + return csg_op(other, blend_k > 0.0f ? "smooth_intersection" : "intersection", "dominant", blend_k); + } + + std::shared_ptr op_difference(const PyGpuSdfStack& other, float blend_k = 0.0f) const { + return csg_op(other, blend_k > 0.0f ? "smooth_difference" : "difference", "dominant", blend_k); + } + + std::shared_ptr offset(float delta_mm) const { + GpuSdfStackHandle res = gpu_sdf_offset(handle_, delta_mm); + if (!res) throw std::runtime_error("SDF offset failed."); + return std::make_shared(res); + } + + std::shared_ptr shell(float thickness_mm) const { + GpuSdfStackHandle res = gpu_sdf_shell(handle_, thickness_mm); + if (!res) throw std::runtime_error("SDF shell failed."); + return std::make_shared(res); + } + + std::shared_ptr generate_infill( + const std::string& pattern_str = "gyroid", + float cell_size_mm = 5.0f, + float density_ratio = 0.20f, + float wall_offset_mm = 1.0f + ) const { + GpuSdfInfillPattern pat = GPU_SDF_INFILL_GYROID; + if (pattern_str == "gyroid") pat = GPU_SDF_INFILL_GYROID; + else if (pattern_str == "schwarz_p" || pattern_str == "schwarz-p") pat = GPU_SDF_INFILL_SCHWARZ_P; + else if (pattern_str == "schwarz_d" || pattern_str == "schwarz-d") pat = GPU_SDF_INFILL_SCHWARZ_D; + else if (pattern_str == "grid") pat = GPU_SDF_INFILL_GRID; + else throw std::invalid_argument("Unknown infill pattern: " + pattern_str); + + GpuSdfStackHandle res = gpu_sdf_generate_infill(handle_, pat, cell_size_mm, density_ratio, wall_offset_mm); + if (!res) throw std::runtime_error("SDF infill generation failed."); + return std::make_shared(res); + } + + py::list vectorize( + int perimeters_count = 2, + float layer_width = 0.45f, + bool export_infill_mask = true, + float simplify_tolerance_mm = 0.02f, + bool subpixel_lerp = true, + const std::vector& offset_distances = std::vector() + ) const { + GpuOffsetConfig config = {0}; + config.perimeters_count = perimeters_count; + config.layer_width = layer_width; + config.export_infill_mask = export_infill_mask ? 1 : 0; + config.simplify_tolerance_mm = simplify_tolerance_mm; + config.subpixel_lerp = subpixel_lerp ? 1 : 0; + config.use_gpu_rasterizer = 1; + config.backend = GPU_OFFSET_BACKEND_CUDA; + + if (!offset_distances.empty()) { + config.offset_distances_mm = offset_distances.data(); + config.num_offsets = offset_distances.size(); + } + + GpuLayerResult* results = nullptr; + size_t out_num_layers = 0; + GpuOffsetStats stats = {0}; + + int ret = gpu_sdf_vectorize(handle_, &config, &results, &out_num_layers, &stats); + if (ret != GPU_OFFSET_SUCCESS || !results) { + throw std::runtime_error("SDF vectorization failed."); + } + + py::list py_layer_results; + for (size_t l = 0; l < out_num_layers; ++l) { + py::dict py_layer_dict; + py_layer_dict["z"] = results[l].z_height; + py_layer_dict["z_um"] = (int32_t)std::round(results[l].z_height * 1000.0f); + + // Perimeters + py::list py_perimeters; + for (size_t o = 0; o < results[l].num_perimeters; ++o) { + py::dict py_perim_dict; + py_perim_dict["offset_mm"] = results[l].perimeters[o].offset_distance_mm; + + py::list py_cnt_list; + for (size_t c = 0; c < results[l].perimeters[o].num_contours; ++c) { + const auto& cnt = results[l].perimeters[o].contours[c]; + py::array_t pts_arr({ (ssize_t)cnt.num_points, (ssize_t)2 }); + auto r = pts_arr.mutable_unchecked<2>(); + for (size_t p = 0; p < cnt.num_points; ++p) { + r(p, 0) = cnt.points[p].x; + r(p, 1) = cnt.points[p].y; + } + py_cnt_list.append(pts_arr); + } + py_perim_dict["contours"] = py_cnt_list; + py_perimeters.append(py_perim_dict); + } + py_layer_dict["perimeters"] = py_perimeters; + + // Infill Contours + py::list py_infill_cnts; + for (size_t c = 0; c < results[l].infill.num_contours; ++c) { + const auto& cnt = results[l].infill.contours[c]; + py::array_t pts_arr({ (ssize_t)cnt.num_points, (ssize_t)2 }); + auto r = pts_arr.mutable_unchecked<2>(); + for (size_t p = 0; p < cnt.num_points; ++p) { + r(p, 0) = cnt.points[p].x; + r(p, 1) = cnt.points[p].y; + } + py_infill_cnts.append(pts_arr); + } + py_layer_dict["infill_contours"] = py_infill_cnts; + py_layer_dict["infill_offset_mm"] = results[l].infill.offset_distance_mm; + + // Infill 2D Mask + if (results[l].infill.mask != nullptr) { + int H = results[l].infill.height; + int W = results[l].infill.width; + py::array_t mask_arr({ (ssize_t)H, (ssize_t)W }); + py::buffer_info mb = mask_arr.request(); + memcpy(mb.ptr, results[l].infill.mask, (size_t)H * W * sizeof(uint8_t)); + py_layer_dict["infill_mask"] = mask_arr; + } + + py_layer_results.append(py_layer_dict); + } + + gpu_offset_free_results(results, out_num_layers); + return py_layer_results; + } + + py::dict cuda_array_interface() const { + GpuSdfMetadata meta; + gpu_sdf_get_metadata(handle_, &meta); + float* d_ptr = gpu_sdf_get_device_ptr(handle_); + + py::dict iface; + iface["shape"] = py::make_tuple(meta.num_layers, meta.height, meta.width); + iface["typestr"] = " py_slice_mesh_to_sdf( + py::array_t triangles_arr, + float min_x, float min_y, float min_z, + float max_x, float max_y, float max_z, + float layer_height_mm = 0.20f, + float first_layer_height_mm = 0.0f, + float pitch_mm = 0.10f, + float pad_mm = 1.5f +) { + py::buffer_info buf = triangles_arr.request(); + size_t num_tris = 0; + if (buf.ndim == 3 && buf.shape[1] == 3 && buf.shape[2] == 3) { + num_tris = buf.shape[0]; + } else if (buf.ndim == 2 && buf.shape[1] == 9) { + num_tris = buf.shape[0]; + } else { + throw std::runtime_error("Triangles array must be of shape (N, 3, 3) or (N, 9) float32."); + } + + if (num_tris == 0) throw std::runtime_error("Triangle array is empty."); + + const GpuTriangle3D* tris_ptr = static_cast(buf.ptr); + + GpuSdfStackHandle handle = gpu_sdf_from_mesh( + tris_ptr, num_tris, + min_x, min_y, min_z, max_x, max_y, max_z, + layer_height_mm, first_layer_height_mm, pitch_mm, pad_mm + ); + + if (!handle) throw std::runtime_error("Failed to slice mesh to GPU SDF stack."); + return std::make_shared(handle); +} + +std::shared_ptr py_masks_to_sdf( + py::array_t in_masks, + py::array_t layer_z_um_arr, + float min_x_mm = 0.0f, + float min_y_mm = 0.0f, + float pitch_mm = 0.10f +) { + py::buffer_info buf = in_masks.request(); + size_t num_layers = 1; + int height = 0; + int width = 0; + + if (buf.ndim == 2) { + num_layers = 1; + height = (int)buf.shape[0]; + width = (int)buf.shape[1]; + } else if (buf.ndim == 3) { + num_layers = (size_t)buf.shape[0]; + height = (int)buf.shape[1]; + width = (int)buf.shape[2]; + } else { + throw std::runtime_error("Input mask must be 2D (H, W) or 3D (N, H, W)."); + } + + py::buffer_info z_buf = layer_z_um_arr.request(); + const int32_t* z_ptr = nullptr; + if (z_buf.size == (ssize_t)num_layers) { + z_ptr = static_cast(z_buf.ptr); + } + + int32_t min_x_um = (int32_t)std::round(min_x_mm * 1000.0f); + int32_t min_y_um = (int32_t)std::round(min_y_mm * 1000.0f); + int32_t pitch_um = (int32_t)std::round(pitch_mm * 1000.0f); + + GpuSdfStackHandle handle = gpu_sdf_from_masks( + static_cast(buf.ptr), + width, height, num_layers, + z_ptr, min_x_um, min_y_um, pitch_um + ); + + if (!handle) throw std::runtime_error("Failed to create GPU SDF stack from masks."); + return std::make_shared(handle); +} + PYBIND11_MODULE(_cudaslicer_backend, m) { - m.doc() = "GPU/CPU Slicing & SDF/JFA Offsetting C++/CUDA Python Extension"; + m.doc() = "GPU/CPU Slicing, Signed Distance Fields (SDF), CSG Booleans & Vectorizer Extension"; m.def("has_cuda", &gpu_offset_has_cuda, "Returns 1 if a functional CUDA GPU is present, 0 otherwise."); m.def("is_cuda_available", &gpu_offset_has_cuda, "Returns 1 if a functional CUDA GPU is present, 0 otherwise."); @@ -417,4 +761,58 @@ PYBIND11_MODULE(_cudaslicer_backend, m) { py::arg("offset_distances") = std::vector(), py::arg("backend") = "auto", "Mesh slicing, perimeter offsetting, and infill mask generation."); + + // GPU SDF Stack Factory Functions + m.def("slice_mesh_to_sdf", &py_slice_mesh_to_sdf, + py::arg("triangles_arr"), py::arg("min_x"), py::arg("min_y"), py::arg("min_z"), + py::arg("max_x"), py::arg("max_y"), py::arg("max_z"), + py::arg("layer_height_mm") = 0.20f, + py::arg("first_layer_height_mm") = 0.0f, + py::arg("pitch_mm") = 0.10f, + py::arg("pad_mm") = 1.5f, + "Slice 3D mesh directly to a GPU-resident Signed Distance Field stack."); + + m.def("masks_to_sdf", &py_masks_to_sdf, + py::arg("in_masks"), + py::arg("layer_z_um_arr") = py::array_t(), + py::arg("min_x_mm") = 0.0f, + py::arg("min_y_mm") = 0.0f, + py::arg("pitch_mm") = 0.10f, + "Compute GPU Signed Distance Field stack from 2D/3D binary raster masks."); + + // GPU SDF Stack Class Binding + py::class_>(m, "GpuSdfStack") + .def("get_metadata", &PyGpuSdfStack::get_metadata, "Get dictionary of spatial metadata.") + .def("get_layer_z_um", &PyGpuSdfStack::get_layer_z_um, "Get 1D array of layer heights in micrometers.") + .def("get_layer_z_mm", &PyGpuSdfStack::get_layer_z_mm, "Get 1D array of layer heights in millimeters.") + .def("to_numpy", &PyGpuSdfStack::to_numpy, "Download SDF float array (N, H, W) to CPU NumPy buffer.") + .def("clone", &PyGpuSdfStack::clone, "Clone this SDF stack on GPU.") + .def("csg_op", &PyGpuSdfStack::csg_op, + py::arg("other"), + py::arg("op_name") = "union", + py::arg("match_mode_name") = "dominant", + py::arg("blend_radius_mm") = 0.0f, + "Execute a GPU CSG operation with another SDF stack.") + .def("op_union", &PyGpuSdfStack::op_union, py::arg("other"), py::arg("blend_radius_mm") = 0.0f, "Compute A U B on GPU.") + .def("op_intersection", &PyGpuSdfStack::op_intersection, py::arg("other"), py::arg("blend_radius_mm") = 0.0f, "Compute A ^ B on GPU.") + .def("op_difference", &PyGpuSdfStack::op_difference, py::arg("other"), py::arg("blend_radius_mm") = 0.0f, "Compute A \\ B on GPU.") + .def("offset", &PyGpuSdfStack::offset, py::arg("delta_mm"), "Offset/dilate (delta > 0) or inset/erode (delta < 0) the SDF.") + .def("shell", &PyGpuSdfStack::shell, py::arg("thickness_mm"), "Hollow/shell the SDF with specified wall thickness.") + .def("generate_infill", &PyGpuSdfStack::generate_infill, + py::arg("pattern_str") = "gyroid", + py::arg("cell_size_mm") = 5.0f, + py::arg("density_ratio") = 0.20f, + py::arg("wall_offset_mm") = 1.0f, + "Intersect interior volume with 3D procedural infill lattice directly on GPU.") + .def("vectorize", &PyGpuSdfStack::vectorize, + py::arg("perimeters_count") = 2, + py::arg("layer_width") = 0.45f, + py::arg("export_infill_mask") = true, + py::arg("simplify_tolerance_mm") = 0.02f, + py::arg("subpixel_lerp") = true, + py::arg("offset_distances") = std::vector(), + "Vectorize SDF stack directly with sub-pixel LERP and RDP polygon simplification.") + .def_property_readonly("__cuda_array_interface__", &PyGpuSdfStack::cuda_array_interface, + "Zero-copy CUDA array interface for PyTorch / CuPy integration."); } + diff --git a/src/gpu_offset.cu b/src/gpu_offset.cu index e002878..fb09087 100644 --- a/src/gpu_offset.cu +++ b/src/gpu_offset.cu @@ -317,6 +317,39 @@ __global__ void gpu_parity_scan_to_mask_kernel( } } +// Helper function for direct mesh slicing to device mask +extern "C" int gpu_slice_mesh_to_device_mask( + const GpuTriangle3D* triangles, + size_t num_triangles, + float min_z, float h_first, float layer_height, int num_layers, + float bounds_min_x, float bounds_min_y, float pitch, + int cols, int rows, + uint8_t* d_masks, + cudaStream_t stream +) { + int row_stride_words = (cols + 31) / 32; + size_t total_words = (size_t)num_layers * rows * row_stride_words; + + unsigned int* d_parity = nullptr; + CUDA_CHECK(cudaMalloc(&d_parity, total_words * sizeof(unsigned int))); + CUDA_CHECK(cudaMemsetAsync(d_parity, 0, total_words * sizeof(unsigned int), stream)); + + int block_sz = 256; + int grid_sz = (num_triangles + block_sz - 1) / block_sz; + gpu_mesh_slice_raster_kernel<<>>( + triangles, num_triangles, min_z, h_first, layer_height, num_layers, + bounds_min_x, bounds_min_y, pitch, cols, rows, row_stride_words, d_parity + ); + + dim3 scan_grid(rows, num_layers); + gpu_parity_scan_to_mask_kernel<<>>( + d_parity, d_masks, cols, rows, num_layers, row_stride_words + ); + + cudaFree(d_parity); + return GPU_OFFSET_SUCCESS; +} + // ----------------------------------------------------------------------------- // Public C API Implementation & Backend Management // ----------------------------------------------------------------------------- diff --git a/src/gpu_sdf.cu b/src/gpu_sdf.cu new file mode 100644 index 0000000..020df77 --- /dev/null +++ b/src/gpu_sdf.cu @@ -0,0 +1,1184 @@ +#include "gpu_sdf.h" +#include "gpu_offset.h" +#include "cpu_offset.h" +#include "vectorizer.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define CUDA_CHECK(call) do { \ + cudaError_t err = call; \ + if (err != cudaSuccess) { \ + std::cerr << "[CUDA Error] " << cudaGetErrorString(err) \ + << " at " << __FILE__ << ":" << __LINE__ << std::endl; \ + return nullptr; \ + } \ +} while (0) + +#define CUDA_CHECK_INT(call) do { \ + cudaError_t err = call; \ + if (err != cudaSuccess) { \ + std::cerr << "[CUDA Error] " << cudaGetErrorString(err) \ + << " at " << __FILE__ << ":" << __LINE__ << std::endl; \ + return GPU_OFFSET_ERR_GENERIC; \ + } \ +} while (0) + +// Seed packing helpers +__device__ __inline__ int pack_seed_sdf(short x, short y) { + return ((int)(unsigned short)y << 16) | (int)(unsigned short)x; +} +__device__ __inline__ short get_seed_x_sdf(int packed) { + return (short)(packed & 0xFFFF); +} +__device__ __inline__ short get_seed_y_sdf(int packed) { + return (short)((packed >> 16) & 0xFFFF); +} + +// ----------------------------------------------------------------------------- +// Internal Implementation of GpuSdfStack +// ----------------------------------------------------------------------------- +struct GpuSdfStack_T { + GpuSdfMetadata meta; + int32_t* d_z_um; // GPU pointer to array of layer heights [num_layers] in micrometers + float* d_sdf; // GPU pointer to 3D continuous distance field [num_layers * height * width] in mm + cudaStream_t stream; + bool owns_stream; +}; + +struct LayerSamplingInfo { + int l0_a; + int l1_a; + float t_a; // < 0 if not present in A + int l0_b; + int l1_b; + float t_b; // < 0 if not present in B +}; + +// ----------------------------------------------------------------------------- +// CUDA Kernels for JFA Full Signed Distance Generation +// ----------------------------------------------------------------------------- +__global__ void jfa_init_interior_exterior_kernel( + const uint8_t* __restrict__ masks, + int* __restrict__ seeds_ext, + int* __restrict__ seeds_int, + int width, int height, int layers +) { + int x = blockIdx.x * blockDim.x + threadIdx.x; + int y = blockIdx.y * blockDim.y + threadIdx.y; + int z = blockIdx.z; + if (x >= width || y >= height || z >= layers) return; + + size_t idx = ((size_t)z * height + y) * width + x; + uint8_t m = masks[idx]; + + // seeds_ext: background pixels (0) are seeds to propagate into foreground (1) + if (m == 0) { + seeds_ext[idx] = pack_seed_sdf((short)x, (short)y); + } else { + seeds_ext[idx] = -1; + } + + // seeds_int: foreground pixels (1) are seeds to propagate into background (0) + if (m > 0) { + seeds_int[idx] = pack_seed_sdf((short)x, (short)y); + } else { + seeds_int[idx] = -1; + } +} + +__global__ void jfa_step_dual_kernel( + const int* __restrict__ in_seeds_ext, + int* __restrict__ out_seeds_ext, + const int* __restrict__ in_seeds_int, + int* __restrict__ out_seeds_int, + int width, int height, int layers, + int step +) { + int x = blockIdx.x * blockDim.x + threadIdx.x; + int y = blockIdx.y * blockDim.y + threadIdx.y; + int z = blockIdx.z; + if (x >= width || y >= height || z >= layers) return; + + size_t layer_offset = (size_t)z * height * width; + size_t idx = layer_offset + y * width + x; + + // 1. Exterior seeds + int best_seed_ext = in_seeds_ext[idx]; + int best_dist2_ext = 0x7FFFFFFF; + if (best_seed_ext != -1) { + int dx = (int)get_seed_x_sdf(best_seed_ext) - x; + int dy = (int)get_seed_y_sdf(best_seed_ext) - y; + best_dist2_ext = dx * dx + dy * dy; + } + + // 2. Interior seeds + int best_seed_int = in_seeds_int[idx]; + int best_dist2_int = 0x7FFFFFFF; + if (best_seed_int != -1) { + int dx = (int)get_seed_x_sdf(best_seed_int) - x; + int dy = (int)get_seed_y_sdf(best_seed_int) - y; + best_dist2_int = dx * dx + dy * dy; + } + + #pragma unroll + for (int dy = -1; dy <= 1; ++dy) { + #pragma unroll + for (int dx = -1; dx <= 1; ++dx) { + if (dx == 0 && dy == 0) continue; + int nx = x + dx * step; + int ny = y + dy * step; + if ((unsigned int)nx < (unsigned int)width && (unsigned int)ny < (unsigned int)height) { + size_t n_idx = layer_offset + ny * width + nx; + + int s_ext = in_seeds_ext[n_idx]; + if (s_ext != -1) { + int dist_x = (int)get_seed_x_sdf(s_ext) - x; + int dist_y = (int)get_seed_y_sdf(s_ext) - y; + int dist2 = dist_x * dist_x + dist_y * dist_y; + if (dist2 < best_dist2_ext) { + best_dist2_ext = dist2; + best_seed_ext = s_ext; + } + } + + int s_int = in_seeds_int[n_idx]; + if (s_int != -1) { + int dist_x = (int)get_seed_x_sdf(s_int) - x; + int dist_y = (int)get_seed_y_sdf(s_int) - y; + int dist2 = dist_x * dist_x + dist_y * dist_y; + if (dist2 < best_dist2_int) { + best_dist2_int = dist2; + best_seed_int = s_int; + } + } + } + } + } + + out_seeds_ext[idx] = best_seed_ext; + out_seeds_int[idx] = best_seed_int; +} + +__global__ void fuse_signed_sdf_kernel( + const uint8_t* __restrict__ masks, + const int* __restrict__ seeds_ext, + const int* __restrict__ seeds_int, + float* __restrict__ out_sdf, + int width, int height, int layers, + float pitch_mm +) { + int x = blockIdx.x * blockDim.x + threadIdx.x; + int y = blockIdx.y * blockDim.y + threadIdx.y; + int z = blockIdx.z; + if (x >= width || y >= height || z >= layers) return; + + size_t idx = ((size_t)z * height + y) * width + x; + uint8_t m = masks[idx]; + + if (m > 0) { + int s = seeds_ext[idx]; + if (s != -1) { + float dx = (float)(get_seed_x_sdf(s) - x); + float dy = (float)(get_seed_y_sdf(s) - y); + out_sdf[idx] = sqrtf(dx * dx + dy * dy) * pitch_mm; + } else { + out_sdf[idx] = 1e5f; + } + } else { + int s = seeds_int[idx]; + if (s != -1) { + float dx = (float)(get_seed_x_sdf(s) - x); + float dy = (float)(get_seed_y_sdf(s) - y); + out_sdf[idx] = -sqrtf(dx * dx + dy * dy) * pitch_mm; + } else { + out_sdf[idx] = -1e5f; + } + } +} + +// ----------------------------------------------------------------------------- +// CSG Mathematical Device Functions +// ----------------------------------------------------------------------------- +__device__ __inline__ float csg_eval_op(float a, float b, int op, float blend_k) { + switch (op) { + case GPU_SDF_OP_UNION: + return fmaxf(a, b); + case GPU_SDF_OP_INTERSECTION: + return fminf(a, b); + case GPU_SDF_OP_DIFFERENCE: + return fminf(a, -b); + case GPU_SDF_OP_SMOOTH_UNION: { + if (blend_k <= 1e-6f) return fmaxf(a, b); + float h = fmaxf(blend_k - fabsf(a - b), 0.0f) / blend_k; + return fmaxf(a, b) + h * h * blend_k * 0.25f; + } + case GPU_SDF_OP_SMOOTH_INTERSECT: { + if (blend_k <= 1e-6f) return fminf(a, b); + float h = fmaxf(blend_k - fabsf(a - b), 0.0f) / blend_k; + return fminf(a, b) - h * h * blend_k * 0.25f; + } + case GPU_SDF_OP_SMOOTH_DIFF: { + float neg_b = -b; + if (blend_k <= 1e-6f) return fminf(a, neg_b); + float h = fmaxf(blend_k - fabsf(a - neg_b), 0.0f) / blend_k; + return fminf(a, neg_b) - h * h * blend_k * 0.25f; + } + default: + return a; + } +} + +// ----------------------------------------------------------------------------- +// Identical Grid Fast CSG Kernel +// ----------------------------------------------------------------------------- +__global__ void gpu_sdf_csg_identical_kernel( + const float* __restrict__ sdf_a, + const float* __restrict__ sdf_b, + float* __restrict__ sdf_out, + size_t total_elements, + int op, + float blend_k +) { + size_t idx = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= total_elements) return; + + float a = sdf_a[idx]; + float b = sdf_b[idx]; + sdf_out[idx] = csg_eval_op(a, b, op, blend_k); +} + +// ----------------------------------------------------------------------------- +// Harmonized Resampling + CSG Kernel +// ----------------------------------------------------------------------------- +__device__ __inline__ float sample_stack_sdf( + const float* __restrict__ sdf, + int w, int h, + int x, int y, + int l0, int l1, float t +) { + if (t < -0.5f) return -1e5f; // Not present in stack (outside) + if (x < 0 || x >= w || y < 0 || y >= h) return -1e5f; + + size_t layer_stride = (size_t)w * h; + if (l0 == l1 || t <= 1e-6f) { + return sdf[(size_t)l0 * layer_stride + y * w + x]; + } else { + float v0 = sdf[(size_t)l0 * layer_stride + y * w + x]; + float v1 = sdf[(size_t)l1 * layer_stride + y * w + x]; + return (1.0f - t) * v0 + t * v1; + } +} + +__global__ void gpu_sdf_csg_resample_kernel( + const float* __restrict__ sdf_a, + GpuSdfMetadata meta_a, + const float* __restrict__ sdf_b, + GpuSdfMetadata meta_b, + float* __restrict__ sdf_out, + GpuSdfMetadata meta_out, + const LayerSamplingInfo* __restrict__ layer_map, + int op, + float blend_k +) { + int x_out = blockIdx.x * blockDim.x + threadIdx.x; + int y_out = blockIdx.y * blockDim.y + threadIdx.y; + int l_out = blockIdx.z; + + if (x_out >= meta_out.width || y_out >= meta_out.height || l_out >= meta_out.num_layers) return; + + LayerSamplingInfo info = layer_map[l_out]; + + // Compute world X, Y in micrometers + int32_t world_x_um = meta_out.min_x_um + x_out * meta_out.pitch_um; + int32_t world_y_um = meta_out.min_y_um + y_out * meta_out.pitch_um; + + // Sample from Stack A + int x_a = (world_x_um - meta_a.min_x_um) / meta_a.pitch_um; + int y_a = (world_y_um - meta_a.min_y_um) / meta_a.pitch_um; + float val_a = sample_stack_sdf(sdf_a, meta_a.width, meta_a.height, x_a, y_a, info.l0_a, info.l1_a, info.t_a); + + // Sample from Stack B + int x_b = (world_x_um - meta_b.min_x_um) / meta_b.pitch_um; + int y_b = (world_y_um - meta_b.min_y_um) / meta_b.pitch_um; + float val_b = sample_stack_sdf(sdf_b, meta_b.width, meta_b.height, x_b, y_b, info.l0_b, info.l1_b, info.t_b); + + float res = csg_eval_op(val_a, val_b, op, blend_k); + + size_t out_idx = ((size_t)l_out * meta_out.height + y_out) * meta_out.width + x_out; + sdf_out[out_idx] = res; +} + +// ----------------------------------------------------------------------------- +// Unary GPU Operations (Offset, Shell, Infill) +// ----------------------------------------------------------------------------- +__global__ void gpu_sdf_offset_kernel( + const float* __restrict__ in_sdf, + float* __restrict__ out_sdf, + size_t total_elements, + float delta_mm +) { + size_t idx = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= total_elements) return; + out_sdf[idx] = in_sdf[idx] + delta_mm; +} + +__global__ void gpu_sdf_shell_kernel( + const float* __restrict__ in_sdf, + float* __restrict__ out_sdf, + size_t total_elements, + float thickness_mm +) { + size_t idx = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= total_elements) return; + + float d = in_sdf[idx]; + float half_t = thickness_mm * 0.5f; + // Positive inside shell, negative outside / hollow center + out_sdf[idx] = half_t - fabsf(d - half_t); +} + +__global__ void gpu_sdf_infill_kernel( + const float* __restrict__ in_sdf, + float* __restrict__ out_sdf, + int width, int height, int num_layers, + int32_t min_x_um, int32_t min_y_um, int32_t pitch_um, + const int32_t* __restrict__ z_heights_um, + int pattern, + float cell_size_mm, + float density_ratio, + float wall_offset_mm +) { + int x = blockIdx.x * blockDim.x + threadIdx.x; + int y = blockIdx.y * blockDim.y + threadIdx.y; + int l = blockIdx.z; + + if (x >= width || y >= height || l >= num_layers) return; + + size_t idx = ((size_t)l * height + y) * width + x; + float dist = in_sdf[idx]; + + // Outside boundary or within shell wall offset -> retain wall distance + if (dist <= wall_offset_mm) { + out_sdf[idx] = dist; + return; + } + + float world_x = (float)(min_x_um + x * pitch_um) * 0.001f; + float world_y = (float)(min_y_um + y * pitch_um) * 0.001f; + float world_z = (float)(z_heights_um[l]) * 0.001f; + + float cell_sz = (cell_size_mm > 0.1f) ? cell_size_mm : 5.0f; + float freq = 6.28318530718f / cell_sz; + + float px = world_x * freq; + float py = world_y * freq; + float pz = world_z * freq; + + float field = 0.0f; + if (pattern == GPU_SDF_INFILL_GYROID) { + field = sinf(px) * cosf(py) + sinf(py) * cosf(pz) + sinf(pz) * cosf(px); + } else if (pattern == GPU_SDF_INFILL_SCHWARZ_P) { + field = cosf(px) + cosf(py) + cosf(pz); + } else if (pattern == GPU_SDF_INFILL_SCHWARZ_D) { + field = cosf(px)*cosf(py)*cosf(pz) - sinf(px)*sinf(py)*sinf(pz); + } else if (pattern == GPU_SDF_INFILL_GRID) { + float fx = fabsf(fmodf(world_x, cell_sz) - cell_sz * 0.5f); + float fy = fabsf(fmodf(world_y, cell_sz) - cell_sz * 0.5f); + field = fminf(fx, fy) / (cell_sz * 0.5f) - 0.5f; + } + + float threshold = (density_ratio - 0.5f) * 1.5f; + // Combine lattice with the infill bounding region + float lattice_sdf = (field >= threshold) ? 1.0f : -1.0f; + out_sdf[idx] = fminf(dist - wall_offset_mm, lattice_sdf); +} + +// ----------------------------------------------------------------------------- +// Direct Bitpack Vectorizer Kernel from SDF +// ----------------------------------------------------------------------------- +__global__ void gpu_sdf_threshold_bitpack_kernel( + const float* __restrict__ sdf, + unsigned int* __restrict__ out_bit_words, + int width, int height, int layers, + const float* __restrict__ thresholds_mm, + int num_offsets, + int row_stride_words +) { + int x_word = blockIdx.x * blockDim.x + threadIdx.x; + int y = blockIdx.y; + int z = blockIdx.z; + + if (x_word >= row_stride_words || y >= height || z >= layers) return; + + size_t in_layer_offset = ((size_t)z * height + y) * width; + size_t total_word_layer = (size_t)layers * height * row_stride_words; + size_t out_layer_offset = (size_t)z * height * row_stride_words + y * row_stride_words; + + int x_base = x_word * 32; + + for (int o = 0; o < num_offsets; ++o) { + float thresh = thresholds_mm[o]; + unsigned int packed = 0; + + #pragma unroll + for (int b = 0; b < 32; ++b) { + int x = x_base + b; + if (x < width) { + float d = sdf[in_layer_offset + x]; + if (d >= thresh) { + packed |= (1u << b); + } + } + } + + size_t out_idx = (size_t)o * total_word_layer + out_layer_offset + x_word; + out_bit_words[out_idx] = packed; + } +} + +// ----------------------------------------------------------------------------- +// Public C API Implementation +// ----------------------------------------------------------------------------- +extern "C" { + +GPU_SDF_API GpuSdfStackHandle gpu_sdf_create( + const GpuSdfMetadata* meta, + const int32_t* layer_z_heights_um +) { + if (!meta || meta->width <= 0 || meta->height <= 0 || meta->num_layers <= 0) return nullptr; + + GpuSdfStackHandle stack = new GpuSdfStack_T(); + stack->meta = *meta; + stack->owns_stream = true; + if (cudaStreamCreate(&stack->stream) != cudaSuccess) { + delete stack; + return nullptr; + } + + size_t total_voxels = (size_t)meta->num_layers * meta->height * meta->width; + size_t sdf_bytes = total_voxels * sizeof(float); + size_t z_bytes = meta->num_layers * sizeof(int32_t); + + if (cudaMalloc(&stack->d_sdf, sdf_bytes) != cudaSuccess || + cudaMalloc(&stack->d_z_um, z_bytes) != cudaSuccess) { + if (stack->d_sdf) cudaFree(stack->d_sdf); + if (stack->d_z_um) cudaFree(stack->d_z_um); + cudaStreamDestroy(stack->stream); + delete stack; + return nullptr; + } + + if (layer_z_heights_um) { + cudaMemcpyAsync(stack->d_z_um, layer_z_heights_um, z_bytes, cudaMemcpyHostToDevice, stack->stream); + } else { + // Default zero heights + cudaMemsetAsync(stack->d_z_um, 0, z_bytes, stack->stream); + } + cudaMemsetAsync(stack->d_sdf, 0, sdf_bytes, stack->stream); + + return stack; +} + +GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_masks( + const uint8_t* h_masks, + int width, + int height, + size_t num_layers, + const int32_t* layer_z_heights_um, + int32_t min_x_um, + int32_t min_y_um, + int32_t pitch_um +) { + if (!h_masks || width <= 0 || height <= 0 || num_layers == 0) return nullptr; + + GpuSdfMetadata meta; + meta.width = width; + meta.height = height; + meta.num_layers = num_layers; + meta.min_x_um = min_x_um; + meta.min_y_um = min_y_um; + meta.pitch_um = pitch_um > 0 ? pitch_um : 100; + + GpuSdfStackHandle stack = gpu_sdf_create(&meta, layer_z_heights_um); + if (!stack) return nullptr; + + size_t total_voxels = (size_t)num_layers * height * width; + size_t mask_bytes = total_voxels * sizeof(uint8_t); + size_t seed_bytes = total_voxels * sizeof(int); + + uint8_t* d_masks = nullptr; + int* d_seeds_ext_a = nullptr; + int* d_seeds_ext_b = nullptr; + int* d_seeds_int_a = nullptr; + int* d_seeds_int_b = nullptr; + + if (cudaMalloc(&d_masks, mask_bytes) != cudaSuccess || + cudaMalloc(&d_seeds_ext_a, seed_bytes) != cudaSuccess || + cudaMalloc(&d_seeds_ext_b, seed_bytes) != cudaSuccess || + cudaMalloc(&d_seeds_int_a, seed_bytes) != cudaSuccess || + cudaMalloc(&d_seeds_int_b, seed_bytes) != cudaSuccess) { + if (d_masks) cudaFree(d_masks); + if (d_seeds_ext_a) cudaFree(d_seeds_ext_a); + if (d_seeds_ext_b) cudaFree(d_seeds_ext_b); + if (d_seeds_int_a) cudaFree(d_seeds_int_a); + if (d_seeds_int_b) cudaFree(d_seeds_int_b); + gpu_sdf_destroy(stack); + return nullptr; + } + + cudaMemcpyAsync(d_masks, h_masks, mask_bytes, cudaMemcpyHostToDevice, stack->stream); + + dim3 block(16, 16, 1); + dim3 grid((width + 15) / 16, (height + 15) / 16, (unsigned int)num_layers); + + jfa_init_interior_exterior_kernel<<stream>>>( + d_masks, d_seeds_ext_a, d_seeds_int_a, width, height, (int)num_layers + ); + + int max_dim = std::max(width, height); + int step = 1; + while (step * 2 < max_dim) step *= 2; + + int* in_seeds_ext = d_seeds_ext_a; + int* out_seeds_ext = d_seeds_ext_b; + int* in_seeds_int = d_seeds_int_a; + int* out_seeds_int = d_seeds_int_b; + + while (step >= 1) { + jfa_step_dual_kernel<<stream>>>( + in_seeds_ext, out_seeds_ext, + in_seeds_int, out_seeds_int, + width, height, (int)num_layers, step + ); + std::swap(in_seeds_ext, out_seeds_ext); + std::swap(in_seeds_int, out_seeds_int); + step /= 2; + } + + // JFA+2 and JFA+1 corrections + jfa_step_dual_kernel<<stream>>>( + in_seeds_ext, out_seeds_ext, + in_seeds_int, out_seeds_int, + width, height, (int)num_layers, 2 + ); + std::swap(in_seeds_ext, out_seeds_ext); + std::swap(in_seeds_int, out_seeds_int); + + jfa_step_dual_kernel<<stream>>>( + in_seeds_ext, out_seeds_ext, + in_seeds_int, out_seeds_int, + width, height, (int)num_layers, 1 + ); + std::swap(in_seeds_ext, out_seeds_ext); + std::swap(in_seeds_int, out_seeds_int); + + float pitch_mm = (float)meta.pitch_um * 0.001f; + fuse_signed_sdf_kernel<<stream>>>( + d_masks, in_seeds_ext, in_seeds_int, stack->d_sdf, + width, height, (int)num_layers, pitch_mm + ); + + cudaStreamSynchronize(stack->stream); + + cudaFree(d_masks); + cudaFree(d_seeds_ext_a); + cudaFree(d_seeds_ext_b); + cudaFree(d_seeds_int_a); + cudaFree(d_seeds_int_b); + + return stack; +} + +// Forward declare GPU mesh rasterizer helper from gpu_offset.cu +extern int gpu_slice_mesh_to_device_mask( + const GpuTriangle3D* triangles, + size_t num_triangles, + float min_z, float h_first, float layer_height, int num_layers, + float bounds_min_x, float bounds_min_y, float pitch, + int cols, int rows, + uint8_t* d_masks, + cudaStream_t stream +); + +GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_mesh( + const GpuTriangle3D* triangles, + size_t num_triangles, + float min_x, float min_y, float min_z, + float max_x, float max_y, float max_z, + float layer_height_mm, + float first_layer_height_mm, + float pitch_mm, + float pad_mm +) { + if (!triangles || num_triangles == 0) return nullptr; + + float pad = (pad_mm > 0.0f) ? pad_mm : 1.5f; + float pitch = (pitch_mm > 0.0f) ? pitch_mm : 0.10f; + float bounds_min_x = min_x - pad; + float bounds_min_y = min_y - pad; + float bounds_max_x = max_x + pad; + float bounds_max_y = max_y + pad; + + int cols = (int)std::ceil((bounds_max_x - bounds_min_x) / pitch); + int rows = (int)std::ceil((bounds_max_y - bounds_min_y) / pitch); + + float h_first = (first_layer_height_mm > 0.0f) ? first_layer_height_mm : layer_height_mm; + float z_first = min_z + h_first; + + int num_layers = 0; + if (max_z >= z_first) { + num_layers = (int)std::floor((max_z - z_first) / layer_height_mm) + 1; + } + if (num_layers <= 0) return nullptr; + + std::vector z_heights_um(num_layers); + for (int l = 0; l < num_layers; ++l) { + float z_val = z_first + l * layer_height_mm; + z_heights_um[l] = (int32_t)std::round(z_val * 1000.0f); + } + + GpuSdfMetadata meta; + meta.width = cols; + meta.height = rows; + meta.num_layers = num_layers; + meta.min_x_um = (int32_t)std::round(bounds_min_x * 1000.0f); + meta.min_y_um = (int32_t)std::round(bounds_min_y * 1000.0f); + meta.pitch_um = (int32_t)std::round(pitch * 1000.0f); + + GpuSdfStackHandle stack = gpu_sdf_create(&meta, z_heights_um.data()); + if (!stack) return nullptr; + + size_t total_volume = (size_t)num_layers * rows * cols; + size_t tri_bytes = num_triangles * sizeof(GpuTriangle3D); + size_t seed_bytes = total_volume * sizeof(int); + + GpuTriangle3D* d_triangles = nullptr; + uint8_t* d_masks = nullptr; + int* d_seeds_ext_a = nullptr; + int* d_seeds_ext_b = nullptr; + int* d_seeds_int_a = nullptr; + int* d_seeds_int_b = nullptr; + + if (cudaMalloc(&d_triangles, tri_bytes) != cudaSuccess || + cudaMalloc(&d_masks, total_volume * sizeof(uint8_t)) != cudaSuccess || + cudaMalloc(&d_seeds_ext_a, seed_bytes) != cudaSuccess || + cudaMalloc(&d_seeds_ext_b, seed_bytes) != cudaSuccess || + cudaMalloc(&d_seeds_int_a, seed_bytes) != cudaSuccess || + cudaMalloc(&d_seeds_int_b, seed_bytes) != cudaSuccess) { + if (d_triangles) cudaFree(d_triangles); + if (d_masks) cudaFree(d_masks); + if (d_seeds_ext_a) cudaFree(d_seeds_ext_a); + if (d_seeds_ext_b) cudaFree(d_seeds_ext_b); + if (d_seeds_int_a) cudaFree(d_seeds_int_a); + if (d_seeds_int_b) cudaFree(d_seeds_int_b); + gpu_sdf_destroy(stack); + return nullptr; + } + + cudaMemcpyAsync(d_triangles, triangles, tri_bytes, cudaMemcpyHostToDevice, stack->stream); + + int slice_ret = gpu_slice_mesh_to_device_mask( + d_triangles, num_triangles, min_z, h_first, layer_height_mm, num_layers, + bounds_min_x, bounds_min_y, pitch, cols, rows, d_masks, stack->stream + ); + if (slice_ret != GPU_OFFSET_SUCCESS) { + cudaFree(d_triangles); + cudaFree(d_masks); + cudaFree(d_seeds_ext_a); + cudaFree(d_seeds_ext_b); + cudaFree(d_seeds_int_a); + cudaFree(d_seeds_int_b); + gpu_sdf_destroy(stack); + return nullptr; + } + + // Run Batched 3D JFA Dual Direction Signed Distance + dim3 block(16, 16, 1); + dim3 grid((cols + 15) / 16, (rows + 15) / 16, num_layers); + + jfa_init_interior_exterior_kernel<<stream>>>( + d_masks, d_seeds_ext_a, d_seeds_int_a, cols, rows, num_layers + ); + + int max_dim = std::max(cols, rows); + int step = 1; + while (step * 2 < max_dim) step *= 2; + + int* in_seeds_ext = d_seeds_ext_a; + int* out_seeds_ext = d_seeds_ext_b; + int* in_seeds_int = d_seeds_int_a; + int* out_seeds_int = d_seeds_int_b; + + while (step >= 1) { + jfa_step_dual_kernel<<stream>>>( + in_seeds_ext, out_seeds_ext, + in_seeds_int, out_seeds_int, + cols, rows, num_layers, step + ); + std::swap(in_seeds_ext, out_seeds_ext); + std::swap(in_seeds_int, out_seeds_int); + step /= 2; + } + + jfa_step_dual_kernel<<stream>>>( + in_seeds_ext, out_seeds_ext, + in_seeds_int, out_seeds_int, + cols, rows, num_layers, 2 + ); + std::swap(in_seeds_ext, out_seeds_ext); + std::swap(in_seeds_int, out_seeds_int); + + jfa_step_dual_kernel<<stream>>>( + in_seeds_ext, out_seeds_ext, + in_seeds_int, out_seeds_int, + cols, rows, num_layers, 1 + ); + std::swap(in_seeds_ext, out_seeds_ext); + std::swap(in_seeds_int, out_seeds_int); + + fuse_signed_sdf_kernel<<stream>>>( + d_masks, in_seeds_ext, in_seeds_int, stack->d_sdf, + cols, rows, num_layers, pitch + ); + + cudaStreamSynchronize(stack->stream); + + cudaFree(d_triangles); + cudaFree(d_masks); + cudaFree(d_seeds_ext_a); + cudaFree(d_seeds_ext_b); + cudaFree(d_seeds_int_a); + cudaFree(d_seeds_int_b); + + return stack; +} + +GPU_SDF_API GpuSdfStackHandle gpu_sdf_clone(GpuSdfStackHandle stack) { + if (!stack) return nullptr; + + std::vector z_um(stack->meta.num_layers); + gpu_sdf_get_layer_z_um(stack, z_um.data()); + + GpuSdfStackHandle clone = gpu_sdf_create(&stack->meta, z_um.data()); + if (!clone) return nullptr; + + size_t total_voxels = (size_t)stack->meta.num_layers * stack->meta.height * stack->meta.width; + cudaMemcpyAsync(clone->d_sdf, stack->d_sdf, total_voxels * sizeof(float), + cudaMemcpyDeviceToDevice, clone->stream); + cudaStreamSynchronize(clone->stream); + + return clone; +} + +GPU_SDF_API void gpu_sdf_destroy(GpuSdfStackHandle stack) { + if (!stack) return; + if (stack->d_sdf) { + cudaFree(stack->d_sdf); + stack->d_sdf = nullptr; + } + if (stack->d_z_um) { + cudaFree(stack->d_z_um); + stack->d_z_um = nullptr; + } + if (stack->owns_stream && stack->stream) { + cudaStreamDestroy(stack->stream); + stack->stream = nullptr; + } + delete stack; +} + +GPU_SDF_API int gpu_sdf_get_metadata(GpuSdfStackHandle stack, GpuSdfMetadata* out_meta) { + if (!stack || !out_meta) return GPU_OFFSET_ERR_INVALID_ARG; + *out_meta = stack->meta; + return GPU_OFFSET_SUCCESS; +} + +GPU_SDF_API int gpu_sdf_get_layer_z_um(GpuSdfStackHandle stack, int32_t* out_z_um) { + if (!stack || !out_z_um) return GPU_OFFSET_ERR_INVALID_ARG; + CUDA_CHECK_INT(cudaMemcpy(out_z_um, stack->d_z_um, stack->meta.num_layers * sizeof(int32_t), cudaMemcpyDeviceToHost)); + return GPU_OFFSET_SUCCESS; +} + +GPU_SDF_API int gpu_sdf_to_host(GpuSdfStackHandle stack, float* h_out_data) { + if (!stack || !h_out_data) return GPU_OFFSET_ERR_INVALID_ARG; + size_t total_voxels = (size_t)stack->meta.num_layers * stack->meta.height * stack->meta.width; + CUDA_CHECK_INT(cudaMemcpy(h_out_data, stack->d_sdf, total_voxels * sizeof(float), cudaMemcpyDeviceToHost)); + return GPU_OFFSET_SUCCESS; +} + +GPU_SDF_API float* gpu_sdf_get_device_ptr(GpuSdfStackHandle stack) { + return stack ? stack->d_sdf : nullptr; +} + +GPU_SDF_API int32_t* gpu_sdf_get_device_z_ptr(GpuSdfStackHandle stack) { + return stack ? stack->d_z_um : nullptr; +} + +// ----------------------------------------------------------------------------- +// CSG Execution Engine +// ----------------------------------------------------------------------------- +GPU_SDF_API GpuSdfStackHandle gpu_sdf_csg_op( + GpuSdfStackHandle a, + GpuSdfStackHandle b, + GpuSdfOp op, + GpuSdfLayerMatchMode match_mode, + float blend_radius_mm +) { + if (!a || !b) return nullptr; + + // Check if grids and layers are 100% identical + std::vector z_a(a->meta.num_layers); + std::vector z_b(b->meta.num_layers); + gpu_sdf_get_layer_z_um(a, z_a.data()); + gpu_sdf_get_layer_z_um(b, z_b.data()); + + bool identical_grid = (a->meta.min_x_um == b->meta.min_x_um) && + (a->meta.min_y_um == b->meta.min_y_um) && + (a->meta.pitch_um == b->meta.pitch_um) && + (a->meta.width == b->meta.width) && + (a->meta.height == b->meta.height) && + (a->meta.num_layers == b->meta.num_layers) && + (z_a == z_b); + + if (identical_grid) { + GpuSdfStackHandle out_stack = gpu_sdf_create(&a->meta, z_a.data()); + if (!out_stack) return nullptr; + + size_t total_voxels = (size_t)a->meta.num_layers * a->meta.height * a->meta.width; + int block_sz = 256; + int grid_sz = (total_voxels + block_sz - 1) / block_sz; + + gpu_sdf_csg_identical_kernel<<stream>>>( + a->d_sdf, b->d_sdf, out_stack->d_sdf, total_voxels, (int)op, blend_radius_mm + ); + cudaStreamSynchronize(out_stack->stream); + return out_stack; + } + + // Harmonize Spatial Grid (XY) + int32_t pitch_um = a->meta.pitch_um; + int32_t min_x_um = std::min(a->meta.min_x_um, b->meta.min_x_um); + int32_t min_y_um = std::min(a->meta.min_y_um, b->meta.min_y_um); + int32_t max_x_um = std::max(a->meta.min_x_um + a->meta.width * a->meta.pitch_um, + b->meta.min_x_um + b->meta.width * b->meta.pitch_um); + int32_t max_y_um = std::max(a->meta.min_y_um + a->meta.height * a->meta.pitch_um, + b->meta.min_y_um + b->meta.height * b->meta.pitch_um); + + int out_w = (int)std::ceil((float)(max_x_um - min_x_um) / pitch_um); + int out_h = (int)std::ceil((float)(max_y_um - min_y_um) / pitch_um); + + // Harmonize Layers (Z) + std::vector out_z_um; + if (match_mode == GPU_SDF_LAYER_MATCH_A_DOMINANT) { + out_z_um = z_a; + } else if (match_mode == GPU_SDF_LAYER_MATCH_UNION) { + std::set z_set(z_a.begin(), z_a.end()); + z_set.insert(z_b.begin(), z_b.end()); + out_z_um.assign(z_set.begin(), z_set.end()); + } else { // INTERSECTION + std::set set_b(z_b.begin(), z_b.end()); + for (int32_t z : z_a) { + if (set_b.count(z)) out_z_um.push_back(z); + } + } + + if (out_z_um.empty()) return nullptr; + + GpuSdfMetadata meta_out; + meta_out.width = out_w; + meta_out.height = out_h; + meta_out.num_layers = out_z_um.size(); + meta_out.min_x_um = min_x_um; + meta_out.min_y_um = min_y_um; + meta_out.pitch_um = pitch_um; + + GpuSdfStackHandle out_stack = gpu_sdf_create(&meta_out, out_z_um.data()); + if (!out_stack) return nullptr; + + // Build Layer Sampling Info table on Host & upload to Device + std::vector layer_infos(meta_out.num_layers); + + auto find_layer_sampling = [](const std::vector& z_list, int32_t target_z, int& l0, int& l1, float& t) { + if (z_list.empty()) { + l0 = l1 = -1; t = -1.0f; return; + } + if (target_z < z_list.front() || target_z > z_list.back()) { + l0 = l1 = -1; t = -1.0f; return; + } + auto it = std::lower_bound(z_list.begin(), z_list.end(), target_z); + if (it != z_list.end() && *it == target_z) { + l0 = l1 = (int)(it - z_list.begin()); + t = 0.0f; + return; + } + // Interpolation between adjacent layers + size_t idx1 = it - z_list.begin(); + size_t idx0 = idx1 - 1; + l0 = (int)idx0; + l1 = (int)idx1; + t = (float)(target_z - z_list[idx0]) / (float)(z_list[idx1] - z_list[idx0]); + }; + + for (size_t l = 0; l < meta_out.num_layers; ++l) { + int32_t tz = out_z_um[l]; + find_layer_sampling(z_a, tz, layer_infos[l].l0_a, layer_infos[l].l1_a, layer_infos[l].t_a); + find_layer_sampling(z_b, tz, layer_infos[l].l0_b, layer_infos[l].l1_b, layer_infos[l].t_b); + } + + LayerSamplingInfo* d_layer_map = nullptr; + cudaMalloc(&d_layer_map, meta_out.num_layers * sizeof(LayerSamplingInfo)); + cudaMemcpyAsync(d_layer_map, layer_infos.data(), meta_out.num_layers * sizeof(LayerSamplingInfo), + cudaMemcpyHostToDevice, out_stack->stream); + + dim3 block(16, 16, 1); + dim3 grid((out_w + 15) / 16, (out_h + 15) / 16, (unsigned int)meta_out.num_layers); + + gpu_sdf_csg_resample_kernel<<stream>>>( + a->d_sdf, a->meta, + b->d_sdf, b->meta, + out_stack->d_sdf, meta_out, + d_layer_map, (int)op, blend_radius_mm + ); + + cudaStreamSynchronize(out_stack->stream); + cudaFree(d_layer_map); + + return out_stack; +} + +GPU_SDF_API GpuSdfStackHandle gpu_sdf_offset( + GpuSdfStackHandle stack, + float delta_mm +) { + if (!stack) return nullptr; + + std::vector z_um(stack->meta.num_layers); + gpu_sdf_get_layer_z_um(stack, z_um.data()); + + GpuSdfStackHandle out_stack = gpu_sdf_create(&stack->meta, z_um.data()); + if (!out_stack) return nullptr; + + size_t total_voxels = (size_t)stack->meta.num_layers * stack->meta.height * stack->meta.width; + int block_sz = 256; + int grid_sz = (total_voxels + block_sz - 1) / block_sz; + + gpu_sdf_offset_kernel<<stream>>>( + stack->d_sdf, out_stack->d_sdf, total_voxels, delta_mm + ); + cudaStreamSynchronize(out_stack->stream); + + return out_stack; +} + +GPU_SDF_API GpuSdfStackHandle gpu_sdf_shell( + GpuSdfStackHandle stack, + float thickness_mm +) { + if (!stack || thickness_mm <= 0.0f) return nullptr; + + std::vector z_um(stack->meta.num_layers); + gpu_sdf_get_layer_z_um(stack, z_um.data()); + + GpuSdfStackHandle out_stack = gpu_sdf_create(&stack->meta, z_um.data()); + if (!out_stack) return nullptr; + + size_t total_voxels = (size_t)stack->meta.num_layers * stack->meta.height * stack->meta.width; + int block_sz = 256; + int grid_sz = (total_voxels + block_sz - 1) / block_sz; + + gpu_sdf_shell_kernel<<stream>>>( + stack->d_sdf, out_stack->d_sdf, total_voxels, thickness_mm + ); + cudaStreamSynchronize(out_stack->stream); + + return out_stack; +} + +GPU_SDF_API GpuSdfStackHandle gpu_sdf_generate_infill( + GpuSdfStackHandle stack, + GpuSdfInfillPattern pattern, + float cell_size_mm, + float density_ratio, + float wall_offset_mm +) { + if (!stack) return nullptr; + + std::vector z_um(stack->meta.num_layers); + gpu_sdf_get_layer_z_um(stack, z_um.data()); + + GpuSdfStackHandle out_stack = gpu_sdf_create(&stack->meta, z_um.data()); + if (!out_stack) return nullptr; + + dim3 block(16, 16, 1); + dim3 grid((stack->meta.width + 15) / 16, (stack->meta.height + 15) / 16, (unsigned int)stack->meta.num_layers); + + gpu_sdf_infill_kernel<<stream>>>( + stack->d_sdf, out_stack->d_sdf, + stack->meta.width, stack->meta.height, (int)stack->meta.num_layers, + stack->meta.min_x_um, stack->meta.min_y_um, stack->meta.pitch_um, + stack->d_z_um, (int)pattern, cell_size_mm, density_ratio, wall_offset_mm + ); + cudaStreamSynchronize(out_stack->stream); + + return out_stack; +} + +GPU_SDF_API int gpu_sdf_vectorize( + GpuSdfStackHandle stack, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + size_t* out_num_layers, + GpuOffsetStats* stats +) { + if (!stack || !config || !out_results) return GPU_OFFSET_ERR_INVALID_ARG; + + auto t_total_start = std::chrono::high_resolution_clock::now(); + + int cols = stack->meta.width; + int rows = stack->meta.height; + size_t num_layers = stack->meta.num_layers; + if (out_num_layers) *out_num_layers = num_layers; + + float bounds_min_x = (float)stack->meta.min_x_um * 0.001f; + float bounds_min_y = (float)stack->meta.min_y_um * 0.001f; + float pitch_mm = (float)stack->meta.pitch_um * 0.001f; + + int P = (config->perimeters_count >= 0) ? config->perimeters_count : 2; + float W = (config->layer_width > 0.0f) ? config->layer_width : 0.45f; + std::vector resolved_offsets; + if (config->offset_distances_mm && config->num_offsets > 0) { + resolved_offsets.assign(config->offset_distances_mm, config->offset_distances_mm + config->num_offsets); + P = (config->perimeters_count > 0) ? config->perimeters_count : (int)resolved_offsets.size() - 1; + if (P < 1) P = (int)resolved_offsets.size(); + } else { + resolved_offsets.resize(P + 1); + for (int i = 0; i < P; ++i) { + resolved_offsets[i] = (i + 0.5f) * W; + } + resolved_offsets[P] = (float)P * W; + } + size_t num_offsets = resolved_offsets.size(); + + float simplify_tol = config->simplify_tolerance_mm; + if (simplify_tol <= 0.0f && config->subpixel_lerp != 0) { + simplify_tol = 0.02f; // default 0.02 mm (20 µm) tolerance + } + + size_t total_voxels = (size_t)num_layers * rows * cols; + size_t sdf_bytes = total_voxels * sizeof(float); + + float* h_pinned_sdf = nullptr; + CUDA_CHECK_INT(cudaMallocHost(&h_pinned_sdf, sdf_bytes)); + CUDA_CHECK_INT(cudaMemcpyAsync(h_pinned_sdf, stack->d_sdf, sdf_bytes, + cudaMemcpyDeviceToHost, stack->stream)); + CUDA_CHECK_INT(cudaStreamSynchronize(stack->stream)); + + std::vector z_heights_um(num_layers); + gpu_sdf_get_layer_z_um(stack, z_heights_um.data()); + + // Multi-threaded C++ Sub-Pixel LERP Vectorization + auto t_vec_0 = std::chrono::high_resolution_clock::now(); + + GpuLayerResult* res_array = (GpuLayerResult*)malloc(num_layers * sizeof(GpuLayerResult)); + memset(res_array, 0, num_layers * sizeof(GpuLayerResult)); + + unsigned int num_threads = std::max(1u, std::thread::hardware_concurrency()); + std::vector workers; + + auto vec_worker = [&](size_t start_layer, size_t end_layer) { + int stride_w = cols + 1; + size_t max_edges = 2 * (size_t)(rows + 1) * stride_w; + std::vector next_edge_buf(max_edges, -1); + std::vector active_edges_buf; + active_edges_buf.reserve(cols * rows / 16); + std::vector visited_buf(max_edges, 0); + + for (size_t l = start_layer; l < end_layer; ++l) { + res_array[l].z_height = (float)z_heights_um[l] * 0.001f; + res_array[l].num_perimeters = P; + res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter)); + + const float* layer_sdf_ptr = h_pinned_sdf + l * ((size_t)rows * cols); + + for (int o = 0; o < P; ++o) { + float iso_thresh = resolved_offsets[o]; + + std::vector contours; + vectorize_sdf_layer(layer_sdf_ptr, cols, rows, bounds_min_x, bounds_min_y, pitch_mm, + iso_thresh, simplify_tol, + next_edge_buf, active_edges_buf, visited_buf, contours); + + res_array[l].perimeters[o].offset_distance_mm = iso_thresh; + res_array[l].perimeters[o].num_contours = contours.size(); + if (!contours.empty()) { + res_array[l].perimeters[o].contours = (GpuContour*)malloc(contours.size() * sizeof(GpuContour)); + memcpy(res_array[l].perimeters[o].contours, contours.data(), contours.size() * sizeof(GpuContour)); + } else { + res_array[l].perimeters[o].contours = nullptr; + } + } + + // Infill + float infill_iso_thresh = resolved_offsets[P]; + std::vector infill_contours; + vectorize_sdf_layer(layer_sdf_ptr, cols, rows, bounds_min_x, bounds_min_y, pitch_mm, + infill_iso_thresh, simplify_tol, + next_edge_buf, active_edges_buf, visited_buf, infill_contours); + + res_array[l].infill.offset_distance_mm = infill_iso_thresh; + res_array[l].infill.width = cols; + res_array[l].infill.height = rows; + res_array[l].infill.num_contours = infill_contours.size(); + if (!infill_contours.empty()) { + res_array[l].infill.contours = (GpuContour*)malloc(infill_contours.size() * sizeof(GpuContour)); + memcpy(res_array[l].infill.contours, infill_contours.data(), infill_contours.size() * sizeof(GpuContour)); + } else { + res_array[l].infill.contours = nullptr; + } + + if (config->export_infill_mask) { + uint8_t* infill_mask = (uint8_t*)malloc((size_t)cols * rows); + for (int i = 0; i < cols * rows; ++i) { + infill_mask[i] = (layer_sdf_ptr[i] >= infill_iso_thresh) ? 1 : 0; + } + res_array[l].infill.mask = infill_mask; + } else { + res_array[l].infill.mask = nullptr; + } + } + }; + + size_t chunk = (num_layers + num_threads - 1) / num_threads; + for (unsigned int t = 0; t < num_threads; ++t) { + size_t s = t * chunk; + size_t e = std::min(s + chunk, num_layers); + if (s < e) workers.emplace_back(vec_worker, s, e); + } + for (auto& w : workers) w.join(); + + auto t_vec_1 = std::chrono::high_resolution_clock::now(); + double t_vec_ms = std::chrono::duration(t_vec_1 - t_vec_0).count(); + + auto t_total_end = std::chrono::high_resolution_clock::now(); + double t_total_ms = std::chrono::duration(t_total_end - t_total_start).count(); + + if (stats) { + stats->t_vectorize_ms = t_vec_ms; + stats->t_total_pipeline_ms = t_total_ms; + } + + cudaFreeHost(h_pinned_sdf); + + *out_results = res_array; + return GPU_OFFSET_SUCCESS; +} + +} // extern "C" diff --git a/src/vectorizer.h b/src/vectorizer.h index b77a865..99b9cf6 100644 --- a/src/vectorizer.h +++ b/src/vectorizer.h @@ -28,6 +28,275 @@ inline GpuPoint2D get_edge_point(int edge_id, int stride_w, float min_x, float m } } +// Sub-pixel linear interpolation on continuous distance field +inline GpuPoint2D get_sdf_edge_point( + int edge_id, + const float* sdf, + int width, int height, + float min_x, float min_y, float pitch_mm, + float iso_val +) { + int type = edge_id & 1; + int cell = edge_id >> 1; + int stride_w = width + 1; + int y = cell / stride_w; + int x = cell % stride_w; + + if (type == 0) { + // Horizontal edge: from (x, y) to (x + 1, y) + float d0 = -1e5f; + float d1 = -1e5f; + if (y >= 0 && y < height) { + if (x >= 0 && x < width) d0 = sdf[y * width + x]; + if (x + 1 >= 0 && x + 1 < width) d1 = sdf[y * width + (x + 1)]; + } + float t = 0.5f; + float denom = d1 - d0; + if (std::abs(denom) > 1e-7f) { + t = (iso_val - d0) / denom; + if (t < 0.0f) t = 0.0f; + else if (t > 1.0f) t = 1.0f; + } + return GpuPoint2D{ min_x + ((float)x + t) * pitch_mm, min_y + (float)y * pitch_mm }; + } else { + // Vertical edge: from (x, y) to (x, y + 1) + float d0 = -1e5f; + float d3 = -1e5f; + if (x >= 0 && x < width) { + if (y >= 0 && y < height) d0 = sdf[y * width + x]; + if (y + 1 >= 0 && y + 1 < height) d3 = sdf[(y + 1) * width + x]; + } + float t = 0.5f; + float denom = d3 - d0; + if (std::abs(denom) > 1e-7f) { + t = (iso_val - d0) / denom; + if (t < 0.0f) t = 0.0f; + else if (t > 1.0f) t = 1.0f; + } + return GpuPoint2D{ min_x + (float)x * pitch_mm, min_y + ((float)y + t) * pitch_mm }; + } +} + +// Ramer-Douglas-Peucker (RDP) Polyline Simplification +static inline void rdp_simplify_recursive( + const std::vector& pts, + size_t i, size_t j, + float eps_sq, + std::vector& keep +) { + if (j <= i + 1) return; + + float p1x = pts[i].x, p1y = pts[i].y; + float p2x = pts[j].x, p2y = pts[j].y; + float dx = p2x - p1x; + float dy = p2y - p1y; + float len_sq = dx * dx + dy * dy; + + float max_dist_sq = 0.0f; + size_t max_idx = i; + + for (size_t k = i + 1; k < j; ++k) { + float px = pts[k].x, py = pts[k].y; + float dist_sq = 0.0f; + + if (len_sq < 1e-10f) { + float dpx = px - p1x; + float dpy = py - p1y; + dist_sq = dpx * dpx + dpy * dpy; + } else { + float cross = dy * (px - p1x) - dx * (py - p1y); + dist_sq = (cross * cross) / len_sq; + } + + if (dist_sq > max_dist_sq) { + max_dist_sq = dist_sq; + max_idx = k; + } + } + + if (max_dist_sq > eps_sq) { + keep[max_idx] = true; + rdp_simplify_recursive(pts, i, max_idx, eps_sq, keep); + rdp_simplify_recursive(pts, max_idx, j, eps_sq, keep); + } +} + +static inline std::vector simplify_contour_rdp( + const std::vector& loop_pts, + float tolerance_mm +) { + size_t n = loop_pts.size(); + if (n < 4 || tolerance_mm <= 1e-6f) return loop_pts; + + float eps_sq = tolerance_mm * tolerance_mm; + + // For closed loop: find point furthest from loop_pts[0] + float max_d_sq = 0.0f; + size_t split_idx = n / 2; + for (size_t i = 1; i < n; ++i) { + float dx = loop_pts[i].x - loop_pts[0].x; + float dy = loop_pts[i].y - loop_pts[0].y; + float d_sq = dx * dx + dy * dy; + if (d_sq > max_d_sq) { + max_d_sq = d_sq; + split_idx = i; + } + } + + std::vector keep(n, false); + keep[0] = true; + keep[split_idx] = true; + + rdp_simplify_recursive(loop_pts, 0, split_idx, eps_sq, keep); + rdp_simplify_recursive(loop_pts, split_idx, n - 1, eps_sq, keep); + + std::vector simplified; + simplified.reserve(n); + for (size_t i = 0; i < n; ++i) { + if (keep[i]) simplified.push_back(loop_pts[i]); + } + return simplified; +} + +static inline void vectorize_sdf_layer( + const float* sdf, + int width, int height, + float min_x, float min_y, float pitch_mm, + float iso_val, + float simplify_tolerance_mm, + std::vector& next_edge_buf, + std::vector& active_edges_buf, + std::vector& visited_buf, + std::vector& out_contours +) { + int stride_w = width + 1; + active_edges_buf.clear(); + + auto add_directed_edge = [&](int e0, int e1) { + next_edge_buf[e0] = e1; + visited_buf[e0] = 0; + active_edges_buf.push_back(e0); + }; + + for (int y = 0; y < height - 1; ++y) { + const float* row0 = sdf + y * width; + const float* row1 = sdf + (y + 1) * width; + + for (int x = 0; x < width - 1; ++x) { + float c0 = row0[x]; + float c1 = row0[x + 1]; + float c2 = row1[x + 1]; + float c3 = row1[x]; + + int b0 = (c0 >= iso_val) ? 1 : 0; + int b1 = (c1 >= iso_val) ? 1 : 0; + int b2 = (c2 >= iso_val) ? 1 : 0; + int b3 = (c3 >= iso_val) ? 1 : 0; + + int case_idx = b0 | (b1 << 1) | (b2 << 2) | (b3 << 3); + if (case_idx == 0 || case_idx == 15) continue; + + int e_bottom = get_h_edge_id(x, y, stride_w); + int e_top = get_h_edge_id(x, y + 1, stride_w); + int e_left = get_v_edge_id(x, y, stride_w); + int e_right = get_v_edge_id(x + 1, y, stride_w); + + switch (case_idx) { + case 1: add_directed_edge(e_left, e_bottom); break; + case 2: add_directed_edge(e_bottom, e_right); break; + case 3: add_directed_edge(e_left, e_right); break; + case 4: add_directed_edge(e_right, e_top); break; + case 5: { + float center = (c0 + c1 + c2 + c3) * 0.25f; + if (center >= iso_val) { + add_directed_edge(e_left, e_top); + add_directed_edge(e_bottom, e_right); + } else { + add_directed_edge(e_left, e_bottom); + add_directed_edge(e_right, e_top); + } + break; + } + case 6: add_directed_edge(e_bottom, e_top); break; + case 7: add_directed_edge(e_left, e_top); break; + case 8: add_directed_edge(e_top, e_left); break; + case 9: add_directed_edge(e_top, e_bottom); break; + case 10: { + float center = (c0 + c1 + c2 + c3) * 0.25f; + if (center >= iso_val) { + add_directed_edge(e_bottom, e_left); + add_directed_edge(e_top, e_right); + } else { + add_directed_edge(e_top, e_left); + add_directed_edge(e_bottom, e_right); + } + break; + } + case 11: add_directed_edge(e_top, e_right); break; + case 12: add_directed_edge(e_right, e_left); break; + case 13: add_directed_edge(e_right, e_bottom); break; + case 14: add_directed_edge(e_bottom, e_left); break; + } + } + } + + if (active_edges_buf.empty()) return; + + std::vector loop_pts; + loop_pts.reserve(1000); + + for (int start_e : active_edges_buf) { + if (visited_buf[start_e]) continue; + + loop_pts.clear(); + int cur = start_e; + + while (cur != -1 && !visited_buf[cur]) { + visited_buf[cur] = 1; + GpuPoint2D pt = get_sdf_edge_point(cur, sdf, width, height, min_x, min_y, pitch_mm, iso_val); + + size_t n = loop_pts.size(); + if (n >= 2) { + float dx1 = loop_pts[n-1].x - loop_pts[n-2].x; + float dy1 = loop_pts[n-1].y - loop_pts[n-2].y; + float dx2 = pt.x - loop_pts[n-1].x; + float dy2 = pt.y - loop_pts[n-1].y; + + float cross = dx1 * dy2 - dy1 * dx2; + if (std::abs(cross) < 1e-7f) { + loop_pts[n-1] = pt; // simplify collinear vertex + } else { + loop_pts.push_back(pt); + } + } else { + loop_pts.push_back(pt); + } + + cur = next_edge_buf[cur]; + if (cur == start_e) break; + if (loop_pts.size() > 50000) break; + } + + if (loop_pts.size() >= 3) { + std::vector simplified = simplify_contour_rdp(loop_pts, simplify_tolerance_mm); + if (simplified.size() >= 3) { + GpuContour c; + c.num_points = simplified.size(); + c.points = (GpuPoint2D*)malloc(c.num_points * sizeof(GpuPoint2D)); + memcpy(c.points, simplified.data(), c.num_points * sizeof(GpuPoint2D)); + c.is_hole = 0; + out_contours.push_back(c); + } + } + } + + // Reset thread-local edge table for zero re-allocation overhead + for (int e : active_edges_buf) { + next_edge_buf[e] = -1; + visited_buf[e] = 0; + } +} + static inline void vectorize_mask_layer( const uint8_t* mask, int width, int height,