commit 1019a91e34952c8bdc269f17d8bde4ce421617e8 Author: Nicolas Predella Date: Fri Sep 25 18:09:19 2026 +0300 new diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..3bc6c43 --- /dev/null +++ b/.gitignore @@ -0,0 +1,34 @@ +# Build and artifact directories +build/ +bin/ +lib/ +dist/ +*.egg-info/ +__pycache__/ +*.pyc +*.pyo +*.pyd + +# Compiled shared and static libraries / objects +*.so +*.a +*.o +*.obj +*.dylib +*.dll +*.lib +*.exp + +# Testing temp tools, AppImages, and experimental benchmarks +tests_and_tools/ + +# Slicing output artifacts +*.gcode +*.png + +# IDE and OS metadata +.vscode/ +.idea/ +.DS_Store +*.swp +*~ diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..c62b8b2 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,109 @@ +cmake_minimum_required(VERSION 3.18) +project(cudaslicer LANGUAGES C CXX CUDA) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CUDA_STANDARD 14) +set(CMAKE_CUDA_STANDARD_REQUIRED ON) + +# Cross-platform compiler settings +if(MSVC) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /O2 /W3 /MP") + set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler=/MD,/W3") +else() + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -pthread -fPIC") + set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler -fPIC,-Wall") +endif() + +# Set broad CUDA GPU architecture support (Turing, Ampere, Ada Lovelace, Hopper/Blackwell) +if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES) + set(CMAKE_CUDA_ARCHITECTURES "75;80;86;89") +endif() + +include_directories( + ${CMAKE_CURRENT_SOURCE_DIR}/include + ${CMAKE_CURRENT_SOURCE_DIR}/src + ${CUDA_INCLUDE_DIRS} +) + +# ----------------------------------------------------------------------------- +# 1. Standalone C/C++ Shared Library (libcudaslicer_offset.so / cudaslicer_offset.dll) +# ----------------------------------------------------------------------------- +add_library(cudaslicer_offset SHARED + src/gpu_offset.cu + src/cpu_offset.cpp +) +set_target_properties(cudaslicer_offset PROPERTIES + OUTPUT_NAME "cudaslicer_offset" + CUDA_SEPARABLE_COMPILATION OFF +) +target_compile_definitions(cudaslicer_offset PRIVATE GPU_OFFSET_EXPORTS) +target_link_libraries(cudaslicer_offset PUBLIC cudart) + +# ----------------------------------------------------------------------------- +# 2. Standalone C/C++ Static Library +# ----------------------------------------------------------------------------- +add_library(cudaslicer_offset_static STATIC + src/gpu_offset.cu + src/cpu_offset.cpp +) +set_target_properties(cudaslicer_offset_static PROPERTIES + OUTPUT_NAME "cudaslicer_offset" + CUDA_SEPARABLE_COMPILATION OFF +) +target_link_libraries(cudaslicer_offset_static PUBLIC cudart) + +# ----------------------------------------------------------------------------- +# 3. C++ Benchmark Binary +# ----------------------------------------------------------------------------- +add_executable(gpu_slicer_c_bench + src/main_bench.cpp + src/stl_slicer.cpp +) +target_link_libraries(gpu_slicer_c_bench PRIVATE cudaslicer_offset cudart) +if(UNIX) + target_link_libraries(gpu_slicer_c_bench PRIVATE pthread) +endif() + +# ----------------------------------------------------------------------------- +# 4. Pybind11 Python Native Extension (_cudaslicer_backend) +# ----------------------------------------------------------------------------- +if(NOT DEFINED Python3_EXECUTABLE) + find_program(Python3_EXECUTABLE NAMES python3 python) +endif() +find_package(Python3 COMPONENTS Interpreter REQUIRED) +execute_process( + COMMAND "${Python3_EXECUTABLE}" -c "import sysconfig; print(sysconfig.get_path('include'))" + OUTPUT_VARIABLE _python3_include_dir + OUTPUT_STRIP_TRAILING_WHITESPACE +) +execute_process( + COMMAND "${Python3_EXECUTABLE}" -c "import pybind11; print(pybind11.get_include())" + OUTPUT_VARIABLE _pybind11_include_dir + OUTPUT_STRIP_TRAILING_WHITESPACE +) +execute_process( + COMMAND "${Python3_EXECUTABLE}" -c "import sysconfig; print(sysconfig.get_config_var('EXT_SUFFIX') or '.so')" + OUTPUT_VARIABLE _ext_suffix + OUTPUT_STRIP_TRAILING_WHITESPACE +) + +add_library(_cudaslicer_backend MODULE + python/cudaslicer/bindings.cpp +) +set_target_properties(_cudaslicer_backend PROPERTIES + PREFIX "" + SUFFIX "${_ext_suffix}" +) +target_include_directories(_cudaslicer_backend PRIVATE "${_python3_include_dir}" "${_pybind11_include_dir}") +target_link_libraries(_cudaslicer_backend PRIVATE cudaslicer_offset_static cudart) + +# ----------------------------------------------------------------------------- +# Output Directories +# ----------------------------------------------------------------------------- +if(NOT CMAKE_LIBRARY_OUTPUT_DIRECTORY) + set_target_properties(cudaslicer_offset PROPERTIES LIBRARY_OUTPUT_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/lib") + set_target_properties(cudaslicer_offset_static PROPERTIES ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/lib") + set_target_properties(gpu_slicer_c_bench PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/bin") + set_target_properties(_cudaslicer_backend PROPERTIES LIBRARY_OUTPUT_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/python/cudaslicer") +endif() diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..5cec1ee --- /dev/null +++ b/Makefile @@ -0,0 +1,45 @@ +NVCC ?= nvcc +CXX ?= g++ +ARCH ?= sm_86 + +CXXFLAGS = -O3 -std=c++17 -Wall -pthread -fPIC -Iinclude -Isrc +NVCCFLAGS = -O3 -std=c++14 -arch=$(ARCH) --use_fast_math -Xcompiler -fPIC,-Wall -Iinclude -Isrc + +LIB_DIR = lib +BIN_DIR = bin +BUILD_DIR = build + +SHARED_LIB = $(LIB_DIR)/libcudaslicer_offset.so +STATIC_LIB = $(LIB_DIR)/libcudaslicer_offset.a +BENCH_BIN = $(BIN_DIR)/gpu_slicer_c_bench + +all: $(SHARED_LIB) $(STATIC_LIB) $(BENCH_BIN) + +$(BUILD_DIR)/gpu_offset.o: src/gpu_offset.cu include/gpu_offset.h src/cpu_offset.h src/vectorizer.h | $(BUILD_DIR) + $(NVCC) $(NVCCFLAGS) -c $< -o $@ + +$(BUILD_DIR)/cpu_offset.o: src/cpu_offset.cpp src/cpu_offset.h src/vectorizer.h include/gpu_offset.h | $(BUILD_DIR) + $(CXX) $(CXXFLAGS) -c $< -o $@ + +$(BUILD_DIR)/stl_slicer.o: src/stl_slicer.cpp src/stl_slicer.h include/gpu_offset.h | $(BUILD_DIR) + $(CXX) $(CXXFLAGS) -c $< -o $@ + +$(BUILD_DIR)/main_bench.o: src/main_bench.cpp src/stl_slicer.h include/gpu_offset.h | $(BUILD_DIR) + $(CXX) $(CXXFLAGS) -c $< -o $@ + +$(SHARED_LIB): $(BUILD_DIR)/gpu_offset.o $(BUILD_DIR)/cpu_offset.o | $(LIB_DIR) + $(NVCC) -shared -arch=$(ARCH) -o $@ $^ -lcudart + +$(STATIC_LIB): $(BUILD_DIR)/gpu_offset.o $(BUILD_DIR)/cpu_offset.o | $(LIB_DIR) + ar rcs $@ $^ + +$(BENCH_BIN): $(BUILD_DIR)/main_bench.o $(BUILD_DIR)/stl_slicer.o $(SHARED_LIB) | $(BIN_DIR) + $(CXX) $(CXXFLAGS) $(BUILD_DIR)/main_bench.o $(BUILD_DIR)/stl_slicer.o -L$(LIB_DIR) -lcudaslicer_offset -lcudart -Wl,-rpath,'$$ORIGIN/../lib' -o $@ + +$(BUILD_DIR) $(LIB_DIR) $(BIN_DIR): + mkdir -p $@ + +clean: + rm -rf $(BUILD_DIR) $(LIB_DIR) $(BIN_DIR) + +.PHONY: all clean diff --git a/README.md b/README.md new file mode 100644 index 0000000..551bf9e --- /dev/null +++ b/README.md @@ -0,0 +1,231 @@ +# CUDA Slicer: GPU-Accelerated Slicing & Polygon Offsetting Library + +[![Python Version](https://img.shields.io/badge/python-3.8%20%7C%203.9%20%7C%203.10%20%7C%203.11%20%7C%203.12-blue)](https://pypi.org/) +[![CUDA](https://img.shields.io/badge/CUDA-11.x%20%7C%2012.x-green.svg)](https://developer.nvidia.com/cuda-toolkit) +[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) + +A high-performance C/C++ and CUDA library with native Python bindings for **Linux & Windows** providing GPU-accelerated 2D/3D polygon offsetting, Jump Flooding Signed Distance Fields (JFA / SDF), and slicing for 3D printing, CAM, and computational geometry. + +--- + +## 🚀 Benchmark Summary (3DBenchy, 239 Layers @ 0.20 mm) + +Tested on **NVIDIA RTX A4500 Laptop GPU (16 GB VRAM)** vs **12th Gen Intel Core i7 CPU**: + +| Pipeline | Total Time (Mesh $\to$ Offset Contours) | Slicing / Offsetting Rate | Speedup vs PrusaSlicer | +| :--- | :---: | :---: | :---: | +| **100% GPU Direct Mesh Slicing (0.10 mm Pitch)** | **74.4 ms** | **3,212 layers/sec** | **57.47x Faster** | +| **100% GPU Direct Mesh Slicing (0.05 mm Pitch)** | **190.9 ms** | **1,252 layers/sec** | **22.55x Faster** | +| **Hybrid CPU Slicing + GPU JFA (0.05 mm Pitch)** | 739.4 ms | 323 layers/sec | 5.82x Faster | +| **PrusaSlicer CLI Headless (CPU Baseline)** | 4,305.7 ms | 55 layers/sec | 1.00x (Baseline) | + +> ⚡ **Key Optimizations Applied:** +> - **Stage 6 (Marching Squares Vectorization):** Reduced from **72.5 ms down to 35.4 ms (2.05x faster / 6,751 layer/s)** via zero-allocation discrete edge-table addressing, thread-local buffer reuse, and on-the-fly collinear point decimation (>95% vertex reduction). +> - **Stage 5 (D2H Transfer):** Reduced from **136.3 ms down to 4.5 ms (30x–48x faster)** via GPU 32-bit bit-packing + page-locked pinned host DMA (`cudaMallocHost`). +> - **Stage 3 (GPU JFA Distance Field):** Reduced from **118.9 ms down to 52.8 ms (2.25x faster)** via perimeter-band limited steps and integer arithmetic. + +--- + +## 🐍 Python Package Integration + +### 1. Installation + +#### A. Install from Source (Linux & Windows) +Make sure you have NVIDIA CUDA Toolkit (11.5+ or 12.x) and a C++ compiler (GCC on Linux, MSVC 2019/2022 on Windows) installed: +```bash +# Clone the repository +cd cudaslicer + +# Install directly with pip +pip install . +``` + +#### B. Build Wheel (.whl) +```bash +python -m pip install --upgrade build wheel +python -m build --wheel +# The wheel will be generated in dist/ +# e.g., dist/cudaslicer-1.0.0-cp310-cp310-linux_x86_64.whl +# or dist/cudaslicer-1.0.0-cp310-cp310-win_amd64.whl +``` + +--- + +### 2. Python API Quickstart + +#### A. Backend Selection & GPU Capability Check +```python +import cudaslicer + +# 1. Check if a CUDA GPU is present +if cudaslicer.has_cuda(): + print("CUDA GPU is active!") +else: + print("Running in multi-threaded CPU compatibility mode.") + +# 2. Control compute backend globally ('auto', 'cuda', 'cpu') +cudaslicer.set_backend("auto") # Auto: use CUDA if available, fallback to CPU +``` + +#### B. Offsetting 2D Polygon Contours (Single Layer or Multi-Layer) +```python +import numpy as np +import cudaslicer + +# Define a polygon (e.g. 50x50 mm square) +poly = np.array([ + [10.0, 10.0], + [60.0, 10.0], + [60.0, 60.0], + [10.0, 60.0] +], dtype=np.float32) + +# Generate 2 inward perimeters (at -0.5*w and -1.5*w) and infill mask (at -2.0*w) +# pitch_mm defaults to 0.10 mm, backend defaults to 'auto' +layer = cudaslicer.offset_polygons( + layers=[poly], + perimeters_count=2, + layer_width=0.45, + pitch_mm=0.10, + backend="auto" # 'auto', 'cuda', or 'cpu' +) + +# Unpack directly +perimeters, infill_mask = layer +print(f"Outer Perimeter 1 (-0.225mm): {len(perimeters[0])} loop(s)") +print(f"Inner Perimeter 2 (-0.675mm): {len(perimeters[1])} loop(s)") +print(f"Infill Region Mask (-0.900mm): shape {infill_mask.shape}, {infill_mask.sum()} px") +``` + +#### C. Offsetting 2D / 3D Occupancy Masks Directly +```python +import numpy as np +import cudaslicer + +# Create a 2D binary grid (e.g., 200 x 200 pixels) +mask = np.zeros((200, 200), dtype=np.uint8) +mask[40:160, 40:160] = 1 + +# Compute Euclidean perimeters and infill mask +perimeters, infill_mask = cudaslicer.offset_masks( + masks=mask, + perimeters_count=2, + layer_width=0.45, + pitch_mm=0.10, + backend="auto" +) +print("Perimeters shape:", perimeters.shape) # (2, 200, 200) +print("Infill mask shape:", infill_mask.shape) # (200, 200) +``` + +#### D. End-to-End Mesh Slicing with Custom First Layer Height +```python +import trimesh +import cudaslicer + +mesh = trimesh.load("benchy.stl") + +# 100% GPU Slicing & Perimeter Offsetting +layers = cudaslicer.slice_mesh( + mesh_or_triangles=mesh, + layer_height=0.20, + first_layer_height=0.25, # custom first layer height + perimeters_count=2, + layer_width=0.45, + pitch_mm=0.10 # default 0.10 +) + +print(f"Sliced {len(layers)} layers in GPU VRAM!") +layer_0 = layers[0] # First layer at Z = 0.25 mm +perims, infill_mask = layer_0 +print(f"Layer 0 (Z = {layer_0.z} mm): {len(perims)} perimeters, infill mask: {infill_mask.shape}") +``` + +--- + +## 🛠️ Building Standalone C/C++ Libraries (Windows & Linux) + +### Linux Build (`libcudaslicer_offset.so` & `.a`) +```bash +mkdir -p build && cd build +cmake .. -DCMAKE_BUILD_TYPE=Release +cmake --build . -j$(nproc) +``` +Outputs: +- Shared Library: `build/lib/libcudaslicer_offset.so` +- Static Library: `build/lib/libcudaslicer_offset.a` +- C Benchmark CLI: `build/bin/gpu_slicer_c_bench` + +### Windows Build (`cudaslicer_offset.dll` & `.lib`) +Open **x64 Native Tools Command Prompt for VS 2022** and run: +```cmd +mkdir build +cd build +cmake .. -G "Visual Studio 17 2022" -A x64 -DCMAKE_BUILD_TYPE=Release +cmake --build . --config Release +``` +Outputs: +- DLL: `build/lib/Release/cudaslicer_offset.dll` +- Import Lib: `build/lib/Release/cudaslicer_offset.lib` +- Benchmark: `build/bin/Release/gpu_slicer_c_bench.exe` + +--- + +## ⚙️ Cross-Platform C API Reference (`include/gpu_offset.h`) + +The public header provides a portable C API with `__declspec(dllexport)` / `__declspec(dllimport)` on Windows and default symbol visibility on Linux: + +```c +#include "gpu_offset.h" + +// 1. Initialize GPU +gpu_offset_init(0); + +// 2. Configure slicing & offsetting parameters +GpuOffsetConfig config = { + .pitch_mm = 0.10f, // default 0.10 mm + .pad_mm = 1.5f, + .perimeters_count = 2, // 2 perimeter passes + .layer_width = 0.45f, // 0.45 mm extrusion width + .first_layer_height = 0.25f, // 0.25 mm first layer height + .export_infill_mask = 1, // generate 2D binary infill mask + .use_gpu_rasterizer = 1 +}; + +// 3. 100% Direct GPU Mesh Slicing & Offsetting +GpuLayerResult* results = NULL; +size_t num_layers = 0; +GpuOffsetStats stats = {0}; + +gpu_offset_slice_mesh_and_offset( + triangles, num_triangles, + min_x, min_y, min_z, max_x, max_y, max_z, + 0.20f, 0.25f, &config, + &results, &num_layers, &stats +); + +// Access perimeters and infill mask +for (size_t l = 0; l < num_layers; ++l) { + float z = results[l].z_height; + size_t num_perims = results[l].num_perimeters; // 2 + GpuContour* infill_cnts = results[l].infill.contours; + uint8_t* infill_mask = results[l].infill.mask; +} + +// 4. Cleanup +gpu_offset_free_results(results, num_layers); +gpu_offset_cleanup(); +``` + +--- + +## 🔄 Automated CI/CD (GitHub Actions) + +- Building Linux Python wheels (`manylinux_2_35_x86_64`) across Python 3.9, 3.10, 3.11, 3.12. +- Building Windows x64 Python wheels (`win_amd64`) across Python 3.9, 3.10, 3.11, 3.12. +- Building standalone Linux `.so` and Windows `.dll` binaries as downloadable CI artifacts. + +--- + +## 📄 License +MIT License. Free for commercial and research use. diff --git a/benchy.stl b/benchy.stl new file mode 100644 index 0000000..5a51bb5 Binary files /dev/null and b/benchy.stl differ diff --git a/examples/c_example.c b/examples/c_example.c new file mode 100644 index 0000000..af9ec0e --- /dev/null +++ b/examples/c_example.c @@ -0,0 +1,87 @@ +#include "gpu_offset.h" +#include +#include +#include + +int main() { + printf("=== C API Example: Linking libcudaslicer_offset.so ===\n"); + + // 1. Initialize GPU + if (gpu_offset_init(0) != 0) { + fprintf(stderr, "Failed to initialize GPU\n"); + return 1; + } + + // 2. Create a test slice layer with a square polygon and a circular hole + GpuContour contours[2]; + + // Square exterior: 20mm x 20mm centered at origin + GpuPoint2D square_pts[4] = { + { -10.0f, -10.0f }, + { 10.0f, -10.0f }, + { 10.0f, 10.0f }, + { -10.0f, 10.0f } + }; + contours[0].points = square_pts; + contours[0].num_points = 4; + contours[0].is_hole = 0; + + // Diamond interior hole + GpuPoint2D hole_pts[4] = { + { 0.0f, -4.0f }, + { 4.0f, 0.0f }, + { 0.0f, 4.0f }, + { -4.0f, 0.0f } + }; + contours[1].points = hole_pts; + contours[1].num_points = 4; + contours[1].is_hole = 1; + + GpuSliceLayer layer; + layer.z_height = 0.20f; + layer.contours = contours; + layer.num_contours = 2; + + // 3. Configure Offsets + float offset_dists[2] = { 0.45f, 0.90f }; + GpuOffsetConfig config; + config.pitch_mm = 0.05f; + config.pad_mm = 1.0f; + config.offset_distances_mm = offset_dists; + config.num_offsets = 2; + config.use_gpu_rasterizer = 1; + + // 4. Run GPU Offset + GpuLayerResult* results = NULL; + GpuOffsetStats stats = {0}; + + int ret = gpu_offset_polygon_layers( + &layer, 1, + -10.0f, -10.0f, 10.0f, 10.0f, + &config, &results, &stats + ); + + if (ret == 0 && results) { + printf("GPU Offset Successful!\n"); + printf(" Rasterization: %.3f ms\n", stats.t_raster_ms); + printf(" H2D Transfer: %.3f ms\n", stats.t_h2d_ms); + printf(" GPU JFA Compute: %.3f ms\n", stats.t_gpu_jfa_ms); + printf(" D2H Transfer: %.3f ms\n", stats.t_d2h_ms); + printf(" Vectorization: %.3f ms\n", stats.t_vectorize_ms); + printf(" Total Time: %.3f ms\n", stats.t_total_pipeline_ms); + + for (size_t o = 0; o < results[0].num_perimeters; ++o) { + printf(" Perimeter %zu (Offset -%.2f mm): %zu vector loops\n", + o + 1, results[0].perimeters[o].offset_distance_mm, + results[0].perimeters[o].num_contours); + } + + gpu_offset_free_results(results, 1); + } else { + fprintf(stderr, "GPU Offset failed!\n"); + } + + gpu_offset_cleanup(); + printf("=== Done ===\n"); + return 0; +} diff --git a/examples/python_ctypes_example.py b/examples/python_ctypes_example.py new file mode 100644 index 0000000..2d3d50a --- /dev/null +++ b/examples/python_ctypes_example.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python3 +""" +Python ctypes wrapper demonstrating how to load and call libcudaslicer_offset.so +""" +import ctypes +import os +import numpy as np + +# Load shared library +lib_path = os.path.abspath(os.path.join(os.path.dirname(__file__), '..', 'lib', 'libcudaslicer_offset.so')) +lib = ctypes.CDLL(lib_path) + +# Define C structures +class GpuPoint2D(ctypes.Structure): + _fields_ = [("x", ctypes.c_float), ("y", ctypes.c_float)] + +class GpuContour(ctypes.Structure): + _fields_ = [ + ("points", ctypes.POINTER(GpuPoint2D)), + ("num_points", ctypes.c_size_t), + ("is_hole", ctypes.c_int) + ] + +class GpuSliceLayer(ctypes.Structure): + _fields_ = [ + ("contours", ctypes.POINTER(GpuContour)), + ("num_contours", ctypes.c_size_t), + ("z_height", ctypes.c_float) + ] + +class GpuOffsetPerimeter(ctypes.Structure): + _fields_ = [ + ("contours", ctypes.POINTER(GpuContour)), + ("num_contours", ctypes.c_size_t), + ("offset_distance_mm", ctypes.c_float) + ] + +class GpuLayerResult(ctypes.Structure): + _fields_ = [ + ("perimeters", ctypes.POINTER(GpuOffsetPerimeter)), + ("num_perimeters", ctypes.c_size_t), + ("z_height", ctypes.c_float) + ] + +class GpuOffsetConfig(ctypes.Structure): + _fields_ = [ + ("pitch_mm", ctypes.c_float), + ("pad_mm", ctypes.c_float), + ("offset_distances_mm", ctypes.POINTER(ctypes.c_float)), + ("num_offsets", ctypes.c_size_t), + ("use_gpu_rasterizer", ctypes.c_int) + ] + +class GpuOffsetStats(ctypes.Structure): + _fields_ = [ + ("t_mesh_slice_ms", ctypes.c_double), + ("t_raster_ms", ctypes.c_double), + ("t_h2d_ms", ctypes.c_double), + ("t_gpu_jfa_ms", ctypes.c_double), + ("t_gpu_thresh_ms", ctypes.c_double), + ("t_d2h_ms", ctypes.c_double), + ("t_vectorize_ms", ctypes.c_double), + ("t_total_pipeline_ms", ctypes.c_double), + ] + +# Setup function signatures +lib.gpu_offset_init.argtypes = [ctypes.c_int] +lib.gpu_offset_init.restype = ctypes.c_int + +lib.gpu_offset_cleanup.argtypes = [] +lib.gpu_offset_cleanup.restype = None + +lib.gpu_offset_polygon_layers.argtypes = [ + ctypes.POINTER(GpuSliceLayer), + ctypes.c_size_t, + ctypes.c_float, ctypes.c_float, ctypes.c_float, ctypes.c_float, + ctypes.POINTER(GpuOffsetConfig), + ctypes.POINTER(ctypes.POINTER(GpuLayerResult)), + ctypes.POINTER(GpuOffsetStats) +] +lib.gpu_offset_polygon_layers.restype = ctypes.c_int + +lib.gpu_offset_free_results.argtypes = [ctypes.POINTER(GpuLayerResult), ctypes.c_size_t] +lib.gpu_offset_free_results.restype = None + +def main(): + print(f"Loaded {lib_path}") + lib.gpu_offset_init(0) + + # Create a 20x20mm square + square_pts = (GpuPoint2D * 4)( + GpuPoint2D(-10.0, -10.0), + GpuPoint2D( 10.0, -10.0), + GpuPoint2D( 10.0, 10.0), + GpuPoint2D(-10.0, 10.0) + ) + cnt = GpuContour(square_pts, 4, 0) + layer = GpuSliceLayer((GpuContour * 1)(cnt), 1, 0.20) + + offsets = (ctypes.c_float * 2)(0.45, 0.90) + cfg = GpuOffsetConfig(0.05, 1.0, offsets, 2, 1) + + results_ptr = ctypes.POINTER(GpuLayerResult)() + stats = GpuOffsetStats() + + ret = lib.gpu_offset_polygon_layers( + ctypes.byref(layer), 1, + -10.0, -10.0, 10.0, 10.0, + ctypes.byref(cfg), + ctypes.byref(results_ptr), + ctypes.byref(stats) + ) + + if ret == 0 and results_ptr: + res = results_ptr[0] + print(f"C Library returned {res.num_perimeters} perimeters in {stats.t_total_pipeline_ms:.2f} ms") + for o in range(res.num_perimeters): + perim = res.perimeters[o] + print(f" Perimeter {o+1} (-{perim.offset_distance_mm:.2f} mm): {perim.num_contours} contours") + lib.gpu_offset_free_results(results_ptr, 1) + + lib.gpu_offset_cleanup() + +if __name__ == '__main__': + main() diff --git a/examples/python_quickstart.py b/examples/python_quickstart.py new file mode 100644 index 0000000..953a0a2 --- /dev/null +++ b/examples/python_quickstart.py @@ -0,0 +1,119 @@ +#!/usr/bin/env python3 +""" +CUDASlicer Python Quickstart & Integration Example +================================================= +Demonstrates how to use `cudaslicer` in a Python project with: +- first_layer_height +- perimeters_count and layer_width (Perimeter 1: -0.5w, Perimeter 2: -1.5w, etc.) +- default pitch_mm = 0.10 +- outputs: perimeters and infill_mask (-perimeters_count * layer_width) +""" + +import time +import numpy as np +import trimesh +import cudaslicer + + +def example_mask_offsetting(): + print("--- 1. GPU Mask Offsetting ---") + # Create a 200x200 test grid with a square + grid = np.zeros((200, 200), dtype=np.uint8) + grid[50:150, 50:150] = 1 + + t0 = time.perf_counter() + # Generate 2 perimeters and infill mask with layer_width=0.45 mm, default pitch_mm=0.10 mm + perimeters, infill_mask = cudaslicer.offset_masks( + masks=grid, + perimeters_count=2, + layer_width=0.45, + pitch_mm=0.10 # default is 0.10 + ) + t1 = time.perf_counter() + + print(f"Mask offset completed in {(t1 - t0)*1000:.3f} ms") + print(f"Original area (pixels): {grid.sum()}") + print(f"Perimeter 1 (-0.5*width = -0.225mm) area (pixels): {perimeters[0].sum()}") + print(f"Perimeter 2 (-1.5*width = -0.675mm) area (pixels): {perimeters[1].sum()}") + print(f"Infill mask (-2.0*width = -0.900mm) area (pixels): {infill_mask.sum()}\n") + + +def example_polygon_offsetting(): + print("--- 2. Single Layer Polygon Offsetting ---") + # Define an L-shaped polygon + l_poly = np.array([ + [0.0, 0.0], + [30.0, 0.0], + [30.0, 10.0], + [10.0, 10.0], + [10.0, 30.0], + [0.0, 30.0] + ], dtype=np.float32) + + t0 = time.perf_counter() + layer_res = cudaslicer.offset_polygons( + layers=[l_poly], + perimeters_count=2, + layer_width=0.45, + pitch_mm=0.10, + pad_mm=1.0, + return_infill_mask=True + ) + t1 = time.perf_counter() + + perimeters, infill_mask = layer_res # Unpack directly! + print(f"Polygon offsetting completed in {(t1 - t0)*1000:.3f} ms") + print(f" Perimeter 1 (-0.225mm): {len(perimeters[0])} loop(s), vertex count: {[len(c) for c in perimeters[0]]}") + print(f" Perimeter 2 (-0.675mm): {len(perimeters[1])} loop(s), vertex count: {[len(c) for c in perimeters[1]]}") + print(f" Infill contours (-0.900mm): {len(layer_res.infill_contours)} loop(s), vertex count: {[len(c) for c in layer_res.infill_contours]}") + if infill_mask is not None: + print(f" Infill mask raster shape: {infill_mask.shape}, active pixels: {infill_mask.sum()}") + print() + + +def example_mesh_slicing_and_offset(): + print("--- 3. 100% GPU Direct 3D Mesh Slicing & Offsetting ---") + mesh = trimesh.load("benchy.stl") + print(f"Loaded mesh with {len(mesh.faces)} triangles. Running 100% GPU Slicing pipeline...") + + t0 = time.perf_counter() + layers = cudaslicer.slice_mesh_and_offset( + mesh_or_triangles=mesh, + layer_height=0.20, + first_layer_height=0.25, # custom first layer height + perimeters_count=2, + layer_width=0.45, + pitch_mm=0.10, # default 0.10 + pad_mm=1.5, + return_infill_mask=True + ) + t1 = time.perf_counter() + + total_ms = (t1 - t0) * 1000 + print(f"100% GPU Slicing & Offsetting done in {total_ms:.2f} ms ({len(layers)} layers)") + print(f"Throughput: {len(layers) / (total_ms / 1000.0):.1f} layers/second!") + + # Inspect first layer and a middle layer + layer_0 = layers[0] + print(f"\nLayer 0 (Z = {layer_0.z:.2f} mm - First Layer Height):") + perims_0, infill_mask_0 = layer_0 + print(f" Perimeters: {len(perims_0)} passes (Outer: {len(perims_0[0])} loops, Inner: {len(perims_0[1])} loops)") + print(f" Infill Mask: {infill_mask_0.shape}, Infill Contours: {len(layer_0.infill_contours)} loop(s)") + + mid_idx = len(layers) // 2 + layer_mid = layers[mid_idx] + print(f"Layer {mid_idx} (Z = {layer_mid.z:.2f} mm):") + perims_mid, infill_mask_mid = layer_mid + print(f" Perimeters: {len(perims_mid)} passes (Outer: {len(perims_mid[0])} loops, Inner: {len(perims_mid[1])} loops)") + print(f" Infill Mask: {infill_mask_mid.shape}, Infill Contours: {len(layer_mid.infill_contours)} loop(s)\n") + + +if __name__ == "__main__": + if not cudaslicer.is_available(): + print("Error: CUDA GPU is not available.") + exit(1) + + print(f"CUDASlicer Version: {cudaslicer.__version__}") + example_mask_offsetting() + example_polygon_offsetting() + example_mesh_slicing_and_offset() diff --git a/include/gpu_offset.h b/include/gpu_offset.h new file mode 100644 index 0000000..1bf67c6 --- /dev/null +++ b/include/gpu_offset.h @@ -0,0 +1,212 @@ +#ifndef GPU_OFFSET_H +#define GPU_OFFSET_H + +#ifdef _WIN32 + #ifdef GPU_OFFSET_EXPORTS + #define GPU_OFFSET_API __declspec(dllexport) + #else + #define GPU_OFFSET_API __declspec(dllimport) + #endif +#else + #define GPU_OFFSET_API __attribute__((visibility("default"))) +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +#include +#include + +/** + * @brief Status and Error Codes + */ +#define GPU_OFFSET_SUCCESS 0 +#define GPU_OFFSET_ERR_GENERIC -1 +#define GPU_OFFSET_ERR_NO_CUDA -2 +#define GPU_OFFSET_ERR_INVALID_ARG -3 + +/** + * @brief Compute Backend Modes + */ +typedef enum { + GPU_OFFSET_BACKEND_AUTO = 0, // Auto: use CUDA if available, fallback to CPU + GPU_OFFSET_BACKEND_CUDA = 1, // Require CUDA: returns GPU_OFFSET_ERR_NO_CUDA if unavailable + GPU_OFFSET_BACKEND_CPU = 2 // Compatibility mode: force multi-threaded CPU execution +} GpuOffsetBackend; + +/** + * @brief 2D Point structure (floating-point millimeters) + */ +typedef struct { + float x; + float y; +} GpuPoint2D; + +/** + * @brief 2D Polygon Contour Loop + */ +typedef struct { + GpuPoint2D* points; + size_t num_points; + int is_hole; // 0 = exterior boundary, 1 = interior hole +} GpuContour; + +/** + * @brief Sliced Layer containing one or more planar polygon loops + */ +typedef struct { + GpuContour* contours; + size_t num_contours; + float z_height; +} GpuSliceLayer; + +/** + * @brief Resulting Offset Contours for a single offset perimeter + */ +typedef struct { + GpuContour* contours; + size_t num_contours; + float offset_distance_mm; +} GpuOffsetPerimeter; + +/** + * @brief Infill region mask and boundary contours for a single sliced layer + */ +typedef struct { + GpuContour* contours; // Infill boundary contours (at -perimeters_count * layer_width) + size_t num_contours; + uint8_t* mask; // 2D binary raster mask (width * height uint8), or NULL if not exported + int width; + int height; + float offset_distance_mm; // Offset distance for infill (perimeters_count * layer_width) +} GpuInfillRegion; + +/** + * @brief All offset perimeters and infill mask for a single sliced layer + */ +typedef struct { + GpuOffsetPerimeter* perimeters; // Array of perimeters (length = num_perimeters) + size_t num_perimeters; // perimeters_count + GpuInfillRegion infill; // Infill boundary contours & 2D mask + float z_height; +} GpuLayerResult; + +/** + * @brief Configuration parameters for Offsetting & Slicing Engine + */ +typedef struct { + float pitch_mm; // Pixel raster pitch (default: 0.10 mm) + float pad_mm; // Bounding box padding margin in mm (default: 1.5 mm) + int perimeters_count; // Number of perimeter passes (e.g. 2, default: 2) + float layer_width; // Line/extrusion width in mm (e.g. 0.45 mm, default: 0.45 mm) + const float* offset_distances_mm; // Optional custom offset distances (if NULL, auto-generated from perimeters_count & layer_width) + size_t num_offsets; // Total offsets (auto-set if offset_distances_mm is NULL) + 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) + int backend; // GpuOffsetBackend: 0 = AUTO, 1 = CUDA, 2 = CPU (default: 0) +} GpuOffsetConfig; + +/** + * @brief Microsecond-accurate timing breakdown + */ +typedef struct { + double t_mesh_slice_ms; + double t_raster_ms; + double t_h2d_ms; + double t_gpu_jfa_ms; + double t_gpu_thresh_ms; + double t_d2h_ms; + double t_vectorize_ms; + double t_total_pipeline_ms; +} GpuOffsetStats; + +/** + * @brief Returns 1 if a functional CUDA-capable device is detected, 0 otherwise. + */ +GPU_OFFSET_API int gpu_offset_has_cuda(void); + +/** + * @brief Set global default backend (0 = AUTO, 1 = CUDA, 2 = CPU). + */ +GPU_OFFSET_API int gpu_offset_set_backend(int backend); + +/** + * @brief Get global default backend (0 = AUTO, 1 = CUDA, 2 = CPU). + */ +GPU_OFFSET_API int gpu_offset_get_backend(void); + +/** + * @brief Initialize GPU device and allocate CUDA stream / warm up kernels. + */ +GPU_OFFSET_API int gpu_offset_init(int device_id); + +/** + * @brief Core Distance Field & Offsetting on raw 2D/3D binary masks (CUDA or CPU fallback). + */ +GPU_OFFSET_API int gpu_offset_compute_masks( + const uint8_t* h_masks, + int width, + int height, + int num_layers, + float pitch_mm, + const float* offset_distances, + int num_offsets, + uint8_t* h_out_masks, + GpuOffsetStats* stats +); + +/** + * @brief 3D Triangle structure for direct mesh slicing (vertices in millimeters) + */ +typedef struct { + float v0[3]; + float v1[3]; + float v2[3]; +} GpuTriangle3D; + +/** + * @brief End-to-end polygon offsetting: Polygons -> Raster -> Distance Field -> Vector Contours & Infill + */ +GPU_OFFSET_API int gpu_offset_polygon_layers( + const GpuSliceLayer* in_layers, + size_t num_layers, + float min_x, float min_y, float max_x, float max_y, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + GpuOffsetStats* stats +); + +/** + * @brief End-to-End Mesh Slicing & Perimeter Offsetting: + * 3D Triangles -> Mesh Slicer -> Distance Field -> Vector Contours & Infill. + */ +GPU_OFFSET_API int gpu_offset_slice_mesh_and_offset( + 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, + float first_layer_height, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + size_t* out_num_layers, + GpuOffsetStats* stats +); + +/** + * @brief Free memory allocated for GpuLayerResult + */ +GPU_OFFSET_API void gpu_offset_free_results(GpuLayerResult* results, size_t num_layers); + +/** + * @brief Release GPU context and cached buffers + */ +GPU_OFFSET_API void gpu_offset_cleanup(); + +#ifdef __cplusplus +} +#endif + +#endif // GPU_OFFSET_H diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..718ec32 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,13 @@ +[build-system] +requires = ["setuptools>=61.0", "wheel", "pybind11>=2.10.0", "cmake>=3.18"] +build-backend = "setuptools.build_meta" + +[project] +name = "cudaslicer" +version = "1.0.0" +description = "GPU-Accelerated 2D/3D Slicing & Polygon Offsetting Library for Linux & Windows" +readme = "README.md" +requires-python = ">=3.8" +dependencies = [ + "numpy>=1.20.0" +] diff --git a/python/cudaslicer/__init__.py b/python/cudaslicer/__init__.py new file mode 100644 index 0000000..cd18bda --- /dev/null +++ b/python/cudaslicer/__init__.py @@ -0,0 +1,363 @@ +""" +CUDASlicer: High-Performance Slicing & Polygon Offsetting Library +================================================================ +Provides ultra-fast GPU distance-field (SDF/JFA) & multi-threaded CPU +polygon offsetting, binary mask transformations, and 3D mesh slicing +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" + +# Try importing the compiled C++/CUDA pybind11 module +try: + from . import _cudaslicer_backend as _backend + _HAS_BACKEND = True +except ImportError: + _backend = None + _HAS_BACKEND = False + + +class LayerResult(dict): + """ + Dictionary holding slicing and offsetting results for a single layer. + Supports dictionary access, attribute access, and tuple unpacking: + `perimeters, infill_mask = layer` + """ + def __getattr__(self, key: str) -> Any: + try: + return self[key] + except KeyError: + raise AttributeError(f"'LayerResult' object has no attribute '{key}'") + + def __setattr__(self, key: str, value: Any) -> None: + self[key] = value + + def __iter__(self): + # Enables direct unpacking: perimeters, infill_mask = layer + yield self.get("perimeters", []) + yield self.get("infill_mask", None) + + +def has_cuda() -> bool: + """ + Check if a functional CUDA GPU device is available on this system. + + Returns + ------- + bool + True if CUDA-capable GPU is available, False otherwise. + """ + if not _HAS_BACKEND or _backend is None: + return False + try: + return bool(_backend.has_cuda()) + except Exception: + return False + + +def is_cuda_available() -> bool: + """Alias for has_cuda(). Returns True if CUDA GPU is functional.""" + return has_cuda() + + +def is_available() -> bool: + """ + Check if the CUDASlicer backend (CUDA or CPU fallback) is available. + + Returns + ------- + bool + True if the native backend is functional. + """ + if not _HAS_BACKEND or _backend is None: + return False + try: + ret = _backend.init(0) + return ret == 0 + except Exception: + return False + + +def get_cuda_device_count() -> int: + """Return the number of available CUDA GPUs.""" + return 1 if has_cuda() else 0 + + +def set_backend(backend: str = "auto") -> None: + """ + Set the global default compute backend. + + Parameters + ---------- + backend : str, default='auto' + Execution backend: 'auto' (use CUDA if available, fallback to CPU), + 'cuda' (require CUDA; throws error if unavailable), or 'cpu' (compatibility mode). + """ + if not _HAS_BACKEND or _backend is None: + raise RuntimeError("CUDASlicer native backend extension is not loaded.") + _backend.set_backend(str(backend).lower()) + + +def get_backend() -> str: + """ + Get the active global compute backend ('auto', 'cuda', or 'cpu'). + """ + if not _HAS_BACKEND or _backend is None: + return "none" + return _backend.get_backend() + + +def offset_masks( + masks: np.ndarray, + perimeters_count: int = 2, + layer_width: float = 0.45, + pitch_mm: float = 0.10, + offset_distances: Optional[Sequence[float]] = None, + backend: Optional[str] = "auto" +) -> Tuple[np.ndarray, np.ndarray]: + """ + Compute Euclidean offset perimeter masks and infill mask (CUDA or CPU). + + Parameters + ---------- + masks : np.ndarray + 2D (H, W) or 3D (N, H, W) binary occupancy grid (uint8 or bool). + perimeters_count : int, default=2 + Number of perimeter passes (Perimeter 1 at -0.5*width, Perimeter 2 at -1.5*width, etc.). + layer_width : float, default=0.45 + Extrusion line width in millimeters. + pitch_mm : float, default=0.10 + Grid pixel pitch in millimeters. + offset_distances : sequence of float, optional + Custom offset distances. If provided, overrides perimeters_count and layer_width. + backend : str, optional, default='auto' + Compute backend: 'auto', 'cuda', or 'cpu'. + + Returns + ------- + perimeters_masks : np.ndarray + Shape (perimeters_count, H, W) or (perimeters_count, N, H, W). + infill_mask : np.ndarray + Shape (H, W) or (N, H, W) at -perimeters_count * layer_width. + """ + if not _HAS_BACKEND or _backend is None: + raise RuntimeError("CUDASlicer backend extension (_cudaslicer_backend) is not available.") + + masks_uint8 = np.ascontiguousarray(masks, dtype=np.uint8) + offsets_float = [float(d) for d in offset_distances] if offset_distances is not None else [] + backend_str = str(backend) if backend is not None else "auto" + + out_masks = _backend.offset_masks( + masks_uint8, + int(perimeters_count), + float(layer_width), + float(pitch_mm), + offsets_float, + backend_str + ) + + if offset_distances is not None: + return out_masks + + # Separate perimeters and infill mask + perimeters_masks = out_masks[:perimeters_count] + infill_mask = out_masks[perimeters_count] + return perimeters_masks, infill_mask + + +def offset_polygons( + layers: Union[List[np.ndarray], List[Dict[str, Any]]], + perimeters_count: int = 2, + layer_width: float = 0.45, + pitch_mm: float = 0.10, + pad_mm: float = 1.5, + return_infill_mask: bool = True, + offset_distances: Optional[Sequence[float]] = None, + backend: Optional[str] = "auto" +) -> Union[LayerResult, List[LayerResult]]: + """ + Offset 2D polygon loops across one or more layers (CUDA or CPU fallback). + + Parameters + ---------- + layers : list of np.ndarray OR list of dict + Input layer data. Can be: + 1. Single layer with a list of (K, 2) NumPy polygon contours: + `[poly1, poly2, ...]` + 2. Multi-layer list of layer dicts: + `[{"z": 0.2, "contours": [poly1, ...]}, {"z": 0.4, "contours": [...]}]` + perimeters_count : int, default=2 + Number of perimeter passes (Perimeter 1: -0.5*width, Perimeter 2: -1.5*width, etc.). + layer_width : float, default=0.45 + Extrusion line width in millimeters. + pitch_mm : float, default=0.10 + Pixel raster pitch in millimeters. + pad_mm : float, default=1.5 + Bounding box padding around the contours in millimeters. + return_infill_mask : bool, default=True + Whether to export 2D binary infill occupancy grid. + offset_distances : sequence of float, optional + Custom offset distances. + backend : str, optional, default='auto' + Compute backend: 'auto', 'cuda', or 'cpu'. + + Returns + ------- + LayerResult or list of LayerResult + If single layer: returns single LayerResult object. + If multi-layer: returns list of LayerResult objects. + """ + if not _HAS_BACKEND or _backend is None: + raise RuntimeError("CUDASlicer backend extension (_cudaslicer_backend) is not available.") + + if not layers: + return [] + + is_single_layer_format = False + if isinstance(layers[0], (np.ndarray, list, tuple)) and not isinstance(layers[0], dict): + is_single_layer_format = True + formatted_layers = [{"z": 0.0, "contours": [np.ascontiguousarray(c, dtype=np.float32) for c in layers]}] + else: + formatted_layers = [] + for l in layers: + z_val = float(l.get("z", l.get("z_height", 0.0))) + cnts = [np.ascontiguousarray(c, dtype=np.float32) for c in l.get("contours", [])] + formatted_layers.append({"z": z_val, "contours": cnts}) + + all_pts = [] + for l in formatted_layers: + for c in l["contours"]: + if len(c) >= 3: + all_pts.append(c) + + if not all_pts: + empty_res = LayerResult({"z": 0.0, "perimeters": [[] for _ in range(perimeters_count)], "infill_contours": [], "infill_mask": None}) + return empty_res if is_single_layer_format else [] + + stacked = np.concatenate(all_pts, axis=0) + min_x, min_y = float(stacked[:, 0].min()), float(stacked[:, 1].min()) + max_x, max_y = float(stacked[:, 0].max()), float(stacked[:, 1].max()) + + offsets_float = [float(d) for d in offset_distances] if offset_distances is not None else [] + backend_str = str(backend) if backend is not None else "auto" + + raw_results = _backend.offset_polygon_layers( + formatted_layers, min_x, min_y, max_x, max_y, + int(perimeters_count), float(layer_width), + float(pitch_mm), float(pad_mm), + bool(return_infill_mask), offsets_float, + backend_str + ) + + layer_results = [LayerResult(r) for r in raw_results] + if is_single_layer_format and layer_results: + return layer_results[0] + + return layer_results + + +def slice_mesh_and_offset( + mesh_or_triangles: Any, + layer_height: float = 0.20, + first_layer_height: Optional[float] = None, + perimeters_count: int = 2, + layer_width: float = 0.45, + pitch_mm: float = 0.10, + pad_mm: float = 1.5, + return_infill_mask: bool = True, + offset_distances: Optional[Sequence[float]] = None, + backend: Optional[str] = "auto" +) -> List[LayerResult]: + """ + 3D mesh slicing, rasterization, distance field, and perimeter offsetting (CUDA or CPU). + + Parameters + ---------- + mesh_or_triangles : trimesh.Trimesh OR np.ndarray + Either a trimesh.Trimesh object, or a NumPy array of shape (N, 3, 3) or (N, 9) of triangle vertices. + 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 or <= 0. + perimeters_count : int, default=2 + Number of perimeter passes (Perimeter 1: -0.5*width, Perimeter 2: -1.5*width, etc.). + layer_width : float, default=0.45 + Extrusion line width in millimeters. + pitch_mm : float, default=0.10 + Grid pixel pitch in millimeters. + pad_mm : float, default=1.5 + Bounding box padding margin in millimeters. + return_infill_mask : bool, default=True + Whether to generate the 2D binary infill occupancy grid. + offset_distances : sequence of float, optional + Custom offset distances. + backend : str, optional, default='auto' + Compute backend: 'auto' (CUDA if available, else CPU), 'cuda' (require CUDA), or 'cpu'. + + Returns + ------- + list of LayerResult + List of LayerResult dictionaries containing: + - 'z': Layer Z height + - 'perimeters': List of perimeters (each containing list of (K, 2) contour arrays) + - 'infill_contours': Infill boundary contour loops + - 'infill_mask': 2D binary NumPy array (H, W) of the infill region + """ + if not _HAS_BACKEND or _backend is None: + raise RuntimeError("CUDASlicer backend extension (_cudaslicer_backend) is not available.") + + # Handle Trimesh object + 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) + offsets_float = [float(d) for d in offset_distances] if offset_distances is not None else [] + backend_str = str(backend) if backend is not None else "auto" + + raw_results = _backend.slice_mesh_and_offset( + tris, min_x, min_y, min_z, max_x, max_y, max_z, + float(layer_height), float(h_first), + int(perimeters_count), float(layer_width), + float(pitch_mm), float(pad_mm), + bool(return_infill_mask), offsets_float, + backend_str + ) + + return [LayerResult(r) for r in raw_results] + + +# Alias for convenience +slice_mesh = slice_mesh_and_offset + +__all__ = [ + "offset_masks", + "offset_polygons", + "slice_mesh_and_offset", + "slice_mesh", + "LayerResult", + "is_available", + "has_cuda", + "is_cuda_available", + "set_backend", + "get_backend", + "get_cuda_device_count", + "__version__", +] diff --git a/python/cudaslicer/bindings.cpp b/python/cudaslicer/bindings.cpp new file mode 100644 index 0000000..9c9e60b --- /dev/null +++ b/python/cudaslicer/bindings.cpp @@ -0,0 +1,420 @@ +#include +#include +#include +#include "gpu_offset.h" +#include +#include +#include + +namespace py = pybind11; + +static int parse_backend_str(const std::string& backend) { + if (backend.empty() || backend == "auto" || backend == "AUTO") return GPU_OFFSET_BACKEND_AUTO; + if (backend == "cuda" || backend == "gpu" || backend == "CUDA" || backend == "GPU") return GPU_OFFSET_BACKEND_CUDA; + if (backend == "cpu" || backend == "compatibility" || backend == "CPU") return GPU_OFFSET_BACKEND_CPU; + throw std::invalid_argument("Invalid backend '" + backend + "'. Allowed options: 'auto', 'cuda', 'cpu'."); +} + +// 1. Direct NumPy Mask Offsetting +py::array_t py_offset_masks( + py::array_t in_masks, + int perimeters_count, + float layer_width, + float pitch_mm, + const std::vector& offset_distances, + const std::string& backend +) { + py::buffer_info buf = in_masks.request(); + int 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 = (int)buf.shape[0]; + height = (int)buf.shape[1]; + width = (int)buf.shape[2]; + } else { + throw std::runtime_error("Input mask must be a 2D (H, W) or 3D (N, H, W) numpy array."); + } + + std::vector resolved_offsets; + if (!offset_distances.empty()) { + resolved_offsets = offset_distances; + } else { + int P = (perimeters_count >= 0) ? perimeters_count : 2; + float W = (layer_width > 0.0f) ? layer_width : 0.45f; + 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; + } + + int num_offsets = (int)resolved_offsets.size(); + if (num_offsets == 0) { + throw std::runtime_error("offset distances cannot be empty."); + } + + std::vector out_shape; + if (buf.ndim == 2) { + out_shape = { (ssize_t)num_offsets, (ssize_t)height, (ssize_t)width }; + } else { + out_shape = { (ssize_t)num_offsets, (ssize_t)num_layers, (ssize_t)height, (ssize_t)width }; + } + + py::array_t out_masks(out_shape); + py::buffer_info out_buf = out_masks.request(); + + const uint8_t* ptr_in = static_cast(buf.ptr); + uint8_t* ptr_out = static_cast(out_buf.ptr); + + int req_backend = parse_backend_str(backend); + int prev_backend = gpu_offset_get_backend(); + if (!backend.empty()) { + gpu_offset_set_backend(req_backend); + } + + GpuOffsetStats stats = {0}; + int ret = gpu_offset_compute_masks( + ptr_in, width, height, num_layers, pitch_mm, + resolved_offsets.data(), num_offsets, ptr_out, &stats + ); + + if (!backend.empty()) { + gpu_offset_set_backend(prev_backend); + } + + if (ret == GPU_OFFSET_ERR_NO_CUDA) { + throw std::runtime_error("CUDA backend was requested, but no functional CUDA-capable device was found."); + } else if (ret != GPU_OFFSET_SUCCESS) { + throw std::runtime_error("Offset computation failed in compute backend."); + } + + return out_masks; +} + +// 2. Direct Python Polygon Offsetting +py::list py_offset_polygon_layers( + py::list py_layers, + float min_x, float min_y, float max_x, float max_y, + int perimeters_count, + float layer_width, + float pitch_mm, + float pad_mm, + bool export_infill_mask, + const std::vector& offset_distances, + const std::string& backend +) { + size_t num_layers = py_layers.size(); + if (num_layers == 0) return py::list(); + + std::vector c_layers(num_layers); + std::vector> c_contours_pool(num_layers); + std::vector>> c_points_pool(num_layers); + + for (size_t l = 0; l < num_layers; ++l) { + py::dict layer_dict = py_layers[l].cast(); + float z_h = layer_dict.contains("z") ? layer_dict["z"].cast() : 0.0f; + py::list cnt_list = layer_dict["contours"].cast(); + + size_t num_cnts = cnt_list.size(); + c_layers[l].z_height = z_h; + c_layers[l].num_contours = num_cnts; + c_contours_pool[l].resize(num_cnts); + c_points_pool[l].resize(num_cnts); + + for (size_t c = 0; c < num_cnts; ++c) { + py::array_t arr = cnt_list[c].cast>(); + py::buffer_info buf = arr.request(); + size_t n_pts = buf.shape[0]; + + auto& pts = c_points_pool[l][c]; + pts.resize(n_pts); + const float* raw_data = static_cast(buf.ptr); + + for (size_t p = 0; p < n_pts; ++p) { + pts[p].x = raw_data[p * 2 + 0]; + pts[p].y = raw_data[p * 2 + 1]; + } + + c_contours_pool[l][c].points = pts.data(); + c_contours_pool[l][c].num_points = n_pts; + c_contours_pool[l][c].is_hole = 0; + } + c_layers[l].contours = c_contours_pool[l].data(); + } + + GpuOffsetConfig config = {0}; + config.pitch_mm = pitch_mm; + config.pad_mm = pad_mm; + config.perimeters_count = perimeters_count; + config.layer_width = layer_width; + config.export_infill_mask = export_infill_mask ? 1 : 0; + config.use_gpu_rasterizer = 1; + config.backend = parse_backend_str(backend); + + if (!offset_distances.empty()) { + config.offset_distances_mm = offset_distances.data(); + config.num_offsets = offset_distances.size(); + } else { + config.offset_distances_mm = nullptr; + config.num_offsets = 0; + } + + GpuLayerResult* results = nullptr; + GpuOffsetStats stats = {0}; + + int ret = gpu_offset_polygon_layers( + c_layers.data(), num_layers, + min_x, min_y, max_x, max_y, + &config, &results, &stats + ); + + if (ret == GPU_OFFSET_ERR_NO_CUDA) { + throw std::runtime_error("CUDA backend was requested, but no functional CUDA-capable device was found."); + } else if (ret != GPU_OFFSET_SUCCESS || !results) { + throw std::runtime_error("Polygon offsetting failed in compute backend."); + } + + py::list py_layer_results; + for (size_t l = 0; l < num_layers; ++l) { + py::dict py_layer_dict; + py_layer_dict["z"] = results[l].z_height; + + // 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, num_layers); + return py_layer_results; +} + +// 3. Direct 3D Mesh Slicing & Offsetting +py::list py_slice_mesh_and_offset( + 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, + float first_layer_height, + int perimeters_count, + float layer_width, + float pitch_mm, + float pad_mm, + bool export_infill_mask, + const std::vector& offset_distances, + const std::string& backend +) { + 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) return py::list(); + + const GpuTriangle3D* tris_ptr = static_cast(buf.ptr); + + GpuOffsetConfig config = {0}; + config.pitch_mm = pitch_mm; + config.pad_mm = pad_mm; + config.perimeters_count = perimeters_count; + config.layer_width = layer_width; + config.first_layer_height = first_layer_height; + config.export_infill_mask = export_infill_mask ? 1 : 0; + config.use_gpu_rasterizer = 1; + config.backend = parse_backend_str(backend); + + if (!offset_distances.empty()) { + config.offset_distances_mm = offset_distances.data(); + config.num_offsets = offset_distances.size(); + } else { + config.offset_distances_mm = nullptr; + config.num_offsets = 0; + } + + GpuLayerResult* results = nullptr; + size_t out_num_layers = 0; + GpuOffsetStats stats = {0}; + + int ret = gpu_offset_slice_mesh_and_offset( + tris_ptr, num_tris, + min_x, min_y, min_z, max_x, max_y, max_z, + layer_height, first_layer_height, &config, + &results, &out_num_layers, &stats + ); + + if (ret == GPU_OFFSET_ERR_NO_CUDA) { + throw std::runtime_error("CUDA backend was requested, but no functional CUDA-capable device was found."); + } else if (ret != GPU_OFFSET_SUCCESS || !results) { + throw std::runtime_error("Slice mesh and offset failed in compute backend."); + } + + 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; + + // 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; +} + +PYBIND11_MODULE(_cudaslicer_backend, m) { + m.doc() = "GPU/CPU Slicing & SDF/JFA Offsetting C++/CUDA Python 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."); + + m.def("set_backend", [](const std::string& b) { + int code = gpu_offset_set_backend(parse_backend_str(b)); + if (code != GPU_OFFSET_SUCCESS) throw std::invalid_argument("Failed to set backend."); + }, py::arg("backend") = "auto", "Set global compute backend ('auto', 'cuda', 'cpu')."); + + m.def("get_backend", []() -> std::string { + int b = gpu_offset_get_backend(); + if (b == GPU_OFFSET_BACKEND_CUDA) return "cuda"; + if (b == GPU_OFFSET_BACKEND_CPU) return "cpu"; + return "auto"; + }, "Get global compute backend ('auto', 'cuda', 'cpu')."); + + m.def("init", &gpu_offset_init, py::arg("device_id") = 0, "Initialize CUDA runtime and allocate stream."); + m.def("cleanup", &gpu_offset_cleanup, "Release backend resources."); + + m.def("offset_masks", &py_offset_masks, + py::arg("in_masks"), + py::arg("perimeters_count") = 2, + py::arg("layer_width") = 0.45f, + py::arg("pitch_mm") = 0.10f, + py::arg("offset_distances") = std::vector(), + py::arg("backend") = "auto", + "Offset 2D (H,W) or 3D (N,H,W) binary occupancy numpy masks."); + + m.def("offset_polygon_layers", &py_offset_polygon_layers, + py::arg("py_layers"), py::arg("min_x"), py::arg("min_y"), + py::arg("max_x"), py::arg("max_y"), + py::arg("perimeters_count") = 2, + py::arg("layer_width") = 0.45f, + py::arg("pitch_mm") = 0.10f, + py::arg("pad_mm") = 1.5f, + py::arg("export_infill_mask") = true, + py::arg("offset_distances") = std::vector(), + py::arg("backend") = "auto", + "Compute polygon perimeter offsets and infill mask."); + + m.def("slice_mesh_and_offset", &py_slice_mesh_and_offset, + 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") = 0.20f, + py::arg("first_layer_height") = 0.0f, + py::arg("perimeters_count") = 2, + py::arg("layer_width") = 0.45f, + py::arg("pitch_mm") = 0.10f, + py::arg("pad_mm") = 1.5f, + py::arg("export_infill_mask") = true, + py::arg("offset_distances") = std::vector(), + py::arg("backend") = "auto", + "Mesh slicing, perimeter offsetting, and infill mask generation."); +} diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..381c82a --- /dev/null +++ b/setup.py @@ -0,0 +1,62 @@ +import os +import sys +import platform +import subprocess +from setuptools import setup, find_packages, Extension +from setuptools.command.build_ext import build_ext + +class CMakeExtension(Extension): + def __init__(self, name, sourcedir=""): + Extension.__init__(self, name, sources=[]) + self.sourcedir = os.path.abspath(sourcedir) + +class CMakeBuild(build_ext): + def build_extension(self, ext): + extdir = os.path.abspath(os.path.dirname(self.get_ext_fullpath(ext.name))) + cmake_args = [ + f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY={extdir}", + f"-DPYTHON_EXECUTABLE={sys.executable}", + f"-DPython3_EXECUTABLE={sys.executable}", + ] + + cfg = "Debug" if self.debug else "Release" + build_args = ["--config", cfg] + + if platform.system() == "Windows": + cmake_args += [ + f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{cfg.upper()}={extdir}", + "-A", "x64" + ] + else: + cmake_args += [f"-DCMAKE_BUILD_TYPE={cfg}"] + build_args += ["--", f"-j{os.cpu_count() or 4}"] + + build_temp = os.path.join(self.build_temp, ext.name) + os.makedirs(build_temp, exist_ok=True) + + subprocess.check_call(["cmake", ext.sourcedir] + cmake_args, cwd=build_temp) + subprocess.check_call(["cmake", "--build", "."] + build_args, cwd=build_temp) + +setup( + name="cudaslicer", + version="1.0.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", + 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 new file mode 100644 index 0000000..6fb39f0 --- /dev/null +++ b/src/cpu_offset.cpp @@ -0,0 +1,766 @@ +#include "cpu_offset.h" +#include "vectorizer.h" +#include +#include +#include +#include +#include +#include +#include + +// ----------------------------------------------------------------------------- +// Fast Separable 2D Exact Euclidean Distance Transform on Binary Grid +// ----------------------------------------------------------------------------- + +// Pass 1 on binary column: fast forward + backward scan (O(N) exact) +static inline void edt_1d_binary_col(const uint8_t* col_mask, int height, int stride, float* col_g) { + const float INF = 1e12f; + // Forward pass + int last = -1000000; + for (int y = 0; y < height; ++y) { + if (col_mask[(size_t)y * stride] == 0) { + last = y; + } + float dy = (float)(y - last); + col_g[y] = (last >= 0) ? (dy * dy) : INF; + } + // Backward pass + last = 1000000; + for (int y = height - 1; y >= 0; --y) { + if (col_mask[(size_t)y * stride] == 0) { + last = y; + } + float dy = (float)(last - y); + float d2 = (last < height) ? (dy * dy) : INF; + if (d2 < col_g[y]) { + col_g[y] = d2; + } + } +} + +// Pass 2 on rows: Parabolic Envelope 1D Squared Distance Transform (Felzenszwalb-Huttenlocher) +static inline void edt_1d_parabolic_row(const float* f, float* d, int* v, float* z, int n) { + int k = 0; + v[0] = 0; + z[0] = -1e20f; + z[1] = 1e20f; + + for (int q = 1; q < n; ++q) { + float f_q = f[q]; + if (f_q > 1e10f) continue; + float s = ((f_q + (float)q * q) - (f[v[k]] + (float)v[k] * v[k])) / (2.0f * q - 2.0f * v[k]); + while (s <= z[k]) { + k--; + s = ((f_q + (float)q * q) - (f[v[k]] + (float)v[k] * v[k])) / (2.0f * q - 2.0f * v[k]); + } + k++; + v[k] = q; + z[k] = s; + z[k + 1] = 1e20f; + } + + k = 0; + for (int q = 0; q < n; ++q) { + while (z[k + 1] < (float)q) { + k++; + } + float dx = (float)(q - v[k]); + d[q] = dx * dx + f[v[k]]; + } +} + +static inline void edt_2d_layer( + const uint8_t* mask, + int width, + int height, + float* out_dist_sq, + std::vector& g_col_buf, + std::vector& row_f, + std::vector& row_d, + std::vector& v_buf, + std::vector& z_buf +) { + // Pass 1: Columns + for (int x = 0; x < width; ++x) { + edt_1d_binary_col(mask + x, height, width, g_col_buf.data() + (size_t)x * height); + } + + // Pass 2: Rows + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + row_f[x] = g_col_buf[(size_t)x * height + y]; + } + edt_1d_parabolic_row(row_f.data(), row_d.data(), v_buf.data(), z_buf.data(), width); + for (int x = 0; x < width; ++x) { + out_dist_sq[(size_t)y * width + x] = row_d[x]; + } + } +} + +// ----------------------------------------------------------------------------- +// CPU Multi-Threaded Polygon Rasterizer +// ----------------------------------------------------------------------------- +struct CpuEdge2D { + float x1, y1; + float x2, y2; +}; + +void cpu_rasterize_polygon_layers_internal( + const GpuSliceLayer* in_layers, + size_t num_layers, + uint8_t* out_masks, + int width, int height, + float min_x, float min_y, float pitch_mm +) { + unsigned int num_threads = std::max(1u, std::thread::hardware_concurrency()); + std::vector workers; + + auto worker_fn = [&](size_t start_layer, size_t end_layer) { + for (size_t l = start_layer; l < end_layer; ++l) { + uint8_t* mask = out_masks + l * ((size_t)width * height); + const GpuSliceLayer& layer = in_layers[l]; + + std::vector edges; + for (size_t c = 0; c < layer.num_contours; ++c) { + const GpuContour& cnt = layer.contours[c]; + if (cnt.num_points < 3) continue; + for (size_t p = 0; p < cnt.num_points; ++p) { + size_t p_next = (p + 1) % cnt.num_points; + edges.push_back({ cnt.points[p].x, cnt.points[p].y, + cnt.points[p_next].x, cnt.points[p_next].y }); + } + } + + for (int y = 0; y < height; ++y) { + float py = min_y + (y + 0.5f) * pitch_mm; + std::vector x_inters; + + for (const auto& e : edges) { + if ((e.y1 <= py && e.y2 > py) || (e.y2 <= py && e.y1 > py)) { + float x_int = e.x1 + (py - e.y1) * (e.x2 - e.x1) / (e.y2 - e.y1); + x_inters.push_back(x_int); + } + } + + std::sort(x_inters.begin(), x_inters.end()); + + for (size_t i = 0; i + 1 < x_inters.size(); i += 2) { + int x_start = std::max(0, (int)std::ceil((x_inters[i] - min_x) / pitch_mm - 0.5f)); + int x_end = std::min(width, (int)std::floor((x_inters[i + 1] - min_x) / pitch_mm - 0.5f) + 1); + + for (int x = x_start; x < x_end; ++x) { + mask[(size_t)y * width + x] = 1; + } + } + } + } + }; + + 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(worker_fn, s, e); + } + } + for (auto& w : workers) w.join(); +} + +// ----------------------------------------------------------------------------- +// CPU Implementation: Core Mask Offsetting +// ----------------------------------------------------------------------------- +int cpu_offset_compute_masks( + const uint8_t* in_masks, + int width, + int height, + int num_layers, + float pitch_mm, + const float* offset_distances, + int num_offsets, + uint8_t* out_masks, + GpuOffsetStats* stats +) { + if (!in_masks || !offset_distances || !out_masks || width <= 0 || height <= 0 || num_layers <= 0 || num_offsets <= 0) { + return GPU_OFFSET_ERR_INVALID_ARG; + } + + auto t_start = std::chrono::high_resolution_clock::now(); + size_t layer_pixels = (size_t)width * height; + size_t total_volume = (size_t)num_layers * layer_pixels; + + unsigned int num_threads = std::max(1u, std::thread::hardware_concurrency()); + std::vector workers; + + auto worker_fn = [&](size_t start_layer, size_t end_layer) { + std::vector dist_sq(layer_pixels); + std::vector g_col_buf(layer_pixels); + std::vector row_f(width), row_d(width); + std::vector v_buf(width + 1); + std::vector z_buf(width + 2); + + for (size_t l = start_layer; l < end_layer; ++l) { + const uint8_t* mask_layer = in_masks + l * layer_pixels; + + // 1. Compute 2D EDT + edt_2d_layer(mask_layer, width, height, dist_sq.data(), + g_col_buf, row_f, row_d, v_buf, z_buf); + + // 2. Thresholding for each offset distance + for (int o = 0; o < num_offsets; ++o) { + float thresh = offset_distances[o]; + float thresh_sq_px = (thresh / pitch_mm) * (thresh / pitch_mm); + uint8_t* out_layer = out_masks + (size_t)o * total_volume + l * layer_pixels; + + for (size_t i = 0; i < layer_pixels; ++i) { + out_layer[i] = (mask_layer[i] > 0 && dist_sq[i] >= thresh_sq_px) ? 1 : 0; + } + } + } + }; + + 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(worker_fn, s, e); + } + for (auto& w : workers) w.join(); + + auto t_end = std::chrono::high_resolution_clock::now(); + double total_ms = std::chrono::duration(t_end - t_start).count(); + + if (stats) { + memset(stats, 0, sizeof(GpuOffsetStats)); + stats->t_gpu_jfa_ms = total_ms; + stats->t_total_pipeline_ms = total_ms; + } + + return GPU_OFFSET_SUCCESS; +} + +// ----------------------------------------------------------------------------- +// CPU Implementation: Polygon Layers Offsetting +// ----------------------------------------------------------------------------- +int cpu_offset_polygon_layers( + const GpuSliceLayer* in_layers, + size_t num_layers, + float min_x, float min_y, float max_x, float max_y, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + GpuOffsetStats* stats +) { + if (!in_layers || num_layers == 0 || !config || !out_results) return GPU_OFFSET_ERR_INVALID_ARG; + + 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); + size_t layer_pixels = (size_t)cols * rows; + size_t total_voxels = (size_t)num_layers * layer_pixels; + + // Resolve Offsets (perimeters + infill) + 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; // -0.5W, -1.5W, ... + } + resolved_offsets[P] = (float)P * W; // -P * W for infill + } + + // 1. Multi-Threaded CPU Rasterization + auto t_rast_0 = std::chrono::high_resolution_clock::now(); + std::vector h_masks(total_voxels, 0); + cpu_rasterize_polygon_layers_internal(in_layers, num_layers, h_masks.data(), + cols, rows, bounds_min_x, bounds_min_y, pitch); + auto t_rast_1 = std::chrono::high_resolution_clock::now(); + double t_raster_ms = std::chrono::duration(t_rast_1 - t_rast_0).count(); + + // 2. Multi-Threaded CPU EDT + Thresholding + Vectorization + 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; + + int row_stride_words = (cols + 31) / 32; + + auto worker_fn = [&](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); + + 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); + + for (size_t l = start_layer; l < end_layer; ++l) { + const uint8_t* mask_ptr = h_masks.data() + l * layer_pixels; + + // Compute exact 2D EDT + edt_2d_layer(mask_ptr, cols, rows, dist_sq.data(), + g_col_buf, row_f, row_d, v_buf, z_buf); + + res_array[l].z_height = in_layers[l].z_height; + res_array[l].num_perimeters = P; + res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter)); + + // Perimeter passes + for (int o = 0; o < P; ++o) { + float thresh = resolved_offsets[o]; + float thresh_sq_px = (thresh / pitch) * (thresh / pitch); + + // Bit-pack thresholded mask + for (int y = 0; y < rows; ++y) { + const uint8_t* row_m = mask_ptr + (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; + } + } + + // Marching squares vectorization + 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); + + 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 pass (index P) + 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 = mask_ptr + (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; + } + } + + 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); + + 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 = 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); + } + } + } + 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(worker_fn, s, e); + } + for (auto& w : workers) w.join(); + + 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) { + memset(stats, 0, sizeof(GpuOffsetStats)); + stats->t_raster_ms = t_raster_ms; + stats->t_gpu_jfa_ms = t_total_ms - t_raster_ms; // EDT + Vectorize + stats->t_total_pipeline_ms = t_total_ms; + } + + *out_results = res_array; + return GPU_OFFSET_SUCCESS; +} + +// ----------------------------------------------------------------------------- +// CPU Implementation: Direct 3D Mesh Slicing & Offsetting +// ----------------------------------------------------------------------------- +int cpu_offset_slice_mesh_and_offset( + 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, + float first_layer_height, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + size_t* out_num_layers, + GpuOffsetStats* stats +) { + if (!triangles || num_triangles == 0 || !config || !out_results) return GPU_OFFSET_ERR_INVALID_ARG; + + 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; + } + + if (num_layers <= 0) return GPU_OFFSET_ERR_INVALID_ARG; + if (out_num_layers) *out_num_layers = num_layers; + + // Resolve Offsets (perimeters + infill) + 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; // -0.5W, -1.5W, ... + } + resolved_offsets[P] = (float)P * W; // -P * W for infill + } + + int row_stride_words = (cols + 31) / 32; + size_t layer_pixels = (size_t)cols * rows; + + // 1. Spatial Triangle Indexing into layers + auto t_idx_0 = std::chrono::high_resolution_clock::now(); + std::vector> layer_tri_indices(num_layers); + for (uint32_t i = 0; i < (uint32_t)num_triangles; ++i) { + const GpuTriangle3D& tri = triangles[i]; + float tz_min = std::min({ tri.v0[2], tri.v1[2], tri.v2[2] }); + float tz_max = std::max({ tri.v0[2], tri.v1[2], tri.v2[2] }); + + int l_start = (tz_min <= z_first) ? 0 : (int)std::ceil((tz_min - z_first) / layer_height); + int l_end = (int)std::floor((tz_max - z_first) / layer_height); + + l_start = std::max(0, l_start); + l_end = std::min(num_layers - 1, l_end); + + for (int l = l_start; l <= l_end; ++l) { + layer_tri_indices[l].push_back(i); + } + } + 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 + 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 worker_fn = [&](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); + + 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); + + 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); + + // Slicing & Parity Rasterization for Layer l + const auto& tri_indices = layer_tri_indices[l]; + for (uint32_t t_idx : tri_indices) { + const GpuTriangle3D& tri = triangles[t_idx]; + const float* v[3] = { tri.v0, tri.v1, tri.v2 }; + + float pts_x[3]; + float pts_y[3]; + int count = 0; + + 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 (std::abs(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; + + 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; + + 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); + + if (y_start < 0) y_start = 0; + if (y_end >= rows) y_end = rows - 1; + + 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); + 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; + parity_words[(size_t)y * row_stride_words + w] ^= (1u << b); + } + } + + // Scan Parity Words into Binary Mask + for (int y = 0; y < rows; ++y) { + const unsigned int* row_p = parity_words.data() + (size_t)y * row_stride_words; + uint8_t* row_m = layer_mask.data() + (size_t)y * cols; + uint8_t inside = 0; + int x_out = 0; + for (int w = 0; w < row_stride_words && x_out < cols; ++w) { + unsigned int word = row_p[w]; + for (int b = 0; b < 32 && x_out < cols; ++b, ++x_out) { + inside ^= (uint8_t)((word >> b) & 1u); + row_m[x_out] = inside; + } + } + } + + // 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); + + res_array[l].z_height = z0; + res_array[l].num_perimeters = P; + res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter)); + + // Perimeters + 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; + } + } + + 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); + + 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 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; + } + } + + 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); + + 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 = 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); + } + } + } + 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(worker_fn, s, e); + } + for (auto& w : workers) w.join(); + + 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) { + memset(stats, 0, sizeof(GpuOffsetStats)); + stats->t_mesh_slice_ms = t_idx_ms; + stats->t_raster_ms = t_idx_ms; + stats->t_total_pipeline_ms = t_total_ms; + } + + *out_results = res_array; + return GPU_OFFSET_SUCCESS; +} diff --git a/src/cpu_offset.h b/src/cpu_offset.h new file mode 100644 index 0000000..4914b70 --- /dev/null +++ b/src/cpu_offset.h @@ -0,0 +1,67 @@ +#ifndef CPU_OFFSET_H +#define CPU_OFFSET_H + +#include "gpu_offset.h" +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif + +void cpu_rasterize_polygon_layers_internal( + const GpuSliceLayer* in_layers, + size_t num_layers, + uint8_t* out_masks, + int width, int height, + float min_x, float min_y, float pitch_mm +); + +/** + * @brief CPU fallback for core Euclidean distance transform and offsetting on raw masks. + */ +int cpu_offset_compute_masks( + const uint8_t* in_masks, + int width, + int height, + int num_layers, + float pitch_mm, + const float* offset_distances, + int num_offsets, + uint8_t* out_masks, + GpuOffsetStats* stats +); + +/** + * @brief CPU fallback for 2D polygon layers offsetting. + */ +int cpu_offset_polygon_layers( + const GpuSliceLayer* in_layers, + size_t num_layers, + float min_x, float min_y, float max_x, float max_y, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + GpuOffsetStats* stats +); + +/** + * @brief CPU fallback for direct 3D mesh slicing, rasterization, distance transform, and offsetting. + */ +int cpu_offset_slice_mesh_and_offset( + 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, + float first_layer_height, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + size_t* out_num_layers, + GpuOffsetStats* stats +); + +#ifdef __cplusplus +} +#endif + +#endif // CPU_OFFSET_H diff --git a/src/gpu_offset.cu b/src/gpu_offset.cu new file mode 100644 index 0000000..e002878 --- /dev/null +++ b/src/gpu_offset.cu @@ -0,0 +1,1059 @@ +#include "gpu_offset.h" +#include "cpu_offset.h" +#include "vectorizer.h" +#include +#include +#include +#include +#include +#include +#include +#include + +// CUDA Error Checking Macro +#define CUDA_CHECK(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 device helper functions +__device__ __inline__ int pack_seed(short x, short y) { + return ((int)(unsigned short)y << 16) | (int)(unsigned short)x; +} +__device__ __inline__ short get_seed_x(int packed) { + return (short)(packed & 0xFFFF); +} +__device__ __inline__ short get_seed_y(int packed) { + return (short)((packed >> 16) & 0xFFFF); +} + +// ----------------------------------------------------------------------------- +// CUDA Kernels: Batched 3D Jump Flooding Algorithm +// ----------------------------------------------------------------------------- +__global__ void batched_jfa_init_kernel( + const uint8_t* __restrict__ masks, + int* __restrict__ seeds, + 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 layer_stride = (size_t)width * height; + size_t idx = (size_t)z * layer_stride + y * width + x; + + // Background pixel (mask == 0) is initialized as seed point + if (masks[idx] == 0) { + seeds[idx] = pack_seed((short)x, (short)y); + } else { + seeds[idx] = -1; + } +} + +__global__ void batched_jfa_step_kernel( + const int* __restrict__ in_seeds, + int* __restrict__ out_seeds, + 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_stride = (size_t)width * height; + size_t layer_offset = (size_t)z * layer_stride; + size_t idx = layer_offset + y * width + x; + + int best_seed = in_seeds[idx]; + int best_dist2 = 0x7FFFFFFF; + if (best_seed != -1) { + int dx = (int)get_seed_x(best_seed) - x; + int dy = (int)get_seed_y(best_seed) - y; + best_dist2 = 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) { + int n_seed = in_seeds[layer_offset + ny * width + nx]; + if (n_seed != -1) { + int dist_x = (int)get_seed_x(n_seed) - x; + int dist_y = (int)get_seed_y(n_seed) - y; + int dist2 = dist_x * dist_x + dist_y * dist_y; + if (dist2 < best_dist2) { + best_dist2 = dist2; + best_seed = n_seed; + } + } + } + } + } + out_seeds[idx] = best_seed; +} + +__global__ void batched_jfa_threshold_kernel( + const int* __restrict__ seeds, + const uint8_t* __restrict__ masks, + uint8_t* __restrict__ out_offsets, + int width, int height, int layers, + float pitch_mm, + const float* __restrict__ thresholds_mm, + int num_offsets +) { + 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_pixels = (size_t)width * height; + size_t total_volume = (size_t)layers * layer_pixels; + size_t in_idx = (size_t)z * layer_pixels + y * width + x; + + uint8_t mask_val = masks[in_idx]; + int s = seeds[in_idx]; + + float dist_mm = 0.0f; + if (mask_val > 0 && s != -1) { + float dx = (float)(get_seed_x(s) - x); + float dy = (float)(get_seed_y(s) - y); + dist_mm = sqrtf(dx * dx + dy * dy) * pitch_mm; + } + + for (int o = 0; o < num_offsets; ++o) { + size_t out_idx = (size_t)o * total_volume + in_idx; + out_offsets[out_idx] = (dist_mm >= thresholds_mm[o]) ? 1 : 0; + } +} + +__global__ void batched_jfa_threshold_bitpack_kernel( + const int* __restrict__ seeds, + const uint8_t* __restrict__ masks, + unsigned int* __restrict__ out_bit_words, + int width, int height, int layers, + float pitch_mm, + 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 layer_pixels = (size_t)width * height; + size_t in_layer_offset = (size_t)z * layer_pixels + 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]; + float thresh_sq_px = (thresh / pitch_mm) * (thresh / pitch_mm); + + unsigned int packed = 0; + + #pragma unroll + for (int b = 0; b < 32; ++b) { + int x = x_base + b; + if (x < width) { + size_t in_idx = in_layer_offset + x; + uint8_t m = masks[in_idx]; + int s = seeds[in_idx]; + + if (m > 0 && s != -1) { + int sx = (int)get_seed_x(s); + int sy = (int)get_seed_y(s); + int dx = sx - x; + int dy = sy - y; + int dist_sq = dx * dx + dy * dy; + if ((float)dist_sq >= thresh_sq_px) { + packed |= (1u << b); + } + } + } + } + + size_t out_idx = (size_t)o * total_word_layer + out_layer_offset + x_word; + out_bit_words[out_idx] = packed; + } +} + +// ----------------------------------------------------------------------------- +// Direct GPU Mesh Slicing & Scanline Parity Kernels +// ----------------------------------------------------------------------------- +__global__ void gpu_mesh_slice_raster_kernel( + const GpuTriangle3D* __restrict__ triangles, + size_t num_triangles, + 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 +) { + 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; + + float gx0 = (pts_x[0] - bounds_min_x) / pitch_mm - 0.5f; + float gy0 = (pts_y[0] - bounds_min_y) / pitch_mm - 0.5f; + float gx1 = (pts_x[1] - bounds_min_x) / pitch_mm - 0.5f; + float gy1 = (pts_y[1] - bounds_min_y) / pitch_mm - 0.5f; + + if (gy0 > gy1) { + float tx = gx0; gx0 = gx1; gx1 = tx; + float ty = gy0; gy0 = gy1; gy1 = ty; + } + + float dy = gy1 - gy0; + if (dy < 1e-6f) continue; + + int y_start = (int)ceilf(gy0); + int y_end = (int)floorf(gy1 - 1e-6f); + + if (y_start < 0) y_start = 0; + if (y_end >= rows) y_end = rows - 1; + + 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)roundf(gx + 0.5f); + 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_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; + } + } +} + +// ----------------------------------------------------------------------------- +// Public C API Implementation & Backend Management +// ----------------------------------------------------------------------------- +extern "C" { + +static cudaStream_t g_cuda_stream = nullptr; +static int g_initialized = 0; +static int g_global_backend = GPU_OFFSET_BACKEND_AUTO; + +GPU_OFFSET_API int gpu_offset_has_cuda(void) { + int count = 0; + cudaError_t err = cudaGetDeviceCount(&count); + if (err != cudaSuccess || count <= 0) { + cudaGetLastError(); // Clear error + return 0; + } + return 1; +} + +GPU_OFFSET_API int gpu_offset_set_backend(int backend) { + if (backend < 0 || backend > 2) return GPU_OFFSET_ERR_INVALID_ARG; + g_global_backend = backend; + return GPU_OFFSET_SUCCESS; +} + +GPU_OFFSET_API int gpu_offset_get_backend(void) { + return g_global_backend; +} + +GPU_OFFSET_API int gpu_offset_init(int device_id) { + if (g_initialized) return GPU_OFFSET_SUCCESS; + + if (!gpu_offset_has_cuda()) { + if (g_global_backend == GPU_OFFSET_BACKEND_CUDA) { + return GPU_OFFSET_ERR_NO_CUDA; + } + return GPU_OFFSET_SUCCESS; // CPU fallback ready + } + + CUDA_CHECK(cudaSetDevice(device_id)); + CUDA_CHECK(cudaStreamCreate(&g_cuda_stream)); + + // Warm up CUDA runtime + void* d_dummy; + CUDA_CHECK(cudaMalloc(&d_dummy, 1024)); + CUDA_CHECK(cudaFree(d_dummy)); + + g_initialized = 1; + return GPU_OFFSET_SUCCESS; +} + +GPU_OFFSET_API void gpu_offset_cleanup() { + if (g_cuda_stream) { + cudaStreamDestroy(g_cuda_stream); + g_cuda_stream = nullptr; + } + g_initialized = 0; +} + +// Core GPU Distance Transform execution on masks +static int gpu_compute_masks_internal( + const uint8_t* h_masks, + int width, + int height, + int num_layers, + float pitch_mm, + const float* offset_distances, + int num_offsets, + uint8_t* h_out_masks, + GpuOffsetStats* stats +) { + if (!g_initialized) { + if (gpu_offset_init(0) != GPU_OFFSET_SUCCESS) return GPU_OFFSET_ERR_GENERIC; + } + + size_t layer_pixels = (size_t)width * height; + size_t total_volume = (size_t)num_layers * layer_pixels; + size_t total_out_bytes = (size_t)num_offsets * total_volume * sizeof(uint8_t); + + cudaEvent_t ev_start, ev_h2d_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_jfa_done)); + CUDA_CHECK(cudaEventCreate(&ev_thresh_done)); + CUDA_CHECK(cudaEventCreate(&ev_d2h_done)); + + // Allocate GPU buffers + uint8_t* d_masks; + int* d_seeds_a; + int* d_seeds_b; + uint8_t* d_out_offsets; + float* d_thresholds; + + 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_offsets, total_out_bytes)); + CUDA_CHECK(cudaMalloc(&d_thresholds, num_offsets * sizeof(float))); + + // 1. H2D Transfer + CUDA_CHECK(cudaEventRecord(ev_start, g_cuda_stream)); + CUDA_CHECK(cudaMemcpyAsync(d_masks, h_masks, total_volume * sizeof(uint8_t), + cudaMemcpyHostToDevice, g_cuda_stream)); + CUDA_CHECK(cudaMemcpyAsync(d_thresholds, offset_distances, num_offsets * sizeof(float), + cudaMemcpyHostToDevice, g_cuda_stream)); + CUDA_CHECK(cudaEventRecord(ev_h2d_done, g_cuda_stream)); + + // 2. Batched 3D JFA Computation + dim3 block(16, 16, 1); + dim3 grid((width + 15) / 16, (height + 15) / 16, num_layers); + + batched_jfa_init_kernel<<>>( + d_masks, d_seeds_a, width, height, num_layers + ); + + int max_dim = std::max(width, height); + int step = 1; + while (step * 2 < max_dim) step *= 2; + + int* in_seeds = d_seeds_a; + int* out_seeds = d_seeds_b; + + while (step >= 1) { + batched_jfa_step_kernel<<>>( + in_seeds, out_seeds, width, height, num_layers, step + ); + std::swap(in_seeds, out_seeds); + step /= 2; + } + + // JFA+2 and JFA+1 correction passes + batched_jfa_step_kernel<<>>( + in_seeds, out_seeds, width, height, num_layers, 2 + ); + std::swap(in_seeds, out_seeds); + + batched_jfa_step_kernel<<>>( + in_seeds, out_seeds, width, height, num_layers, 1 + ); + std::swap(in_seeds, out_seeds); + + CUDA_CHECK(cudaEventRecord(ev_jfa_done, g_cuda_stream)); + + // 3. Multi-Offset Thresholding + batched_jfa_threshold_kernel<<>>( + in_seeds, d_masks, d_out_offsets, + width, height, num_layers, + pitch_mm, d_thresholds, num_offsets + ); + CUDA_CHECK(cudaEventRecord(ev_thresh_done, g_cuda_stream)); + + // 4. D2H Transfer + CUDA_CHECK(cudaMemcpyAsync(h_out_masks, d_out_offsets, total_out_bytes, + cudaMemcpyDeviceToHost, g_cuda_stream)); + CUDA_CHECK(cudaEventRecord(ev_d2h_done, g_cuda_stream)); + CUDA_CHECK(cudaStreamSynchronize(g_cuda_stream)); + + // Timing stats + if (stats) { + float h2d_ms = 0, jfa_ms = 0, thresh_ms = 0, d2h_ms = 0; + cudaEventElapsedTime(&h2d_ms, ev_start, ev_h2d_done); + cudaEventElapsedTime(&jfa_ms, ev_h2d_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_gpu_jfa_ms = jfa_ms; + stats->t_gpu_thresh_ms = thresh_ms; + stats->t_d2h_ms = d2h_ms; + } + + // Cleanup + cudaFree(d_masks); + cudaFree(d_seeds_a); + cudaFree(d_seeds_b); + cudaFree(d_out_offsets); + cudaFree(d_thresholds); + + cudaEventDestroy(ev_start); + cudaEventDestroy(ev_h2d_done); + cudaEventDestroy(ev_jfa_done); + cudaEventDestroy(ev_thresh_done); + cudaEventDestroy(ev_d2h_done); + + return GPU_OFFSET_SUCCESS; +} + +GPU_OFFSET_API int gpu_offset_compute_masks( + const uint8_t* h_masks, + int width, + int height, + int num_layers, + float pitch_mm, + const float* offset_distances, + int num_offsets, + uint8_t* h_out_masks, + GpuOffsetStats* stats +) { + int eff_backend = g_global_backend; + if (eff_backend == GPU_OFFSET_BACKEND_AUTO) { + eff_backend = gpu_offset_has_cuda() ? GPU_OFFSET_BACKEND_CUDA : GPU_OFFSET_BACKEND_CPU; + } + + if (eff_backend == GPU_OFFSET_BACKEND_CUDA) { + if (!gpu_offset_has_cuda()) { + return GPU_OFFSET_ERR_NO_CUDA; + } + return gpu_compute_masks_internal(h_masks, width, height, num_layers, pitch_mm, + offset_distances, num_offsets, h_out_masks, stats); + } else { + return cpu_offset_compute_masks(h_masks, width, height, num_layers, pitch_mm, + offset_distances, num_offsets, h_out_masks, stats); + } +} + +GPU_OFFSET_API int gpu_offset_polygon_layers( + const GpuSliceLayer* in_layers, + size_t num_layers, + float min_x, float min_y, float max_x, float max_y, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + GpuOffsetStats* stats +) { + if (!in_layers || num_layers == 0 || !config || !out_results) return GPU_OFFSET_ERR_INVALID_ARG; + + int eff_backend = (config->backend != GPU_OFFSET_BACKEND_AUTO) ? config->backend : g_global_backend; + if (eff_backend == GPU_OFFSET_BACKEND_AUTO) { + eff_backend = gpu_offset_has_cuda() ? GPU_OFFSET_BACKEND_CUDA : GPU_OFFSET_BACKEND_CPU; + } + + if (eff_backend == GPU_OFFSET_BACKEND_CUDA) { + if (!gpu_offset_has_cuda()) { + return GPU_OFFSET_ERR_NO_CUDA; + } + + // GPU Offsetting pipeline + 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); + size_t total_voxels = (size_t)num_layers * cols * rows; + + 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(); + + // 1. Rasterization + auto t_rast_0 = std::chrono::high_resolution_clock::now(); + std::vector h_masks(total_voxels, 0); + + // Uses CPU polygon rasterizer to populate h_masks + extern void cpu_rasterize_polygon_layers_internal( + const GpuSliceLayer* in_layers, size_t num_layers, uint8_t* out_masks, + int width, int height, float min_x, float min_y, float pitch_mm); + + cpu_rasterize_polygon_layers_internal(in_layers, num_layers, h_masks.data(), + cols, rows, bounds_min_x, bounds_min_y, pitch); + auto t_rast_1 = std::chrono::high_resolution_clock::now(); + double t_raster_ms = std::chrono::duration(t_rast_1 - t_rast_0).count(); + + // 2. GPU Distance Transform & Offsetting + std::vector h_offset_masks(num_offsets * total_voxels); + if (gpu_compute_masks_internal( + h_masks.data(), cols, rows, (int)num_layers, pitch, + resolved_offsets.data(), (int)num_offsets, + h_offset_masks.data(), stats) != GPU_OFFSET_SUCCESS) { + return GPU_OFFSET_ERR_GENERIC; + } + + if (stats) stats->t_raster_ms = t_raster_ms; + + // 3. Multi-Threaded C++ 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 = in_layers[l].z_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_idx = (o * num_layers + l) * ((size_t)cols * rows); + const uint8_t* mask_ptr = h_offset_masks.data() + offset_layer_idx; + + std::vector contours; + vectorize_mask_layer(mask_ptr, cols, rows, 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_idx = (P * num_layers + l) * ((size_t)cols * rows); + const uint8_t* infill_mask_ptr = h_offset_masks.data() + infill_offset_layer_idx; + + std::vector infill_contours; + vectorize_mask_layer(infill_mask_ptr, cols, rows, 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); + memcpy(infill_mask, infill_mask_ptr, (size_t)cols * rows); + 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; + } + + *out_results = res_array; + return GPU_OFFSET_SUCCESS; + } else { + return cpu_offset_polygon_layers(in_layers, num_layers, min_x, min_y, max_x, max_y, config, out_results, stats); + } +} + +GPU_OFFSET_API int gpu_offset_slice_mesh_and_offset( + 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, + float first_layer_height, + const GpuOffsetConfig* config, + GpuLayerResult** out_results, + size_t* out_num_layers, + GpuOffsetStats* stats +) { + if (!triangles || num_triangles == 0 || !config || !out_results) return GPU_OFFSET_ERR_INVALID_ARG; + + int eff_backend = (config->backend != GPU_OFFSET_BACKEND_AUTO) ? config->backend : g_global_backend; + if (eff_backend == GPU_OFFSET_BACKEND_AUTO) { + eff_backend = gpu_offset_has_cuda() ? GPU_OFFSET_BACKEND_CUDA : GPU_OFFSET_BACKEND_CPU; + } + + if (eff_backend == GPU_OFFSET_BACKEND_CUDA) { + if (!gpu_offset_has_cuda()) { + return GPU_OFFSET_ERR_NO_CUDA; + } + + if (!g_initialized) { + 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; + } + + 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; + } 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); + } +} + +GPU_OFFSET_API void gpu_offset_free_results(GpuLayerResult* results, size_t num_layers) { + if (!results) return; + for (size_t l = 0; l < num_layers; ++l) { + if (results[l].perimeters) { + for (size_t o = 0; o < results[l].num_perimeters; ++o) { + if (results[l].perimeters[o].contours) { + for (size_t c = 0; c < results[l].perimeters[o].num_contours; ++c) { + if (results[l].perimeters[o].contours[c].points) { + free(results[l].perimeters[o].contours[c].points); + } + } + free(results[l].perimeters[o].contours); + } + } + free(results[l].perimeters); + } + if (results[l].infill.contours) { + for (size_t c = 0; c < results[l].infill.num_contours; ++c) { + if (results[l].infill.contours[c].points) { + free(results[l].infill.contours[c].points); + } + } + free(results[l].infill.contours); + } + if (results[l].infill.mask) { + free(results[l].infill.mask); + } + } + free(results); +} + +} // extern "C" diff --git a/src/main_bench.cpp b/src/main_bench.cpp new file mode 100644 index 0000000..792c2d5 --- /dev/null +++ b/src/main_bench.cpp @@ -0,0 +1,292 @@ +#include "gpu_offset.h" +#include "stl_slicer.h" +#include +#include +#include +#include +#include +#include +#include + +int main(int argc, char** argv) { + std::string stl_path = "benchy.stl"; + float layer_height = 0.20f; + float first_layer_height = 0.20f; + float pitch_mm = 0.10f; + int perimeters_count = 2; + float layer_width = 0.45f; + int runs = 3; + bool run_prusa = true; + + for (int i = 1; i < argc; ++i) { + std::string arg = argv[i]; + if (arg == "--stl" && i + 1 < argc) stl_path = argv[++i]; + else if (arg == "--layer-height" && i + 1 < argc) layer_height = std::stof(argv[++i]); + else if (arg == "--first-layer-height" && i + 1 < argc) first_layer_height = std::stof(argv[++i]); + else if (arg == "--perimeters" && i + 1 < argc) perimeters_count = std::stoi(argv[++i]); + else if (arg == "--layer-width" && i + 1 < argc) layer_width = std::stof(argv[++i]); + else if (arg == "--pitch" && i + 1 < argc) pitch_mm = std::stof(argv[++i]); + else if (arg == "--runs" && i + 1 < argc) runs = std::stoi(argv[++i]); + else if (arg == "--no-prusa") run_prusa = false; + } + + std::cout << "================================================================================" << std::endl; + std::cout << " HIGH-PERFORMANCE MESH SLICING & OFFSETTING BENCHMARK (CUDA & CPU)" << std::endl; + std::cout << "================================================================================" << std::endl; + + int has_cuda = gpu_offset_has_cuda(); + std::cout << "[Device Capability] CUDA Available: " << (has_cuda ? "YES (GPU Accelerated)" : "NO (CPU Fallback)") << std::endl; + + // 1. Initialize Engine + if (gpu_offset_init(0) != GPU_OFFSET_SUCCESS) { + std::cerr << "[Fatal Error] Failed to initialize backend" << std::endl; + return 1; + } + + // 2. Load STL Mesh + std::cout << "\n[Step 1] Loading Binary STL Mesh: " << stl_path << std::endl; + auto t_load_0 = std::chrono::high_resolution_clock::now(); + StlMesh mesh; + if (!load_stl_binary(stl_path, mesh)) { + return 1; + } + auto t_load_1 = std::chrono::high_resolution_clock::now(); + double t_load_ms = std::chrono::duration(t_load_1 - t_load_0).count(); + + std::cout << " -> Loaded " << mesh.triangles.size() << " triangles in " << t_load_ms << " ms" << std::endl; + std::cout << " -> Bounds: [" << mesh.min_bound.x << ", " << mesh.min_bound.y << ", " << mesh.min_bound.z + << "] to [" << mesh.max_bound.x << ", " << mesh.max_bound.y << ", " << mesh.max_bound.z << "]" << std::endl; + + // 3. Slice Mesh in C++ (Hybrid Baseline) + std::cout << "\n[Step 2] Slicing Mesh with Multi-Threaded CPU Engine at " << layer_height << " mm Layer Height (First Layer: " << first_layer_height << " mm)" << std::endl; + std::vector slice_layers; + std::vector z_levels; + double t_slice_cpu_ms = 0; + slice_stl_mesh(mesh, layer_height, first_layer_height, slice_layers, z_levels, t_slice_cpu_ms); + size_t num_layers = slice_layers.size(); + + std::cout << " -> Sliced " << num_layers << " layers in " << t_slice_cpu_ms << " ms (" + << (t_slice_cpu_ms / num_layers) << " ms/layer)" << std::endl; + + // 4. Configure Offsetting Engine + GpuOffsetConfig config_gpu = {0}; + config_gpu.pitch_mm = pitch_mm; + config_gpu.pad_mm = 1.5f; + config_gpu.perimeters_count = perimeters_count; + config_gpu.layer_width = layer_width; + config_gpu.first_layer_height = first_layer_height; + config_gpu.export_infill_mask = 1; + config_gpu.use_gpu_rasterizer = 1; + config_gpu.backend = GPU_OFFSET_BACKEND_CUDA; + + GpuOffsetConfig config_cpu = config_gpu; + config_cpu.backend = GPU_OFFSET_BACKEND_CPU; + + float bounds_min_x = mesh.min_bound.x - config_gpu.pad_mm; + float bounds_min_y = mesh.min_bound.y - config_gpu.pad_mm; + float bounds_max_x = mesh.max_bound.x + config_gpu.pad_mm; + float bounds_max_y = mesh.max_bound.y + config_gpu.pad_mm; + int cols = (int)std::ceil((bounds_max_x - bounds_min_x) / pitch_mm); + int rows = (int)std::ceil((bounds_max_y - bounds_min_y) / pitch_mm); + size_t total_voxels = (size_t)num_layers * cols * rows; + + std::cout << "\n[Step 3] Grid Configuration (Pitch: " << pitch_mm << " mm, Line Width: " << layer_width << " mm)" << std::endl; + std::cout << " -> Grid Size: " << cols << " x " << rows << " = " << (cols * rows) << " px/layer" << std::endl; + std::cout << " -> Total 3D Voxel Volume: " << total_voxels << " voxels (" + << (total_voxels / (1024.0 * 1024.0)) << " MB uint8)" << std::endl; + + // 5. Benchmark Hybrid Pipeline: CPU Slicing + GPU Offsetting (if CUDA available) + GpuOffsetStats avg_hybrid = {0}; + if (has_cuda) { + std::cout << "\n[Step 4] Benchmarking Pipeline A: Multi-Threaded CPU Slicing + GPU Offsetting (" << runs << " runs)..." << std::endl; + for (int r = 0; r < runs; ++r) { + GpuOffsetStats stats = {0}; + GpuLayerResult* results = nullptr; + + gpu_offset_polygon_layers( + slice_layers.data(), num_layers, + mesh.min_bound.x, mesh.min_bound.y, + mesh.max_bound.x, mesh.max_bound.y, + &config_gpu, &results, &stats + ); + + avg_hybrid.t_raster_ms += stats.t_raster_ms; + avg_hybrid.t_h2d_ms += stats.t_h2d_ms; + avg_hybrid.t_gpu_jfa_ms += stats.t_gpu_jfa_ms; + avg_hybrid.t_gpu_thresh_ms += stats.t_gpu_thresh_ms; + avg_hybrid.t_d2h_ms += stats.t_d2h_ms; + avg_hybrid.t_vectorize_ms += stats.t_vectorize_ms; + avg_hybrid.t_total_pipeline_ms += stats.t_total_pipeline_ms; + + gpu_offset_free_results(results, num_layers); + } + + avg_hybrid.t_raster_ms /= runs; + avg_hybrid.t_h2d_ms /= runs; + avg_hybrid.t_gpu_jfa_ms /= runs; + avg_hybrid.t_gpu_thresh_ms /= runs; + avg_hybrid.t_d2h_ms /= runs; + avg_hybrid.t_vectorize_ms /= runs; + avg_hybrid.t_total_pipeline_ms /= runs; + } + + // 6. Benchmark 100% GPU Slicing Pipeline + std::vector gpu_tris(mesh.triangles.size()); + for (size_t i = 0; i < mesh.triangles.size(); ++i) { + const auto& t = mesh.triangles[i]; + gpu_tris[i].v0[0] = t.v0.x; gpu_tris[i].v0[1] = t.v0.y; gpu_tris[i].v0[2] = t.v0.z; + gpu_tris[i].v1[0] = t.v1.x; gpu_tris[i].v1[1] = t.v1.y; gpu_tris[i].v1[2] = t.v1.z; + gpu_tris[i].v2[0] = t.v2.x; gpu_tris[i].v2[1] = t.v2.y; gpu_tris[i].v2[2] = t.v2.z; + } + + GpuOffsetStats avg_gpu_slice = {0}; + size_t gpu_num_layers = 0; + + if (has_cuda) { + std::cout << "\n[Step 5] Benchmarking Pipeline B: 100% Direct GPU Mesh Slicing & Offsetting (" << runs << " runs)..." << std::endl; + for (int r = 0; r < runs; ++r) { + GpuOffsetStats stats = {0}; + GpuLayerResult* results = nullptr; + + int ret = gpu_offset_slice_mesh_and_offset( + gpu_tris.data(), gpu_tris.size(), + mesh.min_bound.x, mesh.min_bound.y, mesh.min_bound.z, + mesh.max_bound.x, mesh.max_bound.y, mesh.max_bound.z, + layer_height, first_layer_height, &config_gpu, &results, &gpu_num_layers, &stats + ); + + if (ret != GPU_OFFSET_SUCCESS) { + std::cerr << "[Error] GPU direct slice failed" << std::endl; + return 1; + } + + avg_gpu_slice.t_h2d_ms += stats.t_h2d_ms; + avg_gpu_slice.t_mesh_slice_ms += stats.t_mesh_slice_ms; + avg_gpu_slice.t_gpu_jfa_ms += stats.t_gpu_jfa_ms; + avg_gpu_slice.t_gpu_thresh_ms += stats.t_gpu_thresh_ms; + avg_gpu_slice.t_d2h_ms += stats.t_d2h_ms; + avg_gpu_slice.t_vectorize_ms += stats.t_vectorize_ms; + avg_gpu_slice.t_total_pipeline_ms += stats.t_total_pipeline_ms; + + gpu_offset_free_results(results, gpu_num_layers); + } + + avg_gpu_slice.t_h2d_ms /= runs; + avg_gpu_slice.t_mesh_slice_ms /= runs; + avg_gpu_slice.t_gpu_jfa_ms /= runs; + avg_gpu_slice.t_gpu_thresh_ms /= runs; + avg_gpu_slice.t_d2h_ms /= runs; + avg_gpu_slice.t_vectorize_ms /= runs; + avg_gpu_slice.t_total_pipeline_ms /= runs; + } + + // 7. Benchmark 100% CPU Compatibility Pipeline + std::cout << "\n[Step 6] Benchmarking Pipeline C: 100% CPU Compatibility Direct Slicing & Offsetting (" << runs << " runs)..." << std::endl; + GpuOffsetStats avg_cpu_slice = {0}; + size_t cpu_num_layers = 0; + + for (int r = 0; r < runs; ++r) { + GpuOffsetStats stats = {0}; + GpuLayerResult* results = nullptr; + + int ret = gpu_offset_slice_mesh_and_offset( + gpu_tris.data(), gpu_tris.size(), + mesh.min_bound.x, mesh.min_bound.y, mesh.min_bound.z, + mesh.max_bound.x, mesh.max_bound.y, mesh.max_bound.z, + layer_height, first_layer_height, &config_cpu, &results, &cpu_num_layers, &stats + ); + + if (ret != GPU_OFFSET_SUCCESS) { + std::cerr << "[Error] CPU direct slice failed" << std::endl; + return 1; + } + + avg_cpu_slice.t_mesh_slice_ms += stats.t_mesh_slice_ms; + avg_cpu_slice.t_total_pipeline_ms += stats.t_total_pipeline_ms; + + gpu_offset_free_results(results, cpu_num_layers); + } + + avg_cpu_slice.t_mesh_slice_ms /= runs; + avg_cpu_slice.t_total_pipeline_ms /= runs; + + // 8. PrusaSlicer Benchmark + double t_prusa_ms = 0; + if (run_prusa) { + std::cout << "\n[Step 7] Running PrusaSlicer Headless Baseline Benchmark..." << std::endl; + std::string cmd = "./bin/prusa-slicer --slice --layer-height=0.2 --perimeters=2 --fill-density=0% --output=prusa_bench.gcode " + stl_path + " > /dev/null 2>&1"; + auto t_p0 = std::chrono::high_resolution_clock::now(); + int prusa_ret = std::system(cmd.c_str()); + auto t_p1 = std::chrono::high_resolution_clock::now(); + if (prusa_ret == 0) { + t_prusa_ms = std::chrono::duration(t_p1 - t_p0).count(); + std::cout << " -> PrusaSlicer completed in " << t_prusa_ms << " ms (" << (t_prusa_ms / 1000.0) << " s)" << std::endl; + } + } + + // 9. Summary Table + double total_gpu_slice_ms = t_load_ms + avg_gpu_slice.t_total_pipeline_ms; + double total_cpu_slice_ms = t_load_ms + avg_cpu_slice.t_total_pipeline_ms; + + std::cout << "\n" << std::string(88, '=') << std::endl; + std::cout << " BENCHMARK COMPARISON: 100% GPU vs 100% CPU FALLBACK vs PRUSASLICER" << std::endl; + std::cout << std::string(88, '=') << std::endl; + + std::cout << std::left << std::setw(54) << "Pipeline / Execution Stage" + << " | " << std::right << std::setw(11) << "Time (ms)" + << " | " << std::right << std::setw(14) << "Rate (layer/s)" << std::endl; + std::cout << std::string(88, '-') << std::endl; + + auto print_row = [&](const std::string& name, double ms) { + double rate = ms > 0 ? (num_layers / (ms / 1000.0)) : 0.0; + std::cout << std::left << std::setw(54) << name + << " | " << std::right << std::setw(8) << std::fixed << std::setprecision(2) << ms << " ms" + << " | " << std::right << std::setw(14) << std::fixed << std::setprecision(1) << rate << std::endl; + }; + + if (has_cuda) { + std::cout << ">> PIPELINE B: 100% DIRECT GPU MESH SLICING & OFFSETTING <<" << std::endl; + print_row(" 0. STL Mesh Load (Binary C++ Parser)", t_load_ms); + print_row(" 1. Triangle Vertex Buffer PCIe Upload (H2D)", avg_gpu_slice.t_h2d_ms); + print_row(" 2. GPU Direct Mesh Slicing & Rasterization", avg_gpu_slice.t_mesh_slice_ms); + print_row(" 3. GPU Batched 3D CUDA JFA Distance Transform", avg_gpu_slice.t_gpu_jfa_ms); + print_row(" 4. GPU Multi-Perimeter Kernel Thresholding", avg_gpu_slice.t_gpu_thresh_ms); + print_row(" 5. Device-to-Host (D2H) Offset Mask Download", avg_gpu_slice.t_d2h_ms); + print_row(" 6. Multi-Threaded C++ Vectorization (Marching Sq)", avg_gpu_slice.t_vectorize_ms); + std::cout << std::string(88, '-') << std::endl; + print_row("TOTAL 100% GPU PIPELINE (STL -> Sliced Contours)", total_gpu_slice_ms); + std::cout << std::string(88, '-') << std::endl; + } + + std::cout << ">> PIPELINE C: 100% MULTI-THREADED CPU FALLBACK PIPELINE <<" << std::endl; + print_row(" 0. STL Mesh Load (Binary C++ Parser)", t_load_ms); + print_row(" 1. Spatial Indexing + CPU Direct Slicing & Raster", avg_cpu_slice.t_mesh_slice_ms); + print_row(" 2. CPU 2D Fast EDT + Thresholding + Vectorization", avg_cpu_slice.t_total_pipeline_ms - avg_cpu_slice.t_mesh_slice_ms); + std::cout << std::string(88, '-') << std::endl; + print_row("TOTAL 100% CPU PIPELINE (STL -> Sliced Contours)", total_cpu_slice_ms); + std::cout << std::string(88, '-') << std::endl; + + if (t_prusa_ms > 0) { + std::cout << ">> BASELINE: PRUSASLICER HEADLESS CLI <<" << std::endl; + print_row("PRUSASLICER CLI (Mesh -> Perimeters & G-Code)", t_prusa_ms); + std::cout << std::string(88, '=') << std::endl; + std::cout << "\n>>> SPEEDUP SUMMARY:" << std::endl; + if (has_cuda) { + std::cout << " * 100% GPU Slicing Pipeline : " << std::fixed << std::setprecision(2) + << (t_prusa_ms / total_gpu_slice_ms) << "x FASTER than PrusaSlicer (" + << std::fixed << std::setprecision(1) << (num_layers / (total_gpu_slice_ms / 1000.0)) + << " vs " << (num_layers / (t_prusa_ms / 1000.0)) << " layers/sec)" << std::endl; + } + std::cout << " * 100% CPU Fallback Pipeline: " << std::fixed << std::setprecision(2) + << (t_prusa_ms / total_cpu_slice_ms) << "x FASTER than PrusaSlicer (" + << std::fixed << std::setprecision(1) << (num_layers / (total_cpu_slice_ms / 1000.0)) + << " vs " << (num_layers / (t_prusa_ms / 1000.0)) << " layers/sec)" << std::endl; + } + std::cout << std::string(88, '=') << std::endl; + + // Cleanup + free_slice_layers(slice_layers); + gpu_offset_cleanup(); + + return 0; +} diff --git a/src/stl_slicer.cpp b/src/stl_slicer.cpp new file mode 100644 index 0000000..526d3c2 --- /dev/null +++ b/src/stl_slicer.cpp @@ -0,0 +1,273 @@ +#include "stl_slicer.h" +#include +#include +#include +#include +#include +#include +#include + +#pragma pack(push, 1) +struct StlHeader { + uint8_t header[80]; + uint32_t num_triangles; +}; + +struct StlTriangleRaw { + float normal[3]; + float v0[3]; + float v1[3]; + float v2[3]; + uint16_t attribute_byte_count; +}; +#pragma pack(pop) + +bool load_stl_binary(const std::string& filepath, StlMesh& mesh) { + std::ifstream file(filepath, std::ios::binary); + if (!file) { + std::cerr << "[Error] Cannot open STL file: " << filepath << std::endl; + return false; + } + + StlHeader header; + file.read(reinterpret_cast(&header), sizeof(StlHeader)); + if (!file) { + std::cerr << "[Error] Failed to read STL header" << std::endl; + return false; + } + + mesh.triangles.resize(header.num_triangles); + mesh.min_bound = { 1e30f, 1e30f, 1e30f }; + mesh.max_bound = { -1e30f, -1e30f, -1e30f }; + + std::vector raw_tris(header.num_triangles); + file.read(reinterpret_cast(raw_tris.data()), header.num_triangles * sizeof(StlTriangleRaw)); + + for (size_t i = 0; i < header.num_triangles; ++i) { + Triangle3D& tri = mesh.triangles[i]; + tri.normal = { raw_tris[i].normal[0], raw_tris[i].normal[1], raw_tris[i].normal[2] }; + tri.v0 = { raw_tris[i].v0[0], raw_tris[i].v0[1], raw_tris[i].v0[2] }; + tri.v1 = { raw_tris[i].v1[0], raw_tris[i].v1[1], raw_tris[i].v1[2] }; + tri.v2 = { raw_tris[i].v2[0], raw_tris[i].v2[1], raw_tris[i].v2[2] }; + + // Update bounds + const Vec3* verts[3] = { &tri.v0, &tri.v1, &tri.v2 }; + for (int v = 0; v < 3; ++v) { + mesh.min_bound.x = std::min(mesh.min_bound.x, verts[v]->x); + mesh.min_bound.y = std::min(mesh.min_bound.y, verts[v]->y); + mesh.min_bound.z = std::min(mesh.min_bound.z, verts[v]->z); + + mesh.max_bound.x = std::max(mesh.max_bound.x, verts[v]->x); + mesh.max_bound.y = std::max(mesh.max_bound.y, verts[v]->y); + mesh.max_bound.z = std::max(mesh.max_bound.z, verts[v]->z); + } + } + + return true; +} + +// Intersect a triangle with a horizontal plane Z = z0 +static bool intersect_triangle_plane(const Triangle3D& tri, float z0, LineSegment2D& out_seg) { + float z_min = std::min(std::min(tri.v0.z, tri.v1.z), tri.v2.z); + float z_max = std::max(std::max(tri.v0.z, tri.v1.z), tri.v2.z); + if (z0 < z_min || z0 > z_max || (z_min == z_max)) return false; + + const Vec3* v[3] = { &tri.v0, &tri.v1, &tri.v2 }; + GpuPoint2D pts[3]; + int count = 0; + + for (int i = 0; i < 3; ++i) { + int next = (i + 1) % 3; + float za = v[i]->z; + float zb = v[next]->z; + + if ((za <= z0 && zb >= z0) || (zb <= z0 && za >= z0)) { + if (std::abs(zb - za) > 1e-7f) { + float t = (z0 - za) / (zb - za); + pts[count].x = v[i]->x + t * (v[next]->x - v[i]->x); + pts[count].y = v[i]->y + t * (v[next]->y - v[i]->y); + count++; + if (count == 2) break; + } + } + } + + if (count == 2) { + float dx = pts[1].x - pts[0].x; + float dy = pts[1].y - pts[0].y; + if (dx * dx + dy * dy > 1e-8f) { + out_seg.p0 = pts[0]; + out_seg.p1 = pts[1]; + return true; + } + } + return false; +} + +// Quantized 2D point key for spatial hashing +struct PointKey { + int64_t x, y; + bool operator==(const PointKey& o) const { return x == o.x && y == o.y; } +}; + +struct PointKeyHash { + size_t operator()(const PointKey& k) const { + return (size_t)(k.x * 73856093 ^ k.y * 19349663); + } +}; + +static inline PointKey make_point_key(const GpuPoint2D& p, float scale = 1000.0f) { + return PointKey{ (int64_t)std::round(p.x * scale), (int64_t)std::round(p.y * scale) }; +} + +// Link unordered 2D line segments into closed polygon loops +static void link_segments_into_contours(const std::vector& segments, std::vector& out_contours) { + if (segments.empty()) return; + + std::unordered_multimap start_map; + std::vector used(segments.size(), false); + + for (size_t i = 0; i < segments.size(); ++i) { + start_map.insert({ make_point_key(segments[i].p0), i }); + } + + for (size_t i = 0; i < segments.size(); ++i) { + if (used[i]) continue; + + std::vector loop; + loop.push_back(segments[i].p0); + loop.push_back(segments[i].p1); + used[i] = true; + + GpuPoint2D current_end = segments[i].p1; + GpuPoint2D loop_start = segments[i].p0; + + while (true) { + PointKey key = make_point_key(current_end); + auto range = start_map.equal_range(key); + size_t next_idx = (size_t)-1; + + for (auto it = range.first; it != range.second; ++it) { + if (!used[it->second]) { + next_idx = it->second; + break; + } + } + + if (next_idx != (size_t)-1) { + used[next_idx] = true; + current_end = segments[next_idx].p1; + + // Check if loop has closed + PointKey end_key = make_point_key(current_end); + PointKey start_key = make_point_key(loop_start); + if (end_key == start_key || loop.size() > 50000) { + break; + } + loop.push_back(current_end); + } else { + break; + } + } + + if (loop.size() >= 3) { + GpuContour c; + c.num_points = loop.size(); + c.points = (GpuPoint2D*)malloc(c.num_points * sizeof(GpuPoint2D)); + for (size_t p = 0; p < c.num_points; ++p) { + c.points[p] = loop[p]; + } + c.is_hole = 0; // Can be determined by signed area if needed + out_contours.push_back(c); + } + } +} + +bool slice_stl_mesh( + const StlMesh& mesh, + float layer_height, + float first_layer_height, + std::vector& out_layers, + std::vector& z_levels, + double& slice_time_ms +) { + auto t0 = std::chrono::high_resolution_clock::now(); + + float z_min = mesh.min_bound.z; + float z_max = mesh.max_bound.z; + + float h_first = (first_layer_height > 0.0f) ? first_layer_height : layer_height; + z_levels.clear(); + for (float z = z_min + h_first; z <= z_max; z += layer_height) { + z_levels.push_back(z); + } + + size_t num_layers = z_levels.size(); + out_layers.resize(num_layers); + + // Parallel multi-threaded slicing across layers + unsigned int num_threads = std::max(1u, std::thread::hardware_concurrency()); + std::vector workers; + + auto worker_fn = [&](size_t start_idx, size_t end_idx) { + for (size_t l = start_idx; l < end_idx; ++l) { + float z0 = z_levels[l]; + std::vector segments; + segments.reserve(1000); + + for (const auto& tri : mesh.triangles) { + LineSegment2D seg; + if (intersect_triangle_plane(tri, z0, seg)) { + segments.push_back(seg); + } + } + + std::vector contours; + link_segments_into_contours(segments, contours); + + out_layers[l].z_height = z0; + out_layers[l].num_contours = contours.size(); + if (!contours.empty()) { + out_layers[l].contours = (GpuContour*)malloc(contours.size() * sizeof(GpuContour)); + for (size_t c = 0; c < contours.size(); ++c) { + out_layers[l].contours[c] = contours[c]; + } + } else { + out_layers[l].contours = nullptr; + } + } + }; + + size_t chunk = (num_layers + num_threads - 1) / num_threads; + for (unsigned int t = 0; t < num_threads; ++t) { + size_t start_idx = t * chunk; + size_t end_idx = std::min(start_idx + chunk, num_layers); + if (start_idx < end_idx) { + workers.emplace_back(worker_fn, start_idx, end_idx); + } + } + + for (auto& w : workers) { + w.join(); + } + + auto t1 = std::chrono::high_resolution_clock::now(); + slice_time_ms = std::chrono::duration(t1 - t0).count(); + + return true; +} + +void free_slice_layers(std::vector& layers) { + for (auto& layer : layers) { + if (layer.contours) { + for (size_t c = 0; c < layer.num_contours; ++c) { + if (layer.contours[c].points) { + free(layer.contours[c].points); + } + } + free(layer.contours); + layer.contours = nullptr; + } + } + layers.clear(); +} diff --git a/src/stl_slicer.h b/src/stl_slicer.h new file mode 100644 index 0000000..b210040 --- /dev/null +++ b/src/stl_slicer.h @@ -0,0 +1,44 @@ +#ifndef STL_SLICER_H +#define STL_SLICER_H + +#include "gpu_offset.h" +#include +#include + +struct Vec3 { + float x, y, z; +}; + +struct Triangle3D { + Vec3 normal; + Vec3 v0, v1, v2; +}; + +struct LineSegment2D { + GpuPoint2D p0; + GpuPoint2D p1; +}; + +struct StlMesh { + std::vector triangles; + Vec3 min_bound; + Vec3 max_bound; +}; + +// Load binary STL file into memory +bool load_stl_binary(const std::string& filepath, StlMesh& mesh); + +// Slices STL mesh into horizontal layers +bool slice_stl_mesh( + const StlMesh& mesh, + float layer_height, + float first_layer_height, + std::vector& out_layers, + std::vector& z_levels, + double& slice_time_ms +); + +// Free memory allocated for vector layers +void free_slice_layers(std::vector& layers); + +#endif // STL_SLICER_H diff --git a/src/vectorizer.h b/src/vectorizer.h new file mode 100644 index 0000000..b77a865 --- /dev/null +++ b/src/vectorizer.h @@ -0,0 +1,289 @@ +#ifndef VECTORIZER_H +#define VECTORIZER_H + +#include "gpu_offset.h" +#include +#include +#include +#include +#include + +inline int get_h_edge_id(int x, int y, int stride_w) { + return 2 * (y * stride_w + x) + 0; +} +inline int get_v_edge_id(int x, int y, int stride_w) { + return 2 * (y * stride_w + x) + 1; +} + +inline GpuPoint2D get_edge_point(int edge_id, int stride_w, float min_x, float min_y, float pitch_mm) { + int type = edge_id & 1; + int cell = edge_id >> 1; + int y = cell / stride_w; + int x = cell % stride_w; + + if (type == 0) { + return GpuPoint2D{ min_x + (x + 0.5f) * pitch_mm, min_y + (float)y * pitch_mm }; + } else { + return GpuPoint2D{ min_x + (float)x * pitch_mm, min_y + (y + 0.5f) * pitch_mm }; + } +} + +static inline void vectorize_mask_layer( + const uint8_t* mask, + int width, int height, + float min_x, float min_y, float pitch_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 uint8_t* row0 = mask + y * width; + const uint8_t* row1 = mask + (y + 1) * width; + + for (int x = 0; x < width - 1; ++x) { + int c0 = row0[x]; + int c1 = row0[x + 1]; + int c2 = row1[x + 1]; + int c3 = row1[x]; + + int case_idx = (c0) | (c1 << 1) | (c2 << 2) | (c3 << 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: add_directed_edge(e_left, e_top); add_directed_edge(e_bottom, e_right); 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: add_directed_edge(e_bottom, e_left); add_directed_edge(e_top, 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_edge_point(cur, stride_w, min_x, min_y, pitch_mm); + + 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-6f) { + 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) { + size_t n = loop_pts.size(); + float dx1 = loop_pts[0].x - loop_pts[n-1].x; + float dy1 = loop_pts[0].y - loop_pts[n-1].y; + float dx2 = loop_pts[1].x - loop_pts[0].x; + float dy2 = loop_pts[1].y - loop_pts[0].y; + if (std::abs(dx1 * dy2 - dy1 * dx2) < 1e-6f) { + loop_pts.erase(loop_pts.begin()); + } + + if (loop_pts.size() >= 3) { + GpuContour c; + c.num_points = loop_pts.size(); + c.points = (GpuPoint2D*)malloc(c.num_points * sizeof(GpuPoint2D)); + memcpy(c.points, loop_pts.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_bitpack_layer( + const unsigned int* bit_words, + int width, int height, int row_stride_words, + float min_x, float min_y, float pitch_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 unsigned int* row0 = bit_words + (size_t)y * row_stride_words; + const unsigned int* row1 = bit_words + (size_t)(y + 1) * row_stride_words; + + for (int w = 0; w < row_stride_words; ++w) { + unsigned int w0 = row0[w]; + unsigned int w1 = row1[w]; + + unsigned int w0_next = (w + 1 < row_stride_words) ? row0[w + 1] : 0; + unsigned int w1_next = (w + 1 < row_stride_words) ? row1[w + 1] : 0; + + if ((w0 == 0 && w1 == 0 && w0_next == 0 && w1_next == 0) || + (w0 == 0xFFFFFFFF && w1 == 0xFFFFFFFF && (w0_next & 1) && (w1_next & 1))) { + continue; + } + + int x_base = w * 32; + + for (int b = 0; b < 32; ++b) { + int x = x_base + b; + if (x >= width - 1) break; + + int c0 = (w0 >> b) & 1; + int c1 = (b < 31) ? ((w0 >> (b + 1)) & 1) : (w0_next & 1); + int c2 = (b < 31) ? ((w1 >> (b + 1)) & 1) : (w1_next & 1); + int c3 = (w1 >> b) & 1; + + int case_idx = c0 | (c1 << 1) | (c2 << 2) | (c3 << 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: add_directed_edge(e_left, e_top); add_directed_edge(e_bottom, e_right); 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: add_directed_edge(e_bottom, e_left); add_directed_edge(e_top, 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_edge_point(cur, stride_w, min_x, min_y, pitch_mm); + + 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-6f) { + 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) { + size_t n = loop_pts.size(); + float dx1 = loop_pts[0].x - loop_pts[n-1].x; + float dy1 = loop_pts[0].y - loop_pts[n-1].y; + float dx2 = loop_pts[1].x - loop_pts[0].x; + float dy2 = loop_pts[1].y - loop_pts[0].y; + if (std::abs(dx1 * dy2 - dy1 * dx2) < 1e-6f) { + loop_pts.erase(loop_pts.begin()); + } + + if (loop_pts.size() >= 3) { + GpuContour c; + c.num_points = loop_pts.size(); + c.points = (GpuPoint2D*)malloc(c.num_points * sizeof(GpuPoint2D)); + memcpy(c.points, loop_pts.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; + } +} + +#endif // VECTORIZER_H