diff --git a/CMakeLists.txt b/CMakeLists.txt index b3d1485..1ab357c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 3.18) -project(cudaslicer LANGUAGES C CXX CUDA) +project(cudaslicer VERSION 1.0.0 LANGUAGES C CXX CUDA) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) diff --git a/Makefile b/Makefile index 2bc22c2..b8b1238 100644 --- a/Makefile +++ b/Makefile @@ -45,4 +45,7 @@ $(BUILD_DIR) $(LIB_DIR) $(BIN_DIR): clean: rm -rf $(BUILD_DIR) $(LIB_DIR) $(BIN_DIR) -.PHONY: all clean +wheel: + python3 setup.py bdist_wheel + +.PHONY: all clean wheel diff --git a/README.md b/README.md index 551bf9e..5a47845 100644 --- a/README.md +++ b/README.md @@ -141,6 +141,80 @@ perims, infill_mask = layer_0 print(f"Layer 0 (Z = {layer_0.z} mm): {len(perims)} perimeters, infill mask: {infill_mask.shape}") ``` +#### E. GPU Signed Distance Fields (SDF), Boolean CSG & Gyroid Infill +Execute boolean CSG operations (Union, Difference, Intersection, Smooth Blending) and 3D procedural infill directly on GPU SDF stacks with zero CPU memory roundtrips: + +```python +import trimesh +import cudaslicer + +mesh_a = trimesh.load("benchy.stl") +cutter_box = trimesh.creation.box(extents=(12.0, 60.0, 25.0)) + +# 1. Direct GPU slice to continuous SDF stacks +sdf_a = cudaslicer.slice_mesh_to_sdf(mesh_a, layer_height=0.20, pitch_mm=0.10) +sdf_b = cudaslicer.slice_mesh_to_sdf(cutter_box, layer_height=0.20, pitch_mm=0.10) + +# 2. Boolean CSG subtraction entirely in GPU VRAM +sdf_carved = sdf_a.difference(sdf_b.offset(delta_mm=0.5)) + +# 3. GPU procedural 3D gyroid infill generation +sdf_infill = sdf_carved.generate_infill(pattern="gyroid", cell_size_mm=4.0, density=0.25) + +# 4. Direct vectorization with sub-pixel LERP and RDP polyline simplification +layers = sdf_infill.vectorize( + perimeters_count=2, + layer_width=0.45, + subpixel_lerp=True, # Natural arbitrary angles (no 45° quantization) + simplify_tolerance_mm=0.02 # RDP polyline decimation +) +``` + +--- + +## 📐 Sub-Pixel Linear Interpolation (LERP) vs. Standard Marching Squares + +### The Problem: 45° & 90° Angle Quantization in Standard Marching Squares +Standard discrete Marching Squares (and binary mask polygon extractors) evaluates binary cell states (0 or 1 at each voxel corner). When an edge intersects a boundary, it snaps the intersection point to the **geometric midpoint** of the grid cell edge ($t = 0.5$): + +``` +Discrete Marching Squares (Snap to t = 0.5): +(x, y+1) +---------+ (x+1, y+1) + | * | <- Midpoint snapping forces line segments to connect + | / | only at (0.5, 0), (0, 0.5), (1, 0.5), etc. + | * | <- Result: angles restricted to 0°, 45°, 90°, 135°! +(x, y) +---------+ (x+1, y) +``` + +Because segment vectors between grid edge midpoints are strictly constrained to Cartesian diagonals, diagonal walls and curves are quantized into jagged 45° chamfers and stair-stepped facets. + +--- + +### The Solution: Continuous SDF Linear Interpolation (LERP) +In CUDASlicer, the vectorizer uses continuous floating-point Signed Distance Field (SDF) values ($d_0, d_1$) at the cell vertices to solve for the exact sub-pixel boundary crossing where the distance field equals the target isovalue ($iso\_val = 0.0$ for outer perimeter, $-d$ for offset contours): + +$$t = \frac{iso\_val - d_0}{d_1 - d_0}, \quad t \in [0.0, 1.0]$$ + +``` +Continuous SDF LERP (Sub-Pixel Placement): +(x, y+1) +---------+ (x+1, y+1) + | | + | * | <- Exact root crossing t = (0 - d0) / (d1 - d0) + | / | <- Arbitrary natural angles (e.g. 12.3°, 37.8°, smooth curves) +(x, y) +-----*---+ (x+1, y) +``` + +- **Horizontal Edge:** $P_x = \text{min\_x} + (x + t) \cdot \text{pitch}, \quad P_y = \text{min\_y} + y \cdot \text{pitch}$ +- **Vertical Edge:** $P_x = \text{min\_x} + x \cdot \text{pitch}, \quad P_y = \text{min\_y} + (y + t) \cdot \text{pitch}$ + +Since $t$ varies continuously in $[0, 1]$, vertices fall at their true Euclidean geometric boundaries, eliminating jagged 45° artifacts. + +--- + +### Polyline Simplification & Decimation +1. **On-the-Fly Collinear Decimation:** During loop traversal, collinear vertices ($\Delta x_1 \Delta y_2 - \Delta y_1 \Delta x_2 \approx 0$) along straight spans are pruned on the fly with zero memory overhead. +2. **Ramer-Douglas-Peucker (RDP) Polyline Simplification:** Traversed closed loops are post-processed with `simplify_contour_rdp(loop, tolerance_mm)` (default tolerance `0.02 mm`), removing sub-pixel jitter while preserving smooth curved contours and exact sharp corners. + --- ## 🛠️ Building Standalone C/C++ Libraries (Windows & Linux) diff --git a/build.sh b/build.sh index ce82543..f4a2375 100755 --- a/build.sh +++ b/build.sh @@ -36,7 +36,7 @@ NC='\033[0m' # No Color DO_CLEAN=false DO_TEST=false DO_INSTALL=false -DO_WHEEL=false +DO_WHEEL=true BUILD_TYPE="Release" BUILD_C_ONLY=false NUM_JOBS=$(nproc 2>/dev/null || sysctl -n hw.ncpu 2>/dev/null || echo 4) @@ -60,8 +60,13 @@ while [[ $# -gt 0 ]]; do DO_WHEEL=true shift ;; + --no-wheel) + DO_WHEEL=false + shift + ;; --c-only) BUILD_C_ONLY=true + DO_WHEEL=false shift ;; --debug) @@ -84,7 +89,8 @@ while [[ $# -gt 0 ]]; do echo " --clean Remove build directory, libraries, and binaries first" echo " --test Run unit tests and benchmarks upon completion" echo " --install Install Python package into active environment via pip" - echo " --wheel Build standalone Python .whl distribution package in dist/" + echo " --wheel Build standalone Python .whl distribution package in dist/ (default: enabled)" + echo " --no-wheel Skip building Python .whl distribution package" echo " --c-only Build only C/C++ libraries (lib/ and bin/)" echo " --debug Compile with Debug flags" echo " --release Compile with Release flags (default)" @@ -173,9 +179,9 @@ if [ "$BUILD_C_ONLY" = false ]; then echo -e " * Python Backend: ${CYAN}python/cudaslicer/_cudaslicer_backend*.so${NC}" fi -# Optional Python Wheel Package Build +# Python Wheel Package Build if [ "$DO_WHEEL" = true ]; then - echo -e "\n${BLUE}[Optional] Building Python Wheel distribution (.whl)...${NC}" + echo -e "\n${BLUE}[5/5] Building Python Wheel package (.whl)...${NC}" mkdir -p dist python3 setup.py bdist_wheel WHL_FILE=$(ls -t dist/*.whl 2>/dev/null | head -n 1 || true) @@ -193,7 +199,7 @@ fi # Optional Tests & Benchmarks if [ "$DO_TEST" = true ]; then - echo -e "\n${BLUE}[5/5] Running Test Suite & Benchmarks...${NC}" + echo -e "\n${BLUE}[Test Suite] Running Tests & Benchmarks...${NC}" echo -e "\n--- [A] Running Python GPU SDF & CSG Unit Tests ---" python3 tests_and_tools/test_gpu_sdf.py diff --git a/include/gpu_sdf.h b/include/gpu_sdf.h index a6b6717..8d0bd23 100644 --- a/include/gpu_sdf.h +++ b/include/gpu_sdf.h @@ -40,16 +40,7 @@ typedef enum { 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 @@ -180,16 +171,7 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_shell( 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. diff --git a/pyproject.toml b/pyproject.toml index 22b8a92..a82e57b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,10 +4,20 @@ build-backend = "setuptools.build_meta" [project] name = "cudaslicer" -version = "1.1.0" +version = "1.0.0" description = "GPU-Accelerated 2D/3D Slicing, Signed Distance Fields (SDF) & Polygon Offsetting Library" readme = "README.md" requires-python = ">=3.8" +authors = [ + { name = "CUDA Slicer Team" } +] +classifiers = [ + "Programming Language :: Python :: 3", + "Programming Language :: C++", + "Operating System :: POSIX :: Linux", + "Operating System :: Microsoft :: Windows", + "Topic :: Scientific/Engineering", +] dependencies = [ "numpy>=1.20.0" ] diff --git a/python/cudaslicer/__init__.py b/python/cudaslicer/__init__.py index 2e66c68..7e679df 100644 --- a/python/cudaslicer/__init__.py +++ b/python/cudaslicer/__init__.py @@ -9,7 +9,7 @@ with automatic fallback and dynamic backend switching. from typing import List, Dict, Union, Any, Optional, Sequence, Tuple import numpy as np -__version__ = "1.1.0" +__version__ = "1.0.0" # Try importing the compiled C++/CUDA pybind11 module try: @@ -269,7 +269,9 @@ def slice_mesh_and_offset( pad_mm: float = 1.5, return_infill_mask: bool = True, offset_distances: Optional[Sequence[float]] = None, - backend: Optional[str] = "auto" + backend: Optional[str] = "auto", + simplify_tolerance_mm: float = 0.0, + subpixel_lerp: bool = True ) -> List[LayerResult]: """ 3D mesh slicing, rasterization, distance field, and perimeter offsetting (CUDA or CPU). @@ -296,6 +298,10 @@ def slice_mesh_and_offset( Custom offset distances. backend : str, optional, default='auto' Compute backend: 'auto' (CUDA if available, else CPU), 'cuda' (require CUDA), or 'cpu'. + simplify_tolerance_mm : float, default=0.0 + Ramer-Douglas-Peucker (RDP) polygon simplification tolerance in mm (0.0 to disable). + subpixel_lerp : bool, default=True + Continuous distance-field sub-pixel linear interpolation for smooth arbitrary angles. Returns ------- @@ -338,7 +344,9 @@ def slice_mesh_and_offset( int(perimeters_count), float(layer_width), float(pitch_mm), float(pad_mm), bool(return_infill_mask), offsets_float, - backend_str + backend_str, + float(simplify_tolerance_mm), + bool(subpixel_lerp) ) return [LayerResult(r) for r in raw_results] @@ -412,15 +420,7 @@ class 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, diff --git a/python/cudaslicer/bindings.cpp b/python/cudaslicer/bindings.cpp index f315f1e..e922957 100644 --- a/python/cudaslicer/bindings.cpp +++ b/python/cudaslicer/bindings.cpp @@ -252,7 +252,9 @@ py::list py_slice_mesh_and_offset( float pad_mm, bool export_infill_mask, const std::vector& offset_distances, - const std::string& backend + const std::string& backend, + float simplify_tolerance_mm, + bool subpixel_lerp ) { py::buffer_info buf = triangles_arr.request(); size_t num_tris = 0; @@ -275,6 +277,8 @@ py::list py_slice_mesh_and_offset( config.layer_width = layer_width; config.first_layer_height = first_layer_height; 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 = parse_backend_str(backend); @@ -505,24 +509,6 @@ public: 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, @@ -760,6 +746,8 @@ PYBIND11_MODULE(_cudaslicer_backend, m) { py::arg("export_infill_mask") = true, py::arg("offset_distances") = std::vector(), py::arg("backend") = "auto", + py::arg("simplify_tolerance_mm") = 0.0f, + py::arg("subpixel_lerp") = true, "Mesh slicing, perimeter offsetting, and infill mask generation."); // GPU SDF Stack Factory Functions @@ -798,12 +786,6 @@ PYBIND11_MODULE(_cudaslicer_backend, m) { .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, diff --git a/setup.py b/setup.py index 0bc7ead..ec41022 100644 --- a/setup.py +++ b/setup.py @@ -39,24 +39,9 @@ class CMakeBuild(build_ext): setup( name="cudaslicer", - version="1.1.0", - author="CUDA Slicer Team", - description="GPU-Accelerated 2D/3D Slicing & Polygon Offsetting Library", - long_description=open("README.md").read() if os.path.exists("README.md") else "", - long_description_content_type="text/markdown", + version="1.0.0", package_dir={"": "python"}, packages=find_packages(where="python"), ext_modules=[CMakeExtension("cudaslicer._cudaslicer_backend")], cmdclass={"build_ext": CMakeBuild}, - python_requires=">=3.8", - install_requires=[ - "numpy>=1.20.0", - ], - classifiers=[ - "Programming Language :: Python :: 3", - "Programming Language :: C++", - "Operating System :: POSIX :: Linux", - "Operating System :: Microsoft :: Windows", - "Topic :: Scientific/Engineering", - ], ) diff --git a/src/cpu_offset.cpp b/src/cpu_offset.cpp index 6fb39f0..2274db3 100644 --- a/src/cpu_offset.cpp +++ b/src/cpu_offset.cpp @@ -532,7 +532,7 @@ int cpu_offset_slice_mesh_and_offset( auto t_idx_1 = std::chrono::high_resolution_clock::now(); double t_idx_ms = std::chrono::duration(t_idx_1 - t_idx_0).count(); - // 2. Parallel Layer Slicing, Parity Scan, EDT, Thresholding, Vectorization + // 2. Parallel Layer Slicing, Parity Scan, Exact Segment SDF, Thresholding, Vectorization GpuLayerResult* res_array = (GpuLayerResult*)malloc(num_layers * sizeof(GpuLayerResult)); memset(res_array, 0, num_layers * sizeof(GpuLayerResult)); @@ -549,16 +549,18 @@ int cpu_offset_slice_mesh_and_offset( std::vector parity_words((size_t)rows * row_stride_words, 0); std::vector layer_mask(layer_pixels, 0); - std::vector dist_sq(layer_pixels); - std::vector g_col_buf(layer_pixels); - std::vector row_f(cols), row_d(cols); - std::vector v_buf(cols + 1); - std::vector z_buf(cols + 2); - std::vector bit_words((size_t)rows * row_stride_words); + std::vector layer_sdf(layer_pixels); + + struct CpuSegment2D { + float x0, y0, x1, y1; + }; + std::vector layer_segs; + layer_segs.reserve(256); for (size_t l = start_layer; l < end_layer; ++l) { float z0 = z_first + l * layer_height; std::fill(parity_words.begin(), parity_words.end(), 0u); + layer_segs.clear(); // Slicing & Parity Rasterization for Layer l const auto& tri_indices = layer_tri_indices[l]; @@ -588,29 +590,33 @@ int cpu_offset_slice_mesh_and_offset( if (count != 2) continue; - float gx0 = (pts_x[0] - bounds_min_x) / pitch - 0.5f; - float gy0 = (pts_y[0] - bounds_min_y) / pitch - 0.5f; - float gx1 = (pts_x[1] - bounds_min_x) / pitch - 0.5f; - float gy1 = (pts_y[1] - bounds_min_y) / pitch - 0.5f; + layer_segs.push_back({ pts_x[0], pts_y[0], pts_x[1], pts_y[1] }); + + float gx0 = (pts_x[0] - bounds_min_x) / pitch; + float gy0 = (pts_y[0] - bounds_min_y) / pitch; + float gx1 = (pts_x[1] - bounds_min_x) / pitch; + float gy1 = (pts_y[1] - bounds_min_y) / pitch; + + if (gy0 == gy1) continue; if (gy0 > gy1) { std::swap(gx0, gx1); std::swap(gy0, gy1); } - float dy = gy1 - gy0; - if (dy < 1e-6f) continue; - int y_start = (int)std::ceil(gy0); - int y_end = (int)std::floor(gy1 - 1e-6f); + int y_end = (int)std::ceil(gy1) - 1; if (y_start < 0) y_start = 0; if (y_end >= rows) y_end = rows - 1; + float dy = gy1 - gy0; + if (dy < 1e-7f) continue; + for (int y = y_start; y <= y_end; ++y) { float t = ((float)y - gy0) / dy; float gx = gx0 + t * (gx1 - gx0); - int x_int = (int)std::round(gx + 0.5f); + int x_int = (int)std::ceil(gx); if (x_int < 0) x_int = 0; if (x_int >= cols) x_int = cols - 1; @@ -635,43 +641,46 @@ int cpu_offset_slice_mesh_and_offset( } } - // Compute exact 2D EDT - edt_2d_layer(layer_mask.data(), cols, rows, dist_sq.data(), - g_col_buf, row_f, row_d, v_buf, z_buf); + // Compute exact continuous signed distance field from analytical layer segments + int num_segs = (int)layer_segs.size(); + for (int y = 0; y < rows; ++y) { + float py = bounds_min_y + (float)y * pitch; + for (int x = 0; x < cols; ++x) { + float px = bounds_min_x + (float)x * pitch; + float min_d2 = 1e20f; + for (int s = 0; s < num_segs; ++s) { + float sx0 = layer_segs[s].x0; + float sy0 = layer_segs[s].y0; + float sx1 = layer_segs[s].x1; + float sy1 = layer_segs[s].y1; + float dx = sx1 - sx0; + float dy = sy1 - sy0; + float l2 = dx * dx + dy * dy; + float t = (l2 > 1e-12f) ? ((px - sx0) * dx + (py - sy0) * dy) / l2 : 0.0f; + t = std::max(0.0f, std::min(1.0f, t)); + float qx = sx0 + t * dx; + float qy = sy0 + t * dy; + float d2 = (px - qx) * (px - qx) + (py - qy) * (py - qy); + if (d2 < min_d2) min_d2 = d2; + } + float d = sqrtf(min_d2); + uint8_t m = layer_mask[(size_t)y * cols + x]; + layer_sdf[(size_t)y * cols + x] = (m > 0) ? d : -d; + } + } res_array[l].z_height = z0; res_array[l].num_perimeters = P; res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter)); - // Perimeters + // Perimeters with continuous SDF and sub-pixel LERP for (int o = 0; o < P; ++o) { - float thresh = resolved_offsets[o]; - float thresh_sq_px = (thresh / pitch) * (thresh / pitch); - - for (int y = 0; y < rows; ++y) { - const uint8_t* row_m = layer_mask.data() + (size_t)y * cols; - const float* row_d2 = dist_sq.data() + (size_t)y * cols; - unsigned int* row_w = bit_words.data() + (size_t)y * row_stride_words; - - for (int w = 0; w < row_stride_words; ++w) { - unsigned int packed = 0; - int x_base = w * 32; - for (int b = 0; b < 32; ++b) { - int x = x_base + b; - if (x < cols) { - if (row_m[x] > 0 && row_d2[x] >= thresh_sq_px) { - packed |= (1u << b); - } - } - } - row_w[w] = packed; - } - } + float iso_thresh = resolved_offsets[o]; std::vector contours; - vectorize_bitpack_layer(bit_words.data(), cols, rows, row_stride_words, - bounds_min_x, bounds_min_y, pitch, - next_edge_buf, active_edges_buf, visited_buf, contours); + vectorize_sdf_layer(layer_sdf.data(), cols, rows, bounds_min_x, bounds_min_y, pitch, + iso_thresh, config->simplify_tolerance_mm, config->subpixel_lerp, + next_edge_buf, active_edges_buf, visited_buf, contours); res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o]; res_array[l].perimeters[o].num_contours = contours.size(); @@ -684,35 +693,13 @@ int cpu_offset_slice_mesh_and_offset( } // Infill pass - float infill_thresh = resolved_offsets[P]; - float infill_thresh_sq_px = (infill_thresh / pitch) * (infill_thresh / pitch); - - for (int y = 0; y < rows; ++y) { - const uint8_t* row_m = layer_mask.data() + (size_t)y * cols; - const float* row_d2 = dist_sq.data() + (size_t)y * cols; - unsigned int* row_w = bit_words.data() + (size_t)y * row_stride_words; - - for (int w = 0; w < row_stride_words; ++w) { - unsigned int packed = 0; - int x_base = w * 32; - for (int b = 0; b < 32; ++b) { - int x = x_base + b; - if (x < cols) { - if (row_m[x] > 0 && row_d2[x] >= infill_thresh_sq_px) { - packed |= (1u << b); - } - } - } - row_w[w] = packed; - } - } - + float infill_iso_thresh = resolved_offsets[P]; std::vector infill_contours; - vectorize_bitpack_layer(bit_words.data(), cols, rows, row_stride_words, - bounds_min_x, bounds_min_y, pitch, - next_edge_buf, active_edges_buf, visited_buf, infill_contours); + vectorize_sdf_layer(layer_sdf.data(), cols, rows, bounds_min_x, bounds_min_y, pitch, + infill_iso_thresh, config->simplify_tolerance_mm, config->subpixel_lerp, + next_edge_buf, active_edges_buf, visited_buf, infill_contours); - res_array[l].infill.offset_distance_mm = resolved_offsets[P]; + 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(); @@ -724,17 +711,9 @@ int cpu_offset_slice_mesh_and_offset( } if (config->export_infill_mask) { - uint8_t* infill_mask = (uint8_t*)malloc((size_t)cols * rows); - for (int y = 0; y < rows; ++y) { - const unsigned int* row_words = bit_words.data() + (size_t)y * row_stride_words; - uint8_t* row_out = infill_mask + (size_t)y * cols; - int x_out = 0; - for (int w = 0; w < row_stride_words && x_out < cols; ++w) { - unsigned int word = row_words[w]; - for (int b = 0; b < 32 && x_out < cols; ++b, ++x_out) { - row_out[x_out] = (uint8_t)((word >> b) & 1u); - } - } + uint8_t* infill_mask = (uint8_t*)malloc(layer_pixels); + for (size_t i = 0; i < layer_pixels; ++i) { + infill_mask[i] = (layer_sdf[i] >= infill_iso_thresh) ? 1 : 0; } res_array[l].infill.mask = infill_mask; } else { diff --git a/src/gpu_offset.cu b/src/gpu_offset.cu index fb09087..1324e12 100644 --- a/src/gpu_offset.cu +++ b/src/gpu_offset.cu @@ -1,4 +1,5 @@ #include "gpu_offset.h" +#include "gpu_sdf.h" #include "cpu_offset.h" #include "vectorizer.h" #include @@ -763,295 +764,25 @@ GPU_OFFSET_API int gpu_offset_slice_mesh_and_offset( if (gpu_offset_init(0) != GPU_OFFSET_SUCCESS) return GPU_OFFSET_ERR_GENERIC; } - auto t_total_start = std::chrono::high_resolution_clock::now(); - float pad = (config->pad_mm > 0.0f) ? config->pad_mm : 1.5f; float pitch = (config->pitch_mm > 0.0f) ? config->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; if (h_first <= 0.0f) h_first = config->first_layer_height; if (h_first <= 0.0f) h_first = layer_height; - 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) + 1; + GpuSdfStackHandle stack = gpu_sdf_from_mesh( + triangles, num_triangles, + min_x, min_y, min_z, max_x, max_y, max_z, + layer_height, h_first, + pitch, pad + ); + if (!stack) { + return GPU_OFFSET_ERR_GENERIC; } - if (num_layers <= 0) return GPU_OFFSET_ERR_INVALID_ARG; - if (out_num_layers) *out_num_layers = num_layers; - - 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(); - - int row_stride_words = (cols + 31) / 32; - size_t total_words = (size_t)num_layers * rows * row_stride_words; - size_t layer_pixels = (size_t)cols * rows; - size_t total_volume = (size_t)num_layers * layer_pixels; - size_t total_bit_words = num_offsets * total_words; - size_t total_bit_bytes = total_bit_words * sizeof(unsigned int); - - cudaEvent_t ev_start, ev_h2d_done, ev_slice_done, ev_jfa_done, ev_thresh_done, ev_d2h_done; - CUDA_CHECK(cudaEventCreate(&ev_start)); - CUDA_CHECK(cudaEventCreate(&ev_h2d_done)); - CUDA_CHECK(cudaEventCreate(&ev_slice_done)); - CUDA_CHECK(cudaEventCreate(&ev_jfa_done)); - CUDA_CHECK(cudaEventCreate(&ev_thresh_done)); - CUDA_CHECK(cudaEventCreate(&ev_d2h_done)); - - GpuTriangle3D* d_triangles; - unsigned int* d_parity; - uint8_t* d_masks; - int* d_seeds_a; - int* d_seeds_b; - unsigned int* d_out_bit_words; - float* d_thresholds; - - size_t tri_bytes = num_triangles * sizeof(GpuTriangle3D); - CUDA_CHECK(cudaMalloc(&d_triangles, tri_bytes)); - CUDA_CHECK(cudaMalloc(&d_parity, total_words * sizeof(unsigned int))); - CUDA_CHECK(cudaMalloc(&d_masks, total_volume * sizeof(uint8_t))); - CUDA_CHECK(cudaMalloc(&d_seeds_a, total_volume * sizeof(int))); - CUDA_CHECK(cudaMalloc(&d_seeds_b, total_volume * sizeof(int))); - CUDA_CHECK(cudaMalloc(&d_out_bit_words, total_bit_bytes)); - CUDA_CHECK(cudaMalloc(&d_thresholds, num_offsets * sizeof(float))); - - unsigned int* h_pinned_bit_words = nullptr; - CUDA_CHECK(cudaMallocHost(&h_pinned_bit_words, total_bit_bytes)); - - // 1. Upload Triangles & Thresholds (H2D) - CUDA_CHECK(cudaEventRecord(ev_start, g_cuda_stream)); - CUDA_CHECK(cudaMemcpyAsync(d_triangles, triangles, tri_bytes, cudaMemcpyHostToDevice, g_cuda_stream)); - CUDA_CHECK(cudaMemcpyAsync(d_thresholds, resolved_offsets.data(), - num_offsets * sizeof(float), cudaMemcpyHostToDevice, g_cuda_stream)); - CUDA_CHECK(cudaEventRecord(ev_h2d_done, g_cuda_stream)); - - // 2. GPU Direct Mesh Slicing & Scanline Rasterization - CUDA_CHECK(cudaMemsetAsync(d_parity, 0, total_words * sizeof(unsigned int), g_cuda_stream)); - - int block_sz = 256; - int grid_sz = (num_triangles + block_sz - 1) / block_sz; - gpu_mesh_slice_raster_kernel<<>>( - d_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 - ); - CUDA_CHECK(cudaEventRecord(ev_slice_done, g_cuda_stream)); - - // 3. Batched 3D GPU JFA Distance Transform - dim3 jfa_block(16, 16, 1); - dim3 jfa_grid((cols + 15) / 16, (rows + 15) / 16, num_layers); - - batched_jfa_init_kernel<<>>( - d_masks, d_seeds_a, cols, rows, num_layers - ); - - float max_offset_d = 0.0f; - for (size_t o = 0; o < num_offsets; ++o) { - if (resolved_offsets[o] > max_offset_d) { - max_offset_d = resolved_offsets[o]; - } - } - int max_radius_px = (int)std::ceil(max_offset_d / pitch) + 2; - - int step = 1; - while (step * 2 <= max_radius_px && step < 512) step *= 2; - if (step < 8) step = 8; - - int* in_seeds = d_seeds_a; - int* out_seeds = d_seeds_b; - - while (step >= 1) { - batched_jfa_step_kernel<<>>( - in_seeds, out_seeds, cols, rows, num_layers, step - ); - std::swap(in_seeds, out_seeds); - step /= 2; - } - - batched_jfa_step_kernel<<>>( - in_seeds, out_seeds, cols, rows, num_layers, 2 - ); - std::swap(in_seeds, out_seeds); - - batched_jfa_step_kernel<<>>( - in_seeds, out_seeds, cols, rows, num_layers, 1 - ); - std::swap(in_seeds, out_seeds); - - CUDA_CHECK(cudaEventRecord(ev_jfa_done, g_cuda_stream)); - - // 4. GPU Fused Multi-Perimeter Offset Thresholding & 32-bit Bit-Packing - dim3 thresh_block(32, 1, 1); - dim3 thresh_grid((row_stride_words + 31) / 32, rows, num_layers); - - batched_jfa_threshold_bitpack_kernel<<>>( - in_seeds, d_masks, d_out_bit_words, - cols, rows, num_layers, - pitch, d_thresholds, (int)num_offsets, row_stride_words - ); - CUDA_CHECK(cudaEventRecord(ev_thresh_done, g_cuda_stream)); - - // 5. Ultra-Fast D2H Transfer using Pinned Host Memory DMA - CUDA_CHECK(cudaMemcpyAsync(h_pinned_bit_words, d_out_bit_words, total_bit_bytes, - cudaMemcpyDeviceToHost, g_cuda_stream)); - CUDA_CHECK(cudaEventRecord(ev_d2h_done, g_cuda_stream)); - CUDA_CHECK(cudaStreamSynchronize(g_cuda_stream)); - - // 6. Multi-Threaded C++ Bit-Packed 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 = z_first + l * layer_height; - res_array[l].num_perimeters = P; - res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter)); - - for (int o = 0; o < P; ++o) { - size_t offset_layer_word_idx = (o * num_layers + l) * ((size_t)rows * row_stride_words); - const unsigned int* bit_ptr = h_pinned_bit_words + offset_layer_word_idx; - - std::vector contours; - vectorize_bitpack_layer(bit_ptr, cols, rows, row_stride_words, bounds_min_x, bounds_min_y, pitch, - next_edge_buf, active_edges_buf, visited_buf, contours); - - res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o]; - 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 - size_t infill_offset_layer_word_idx = (P * num_layers + l) * ((size_t)rows * row_stride_words); - const unsigned int* infill_bit_ptr = h_pinned_bit_words + infill_offset_layer_word_idx; - - std::vector infill_contours; - vectorize_bitpack_layer(infill_bit_ptr, cols, rows, row_stride_words, bounds_min_x, bounds_min_y, pitch, - next_edge_buf, active_edges_buf, visited_buf, infill_contours); - - res_array[l].infill.offset_distance_mm = resolved_offsets[P]; - 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 y = 0; y < rows; ++y) { - const unsigned int* row_words = infill_bit_ptr + (size_t)y * row_stride_words; - uint8_t* row_out = infill_mask + (size_t)y * cols; - int x_out = 0; - for (int w = 0; w < row_stride_words && x_out < cols; ++w) { - unsigned int word = row_words[w]; - for (int b = 0; b < 32 && x_out < cols; ++b, ++x_out) { - row_out[x_out] = (uint8_t)((word >> b) & 1u); - } - } - } - 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, (size_t)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) { - float h2d_ms = 0, slice_ms = 0, jfa_ms = 0, thresh_ms = 0, d2h_ms = 0; - cudaEventElapsedTime(&h2d_ms, ev_start, ev_h2d_done); - cudaEventElapsedTime(&slice_ms, ev_h2d_done, ev_slice_done); - cudaEventElapsedTime(&jfa_ms, ev_slice_done, ev_jfa_done); - cudaEventElapsedTime(&thresh_ms, ev_jfa_done, ev_thresh_done); - cudaEventElapsedTime(&d2h_ms, ev_thresh_done, ev_d2h_done); - - stats->t_h2d_ms = h2d_ms; - stats->t_mesh_slice_ms = slice_ms; - stats->t_raster_ms = slice_ms; - stats->t_gpu_jfa_ms = jfa_ms; - stats->t_gpu_thresh_ms = thresh_ms; - stats->t_d2h_ms = d2h_ms; - stats->t_vectorize_ms = t_vec_ms; - stats->t_total_pipeline_ms = t_total_ms; - } - - // Free CUDA & Pinned resources - cudaFree(d_triangles); - cudaFree(d_parity); - cudaFree(d_masks); - cudaFree(d_seeds_a); - cudaFree(d_seeds_b); - cudaFree(d_out_bit_words); - cudaFree(d_thresholds); - cudaFreeHost(h_pinned_bit_words); - - cudaEventDestroy(ev_start); - cudaEventDestroy(ev_h2d_done); - cudaEventDestroy(ev_slice_done); - cudaEventDestroy(ev_jfa_done); - cudaEventDestroy(ev_thresh_done); - cudaEventDestroy(ev_d2h_done); - - *out_results = res_array; - return GPU_OFFSET_SUCCESS; + int ret = gpu_sdf_vectorize(stack, config, out_results, out_num_layers, stats); + gpu_sdf_destroy(stack); + return ret; } else { return cpu_offset_slice_mesh_and_offset(triangles, num_triangles, min_x, min_y, min_z, max_x, max_y, max_z, layer_height, first_layer_height, config, out_results, out_num_layers, stats); diff --git a/src/gpu_sdf.cu b/src/gpu_sdf.cu index 020df77..b507893 100644 --- a/src/gpu_sdf.cu +++ b/src/gpu_sdf.cu @@ -345,61 +345,7 @@ __global__ void gpu_sdf_shell_kernel( 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 @@ -597,16 +543,210 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_masks( return stack; } -// Forward declare GPU mesh rasterizer helper from gpu_offset.cu -extern int gpu_slice_mesh_to_device_mask( - const GpuTriangle3D* triangles, +struct GpuSegment2D { + float x0, y0; + float x1, y1; +}; + +__global__ void gpu_mesh_slice_raster_and_segments_kernel( + const GpuTriangle3D* __restrict__ 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 -); + float z_min_base, + float first_layer_height, + float layer_height, + int num_layers, + float bounds_min_x, + float bounds_min_y, + float pitch_mm, + int cols, + int rows, + int row_stride_words, + unsigned int* __restrict__ parity_words, + GpuSegment2D* __restrict__ layer_segments, + int* __restrict__ layer_seg_counts, + int max_segs_per_layer +) { + size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x; + if (i >= num_triangles) return; + + GpuTriangle3D tri = triangles[i]; + + float tz_min = fminf(tri.v0[2], fminf(tri.v1[2], tri.v2[2])); + float tz_max = fmaxf(tri.v0[2], fmaxf(tri.v1[2], tri.v2[2])); + + float z_first = z_min_base + first_layer_height; + int l_start = (tz_min <= z_first) ? 0 : (int)ceilf((tz_min - z_first) / layer_height); + int l_end = (int)floorf((tz_max - z_first) / layer_height); + + if (l_start < 0) l_start = 0; + if (l_end >= num_layers) l_end = num_layers - 1; + if (l_start > l_end) return; + + const float* v[3] = { tri.v0, tri.v1, tri.v2 }; + + for (int l = l_start; l <= l_end; ++l) { + float z0 = z_first + l * layer_height; + if (z0 < tz_min || z0 > tz_max) continue; + + float pts_x[3]; + float pts_y[3]; + int count = 0; + + #pragma unroll + for (int e = 0; e < 3; ++e) { + int next = (e + 1) % 3; + float za = v[e][2]; + float zb = v[next][2]; + + if ((za <= z0 && zb >= z0) || (zb <= z0 && za >= z0)) { + if (fabsf(zb - za) > 1e-7f) { + float t = (z0 - za) / (zb - za); + pts_x[count] = v[e][0] + t * (v[next][0] - v[e][0]); + pts_y[count] = v[e][1] + t * (v[next][1] - v[e][1]); + count++; + if (count == 2) break; + } + } + } + + if (count != 2) continue; + + if (layer_segments && layer_seg_counts) { + int seg_idx = atomicAdd(&layer_seg_counts[l], 1); + if (seg_idx < max_segs_per_layer) { + GpuSegment2D seg; + seg.x0 = pts_x[0]; + seg.y0 = pts_y[0]; + seg.x1 = pts_x[1]; + seg.y1 = pts_y[1]; + layer_segments[(size_t)l * max_segs_per_layer + seg_idx] = seg; + } + } + + float gx0 = (pts_x[0] - bounds_min_x) / pitch_mm; + float gy0 = (pts_y[0] - bounds_min_y) / pitch_mm; + float gx1 = (pts_x[1] - bounds_min_x) / pitch_mm; + float gy1 = (pts_y[1] - bounds_min_y) / pitch_mm; + + if (gy0 == gy1) continue; + + if (gy0 > gy1) { + float tx = gx0; gx0 = gx1; gx1 = tx; + float ty = gy0; gy0 = gy1; gy1 = ty; + } + + int y_start = (int)ceilf(gy0); + int y_end = (int)ceilf(gy1) - 1; + + if (y_start < 0) y_start = 0; + if (y_end >= rows) y_end = rows - 1; + + float dy = gy1 - gy0; + if (dy < 1e-7f) continue; + + size_t layer_word_offset = (size_t)l * rows * row_stride_words; + + for (int y = y_start; y <= y_end; ++y) { + float t = ((float)y - gy0) / dy; + float gx = gx0 + t * (gx1 - gx0); + int x_int = (int)ceilf(gx); + if (x_int < 0) x_int = 0; + if (x_int >= cols) x_int = cols - 1; + + int w = x_int / 32; + int b = x_int % 32; + size_t w_idx = layer_word_offset + (size_t)y * row_stride_words + w; + atomicXor(&parity_words[w_idx], (1u << b)); + } + } +} + +__global__ void gpu_parity_scan_to_mask_sdf_kernel( + const unsigned int* __restrict__ parity_words, + uint8_t* __restrict__ masks, + int cols, + int rows, + int num_layers, + int row_stride_words +) { + int y = blockIdx.x; + int l = blockIdx.y; + if (y >= rows || l >= num_layers) return; + + size_t layer_word_offset = (size_t)l * rows * row_stride_words; + size_t layer_pixel_offset = (size_t)l * rows * cols; + size_t row_pixel_offset = layer_pixel_offset + (size_t)y * cols; + + uint8_t inside = 0; + int x_out = 0; + + for (int w = 0; w < row_stride_words; ++w) { + unsigned int word = parity_words[layer_word_offset + (size_t)y * row_stride_words + w]; + for (int b = 0; b < 32 && x_out < cols; ++b, ++x_out) { + inside ^= (uint8_t)((word >> b) & 1u); + masks[row_pixel_offset + x_out] = inside; + } + } +} + +__global__ void gpu_compute_exact_mesh_sdf_kernel( + const uint8_t* __restrict__ masks, + const GpuSegment2D* __restrict__ layer_segments, + const int* __restrict__ layer_seg_counts, + int max_segs_per_layer, + float* __restrict__ out_sdf, + int cols, + int rows, + int num_layers, + float bounds_min_x, + float bounds_min_y, + float pitch_mm +) { + int x = blockIdx.x * blockDim.x + threadIdx.x; + int y = blockIdx.y * blockDim.y + threadIdx.y; + int l = blockIdx.z; + + if (x >= cols || y >= rows || l >= num_layers) return; + + size_t layer_pixel_offset = (size_t)l * rows * cols; + size_t idx = layer_pixel_offset + (size_t)y * cols + x; + + int seg_count = layer_seg_counts[l]; + if (seg_count > max_segs_per_layer) seg_count = max_segs_per_layer; + + if (seg_count == 0) { + out_sdf[idx] = -1e5f; + return; + } + + float px = bounds_min_x + (float)x * pitch_mm; + float py = bounds_min_y + (float)y * pitch_mm; + + float min_d2 = 1e10f; + const GpuSegment2D* segs = layer_segments + (size_t)l * max_segs_per_layer; + + for (int s = 0; s < seg_count; ++s) { + GpuSegment2D seg = segs[s]; + float vx = seg.x1 - seg.x0; + float vy = seg.y1 - seg.y0; + float wx = px - seg.x0; + float wy = py - seg.y0; + float len_sq = vx * vx + vy * vy; + float t = (len_sq > 1e-10f) ? (wx * vx + wy * vy) / len_sq : 0.0f; + if (t < 0.0f) t = 0.0f; + else if (t > 1.0f) t = 1.0f; + float qx = seg.x0 + t * vx; + float qy = seg.y0 + t * vy; + float dx = px - qx; + float dy = py - qy; + float d2 = dx * dx + dy * dy; + if (d2 < min_d2) min_d2 = d2; + } + + float d = sqrtf(min_d2); + uint8_t m = masks[idx]; + out_sdf[idx] = (m > 0) ? d : -d; +} GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_mesh( const GpuTriangle3D* triangles, @@ -656,107 +796,67 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_mesh( GpuSdfStackHandle stack = gpu_sdf_create(&meta, z_heights_um.data()); if (!stack) return nullptr; + int row_stride_words = (cols + 31) / 32; + size_t total_words = (size_t)num_layers * rows * row_stride_words; 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); + + int max_segs_per_layer = std::min(8192, (int)num_triangles); + if (max_segs_per_layer < 256) max_segs_per_layer = 256; + size_t total_segs = (size_t)num_layers * max_segs_per_layer; GpuTriangle3D* d_triangles = nullptr; + unsigned int* d_parity = 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; + GpuSegment2D* d_layer_segments = nullptr; + int* d_layer_seg_counts = nullptr; if (cudaMalloc(&d_triangles, tri_bytes) != cudaSuccess || + cudaMalloc(&d_parity, total_words * sizeof(unsigned int)) != 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) { + cudaMalloc(&d_layer_segments, total_segs * sizeof(GpuSegment2D)) != cudaSuccess || + cudaMalloc(&d_layer_seg_counts, num_layers * sizeof(int)) != cudaSuccess) { if (d_triangles) cudaFree(d_triangles); + if (d_parity) cudaFree(d_parity); 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); + if (d_layer_segments) cudaFree(d_layer_segments); + if (d_layer_seg_counts) cudaFree(d_layer_seg_counts); gpu_sdf_destroy(stack); return nullptr; } + cudaMemsetAsync(d_parity, 0, total_words * sizeof(unsigned int), stack->stream); + cudaMemsetAsync(d_layer_seg_counts, 0, num_layers * sizeof(int), stack->stream); cudaMemcpyAsync(d_triangles, triangles, tri_bytes, cudaMemcpyHostToDevice, stack->stream); - int slice_ret = gpu_slice_mesh_to_device_mask( + int block_sz = 256; + int grid_sz = (num_triangles + block_sz - 1) / block_sz; + gpu_mesh_slice_raster_and_segments_kernel<<stream>>>( 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 + bounds_min_x, bounds_min_y, pitch, cols, rows, row_stride_words, + d_parity, d_layer_segments, d_layer_seg_counts, max_segs_per_layer ); - 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 + dim3 scan_grid(rows, num_layers); + gpu_parity_scan_to_mask_sdf_kernel<<stream>>>( + d_parity, d_masks, cols, rows, num_layers, row_stride_words ); - 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 + dim3 exact_block(16, 16, 1); + dim3 exact_grid((cols + 15) / 16, (rows + 15) / 16, num_layers); + gpu_compute_exact_mesh_sdf_kernel<<stream>>>( + d_masks, d_layer_segments, d_layer_seg_counts, max_segs_per_layer, + stack->d_sdf, cols, rows, num_layers, + bounds_min_x, bounds_min_y, pitch ); cudaStreamSynchronize(stack->stream); cudaFree(d_triangles); + cudaFree(d_parity); cudaFree(d_masks); - cudaFree(d_seeds_ext_a); - cudaFree(d_seeds_ext_b); - cudaFree(d_seeds_int_a); - cudaFree(d_seeds_int_b); + cudaFree(d_layer_segments); + cudaFree(d_layer_seg_counts); return stack; } @@ -1002,35 +1102,6 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_shell( 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, @@ -1068,9 +1139,7 @@ GPU_SDF_API int gpu_sdf_vectorize( 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 - } + int subpixel_lerp = config->subpixel_lerp; size_t total_voxels = (size_t)num_layers * rows * cols; size_t sdf_bytes = total_voxels * sizeof(float); @@ -1113,7 +1182,7 @@ GPU_SDF_API int gpu_sdf_vectorize( std::vector contours; vectorize_sdf_layer(layer_sdf_ptr, cols, rows, bounds_min_x, bounds_min_y, pitch_mm, - iso_thresh, simplify_tol, + iso_thresh, simplify_tol, subpixel_lerp, next_edge_buf, active_edges_buf, visited_buf, contours); res_array[l].perimeters[o].offset_distance_mm = iso_thresh; @@ -1126,11 +1195,11 @@ GPU_SDF_API int gpu_sdf_vectorize( } } - // Infill + // Infill Contours (Innermost offset boundary) 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, + infill_iso_thresh, simplify_tol, subpixel_lerp, next_edge_buf, active_edges_buf, visited_buf, infill_contours); res_array[l].infill.offset_distance_mm = infill_iso_thresh; diff --git a/src/vectorizer.h b/src/vectorizer.h index 99b9cf6..d579af4 100644 --- a/src/vectorizer.h +++ b/src/vectorizer.h @@ -143,12 +143,16 @@ static inline std::vector simplify_contour_rdp( } } - std::vector keep(n, false); + std::vector closed_pts = loop_pts; + closed_pts.push_back(loop_pts[0]); + + std::vector keep(n + 1, false); keep[0] = true; keep[split_idx] = true; + keep[n] = true; - rdp_simplify_recursive(loop_pts, 0, split_idx, eps_sq, keep); - rdp_simplify_recursive(loop_pts, split_idx, n - 1, eps_sq, keep); + rdp_simplify_recursive(closed_pts, 0, split_idx, eps_sq, keep); + rdp_simplify_recursive(closed_pts, split_idx, n, eps_sq, keep); std::vector simplified; simplified.reserve(n); @@ -164,6 +168,7 @@ static inline void vectorize_sdf_layer( float min_x, float min_y, float pitch_mm, float iso_val, float simplify_tolerance_mm, + int subpixel_lerp, std::vector& next_edge_buf, std::vector& active_edges_buf, std::vector& visited_buf, @@ -253,7 +258,12 @@ static inline void vectorize_sdf_layer( 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); + GpuPoint2D pt; + if (subpixel_lerp) { + pt = get_sdf_edge_point(cur, sdf, width, height, min_x, min_y, pitch_mm, iso_val); + } else { + pt = get_edge_point(cur, stride_w, min_x, min_y, pitch_mm); + } size_t n = loop_pts.size(); if (n >= 2) { @@ -274,11 +284,13 @@ static inline void vectorize_sdf_layer( cur = next_edge_buf[cur]; if (cur == start_e) break; - if (loop_pts.size() > 50000) break; + if (loop_pts.size() > 100000) break; } if (loop_pts.size() >= 3) { - std::vector simplified = simplify_contour_rdp(loop_pts, simplify_tolerance_mm); + std::vector simplified = (simplify_tolerance_mm > 1e-6f) + ? simplify_contour_rdp(loop_pts, simplify_tolerance_mm) + : loop_pts; if (simplified.size() >= 3) { GpuContour c; c.num_points = simplified.size();