diff --git a/CMakeLists.txt b/CMakeLists.txt index 1ab357c..72b4664 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,12 +7,13 @@ set(CMAKE_CUDA_STANDARD 14) set(CMAKE_CUDA_STANDARD_REQUIRED ON) # Cross-platform compiler settings +find_package(OpenMP REQUIRED) 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") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /O2 /W3 /MP /openmp") + set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler=/MD,/W3,/openmp") 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") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -pthread -fopenmp -fPIC") + set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler -fPIC,-Wall,-fopenmp") endif() # Set broad CUDA GPU architecture support (Turing, Ampere, Ada Lovelace, Hopper/Blackwell) @@ -83,7 +84,16 @@ execute_process( COMMAND "${Python3_EXECUTABLE}" -c "import pybind11; print(pybind11.get_include())" OUTPUT_VARIABLE _pybind11_include_dir OUTPUT_STRIP_TRAILING_WHITESPACE + ERROR_QUIET ) +if(NOT _pybind11_include_dir OR NOT EXISTS "${_pybind11_include_dir}") + execute_process( + COMMAND python3 -c "import pybind11; print(pybind11.get_include())" + OUTPUT_VARIABLE _pybind11_include_dir + OUTPUT_STRIP_TRAILING_WHITESPACE + ERROR_QUIET + ) +endif() execute_process( COMMAND "${Python3_EXECUTABLE}" -c "import sysconfig; print(sysconfig.get_config_var('EXT_SUFFIX') or '.so')" OUTPUT_VARIABLE _ext_suffix diff --git a/Makefile b/Makefile index b8b1238..b792438 100644 --- a/Makefile +++ b/Makefile @@ -2,8 +2,8 @@ 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 +CXXFLAGS = -O3 -std=c++17 -Wall -pthread -fopenmp -fPIC -Iinclude -Isrc +NVCCFLAGS = -O3 -std=c++14 -arch=$(ARCH) --use_fast_math -Xcompiler -fPIC,-Wall,-fopenmp -Iinclude -Isrc LIB_DIR = lib BIN_DIR = bin diff --git a/python/cudaslicer/__init__.py b/python/cudaslicer/__init__.py index 7e679df..2b86d1b 100644 --- a/python/cudaslicer/__init__.py +++ b/python/cudaslicer/__init__.py @@ -176,7 +176,8 @@ def offset_polygons( pad_mm: float = 1.5, return_infill_mask: bool = True, offset_distances: Optional[Sequence[float]] = None, - backend: Optional[str] = "auto" + backend: Optional[str] = "auto", + output_3d: bool = False ) -> Union[LayerResult, List[LayerResult]]: """ Offset 2D polygon loops across one or more layers (CUDA or CPU fallback). @@ -249,7 +250,8 @@ def offset_polygons( int(perimeters_count), float(layer_width), float(pitch_mm), float(pad_mm), bool(return_infill_mask), offsets_float, - backend_str + backend_str, + bool(output_3d) ) layer_results = [LayerResult(r) for r in raw_results] diff --git a/python/cudaslicer/bindings.cpp b/python/cudaslicer/bindings.cpp index e922957..12e17cf 100644 --- a/python/cudaslicer/bindings.cpp +++ b/python/cudaslicer/bindings.cpp @@ -107,7 +107,8 @@ py::list py_offset_polygon_layers( float pad_mm, bool export_infill_mask, const std::vector& offset_distances, - const std::string& backend + const std::string& backend, + bool output_3d = false ) { size_t num_layers = py_layers.size(); if (num_layers == 0) return py::list(); @@ -194,13 +195,28 @@ py::list py_offset_polygon_layers( 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; + if (output_3d && cnt.num_points > 0) { + py::array_t pts_arr({ (ssize_t)(cnt.num_points + 1), (ssize_t)3 }); + auto r = pts_arr.mutable_unchecked<2>(); + float z = results[l].z_height; + for (size_t p = 0; p < cnt.num_points; ++p) { + r(p, 0) = cnt.points[p].x; + r(p, 1) = cnt.points[p].y; + r(p, 2) = z; + } + r(cnt.num_points, 0) = cnt.points[0].x; + r(cnt.num_points, 1) = cnt.points[0].y; + r(cnt.num_points, 2) = z; + py_cnt_list.append(pts_arr); + } else { + 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_cnt_list.append(pts_arr); } py_perim_dict["contours"] = py_cnt_list; py_perimeters.append(py_perim_dict); @@ -211,13 +227,28 @@ py::list py_offset_polygon_layers( 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; + if (output_3d && cnt.num_points > 0) { + py::array_t pts_arr({ (ssize_t)(cnt.num_points + 1), (ssize_t)3 }); + auto r = pts_arr.mutable_unchecked<2>(); + float z = results[l].z_height; + for (size_t p = 0; p < cnt.num_points; ++p) { + r(p, 0) = cnt.points[p].x; + r(p, 1) = cnt.points[p].y; + r(p, 2) = z; + } + r(cnt.num_points, 0) = cnt.points[0].x; + r(cnt.num_points, 1) = cnt.points[0].y; + r(cnt.num_points, 2) = z; + py_infill_cnts.append(pts_arr); + } else { + 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_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; @@ -254,7 +285,8 @@ py::list py_slice_mesh_and_offset( const std::vector& offset_distances, const std::string& backend, float simplify_tolerance_mm, - bool subpixel_lerp + bool subpixel_lerp, + bool output_3d = false ) { py::buffer_info buf = triangles_arr.request(); size_t num_tris = 0; @@ -321,13 +353,28 @@ py::list py_slice_mesh_and_offset( 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; + if (output_3d && cnt.num_points > 0) { + py::array_t pts_arr({ (ssize_t)(cnt.num_points + 1), (ssize_t)3 }); + auto r = pts_arr.mutable_unchecked<2>(); + float z = results[l].z_height; + for (size_t p = 0; p < cnt.num_points; ++p) { + r(p, 0) = cnt.points[p].x; + r(p, 1) = cnt.points[p].y; + r(p, 2) = z; + } + r(cnt.num_points, 0) = cnt.points[0].x; + r(cnt.num_points, 1) = cnt.points[0].y; + r(cnt.num_points, 2) = z; + py_cnt_list.append(pts_arr); + } else { + 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_cnt_list.append(pts_arr); } py_perim_dict["contours"] = py_cnt_list; py_perimeters.append(py_perim_dict); @@ -338,13 +385,28 @@ py::list py_slice_mesh_and_offset( 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; + if (output_3d && cnt.num_points > 0) { + py::array_t pts_arr({ (ssize_t)(cnt.num_points + 1), (ssize_t)3 }); + auto r = pts_arr.mutable_unchecked<2>(); + float z = results[l].z_height; + for (size_t p = 0; p < cnt.num_points; ++p) { + r(p, 0) = cnt.points[p].x; + r(p, 1) = cnt.points[p].y; + r(p, 2) = z; + } + r(cnt.num_points, 0) = cnt.points[0].x; + r(cnt.num_points, 1) = cnt.points[0].y; + r(cnt.num_points, 2) = z; + py_infill_cnts.append(pts_arr); + } else { + 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_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; @@ -515,7 +577,8 @@ public: bool export_infill_mask = true, float simplify_tolerance_mm = 0.02f, bool subpixel_lerp = true, - const std::vector& offset_distances = std::vector() + const std::vector& offset_distances = std::vector(), + bool output_3d = false ) const { GpuOffsetConfig config = {0}; config.perimeters_count = perimeters_count; @@ -555,13 +618,28 @@ public: 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; + if (output_3d && cnt.num_points > 0) { + py::array_t pts_arr({ (ssize_t)(cnt.num_points + 1), (ssize_t)3 }); + auto r = pts_arr.mutable_unchecked<2>(); + float z = results[l].z_height; + for (size_t p = 0; p < cnt.num_points; ++p) { + r(p, 0) = cnt.points[p].x; + r(p, 1) = cnt.points[p].y; + r(p, 2) = z; + } + r(cnt.num_points, 0) = cnt.points[0].x; + r(cnt.num_points, 1) = cnt.points[0].y; + r(cnt.num_points, 2) = z; + py_cnt_list.append(pts_arr); + } else { + 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_cnt_list.append(pts_arr); } py_perim_dict["contours"] = py_cnt_list; py_perimeters.append(py_perim_dict); @@ -572,13 +650,28 @@ public: 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; + if (output_3d && cnt.num_points > 0) { + py::array_t pts_arr({ (ssize_t)(cnt.num_points + 1), (ssize_t)3 }); + auto r = pts_arr.mutable_unchecked<2>(); + float z = results[l].z_height; + for (size_t p = 0; p < cnt.num_points; ++p) { + r(p, 0) = cnt.points[p].x; + r(p, 1) = cnt.points[p].y; + r(p, 2) = z; + } + r(cnt.num_points, 0) = cnt.points[0].x; + r(cnt.num_points, 1) = cnt.points[0].y; + r(cnt.num_points, 2) = z; + py_infill_cnts.append(pts_arr); + } else { + 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_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; @@ -732,6 +825,7 @@ PYBIND11_MODULE(_cudaslicer_backend, m) { py::arg("export_infill_mask") = true, py::arg("offset_distances") = std::vector(), py::arg("backend") = "auto", + py::arg("output_3d") = false, "Compute polygon perimeter offsets and infill mask."); m.def("slice_mesh_and_offset", &py_slice_mesh_and_offset, @@ -748,6 +842,7 @@ PYBIND11_MODULE(_cudaslicer_backend, m) { py::arg("backend") = "auto", py::arg("simplify_tolerance_mm") = 0.0f, py::arg("subpixel_lerp") = true, + py::arg("output_3d") = false, "Mesh slicing, perimeter offsetting, and infill mask generation."); // GPU SDF Stack Factory Functions @@ -793,6 +888,7 @@ PYBIND11_MODULE(_cudaslicer_backend, m) { py::arg("simplify_tolerance_mm") = 0.02f, py::arg("subpixel_lerp") = true, py::arg("offset_distances") = std::vector(), + py::arg("output_3d") = false, "Vectorize SDF stack directly with sub-pixel LERP and RDP polygon simplification.") .def_property_readonly("__cuda_array_interface__", &PyGpuSdfStack::cuda_array_interface, "Zero-copy CUDA array interface for PyTorch / CuPy integration."); diff --git a/src/cpu_offset.cpp b/src/cpu_offset.cpp index 2274db3..ff6fd82 100644 --- a/src/cpu_offset.cpp +++ b/src/cpu_offset.cpp @@ -641,31 +641,82 @@ int cpu_offset_slice_mesh_and_offset( } } - // Compute exact continuous signed distance field from analytical layer segments + // Fast 2D Tile Spatial Binning for analytical segment distance on CPU + const int TILE_SZ = 16; + int num_tiles_x = (cols + TILE_SZ - 1) / TILE_SZ; + int num_tiles_y = (rows + TILE_SZ - 1) / TILE_SZ; + std::vector> tile_segs(num_tiles_x * num_tiles_y); + + const float SEARCH_R = 4.0f; + float inv_pitch = 1.0f / pitch; int num_segs = (int)layer_segs.size(); - for (int y = 0; y < rows; ++y) { - float py = bounds_min_y + (float)y * pitch; - for (int x = 0; x < cols; ++x) { - float px = bounds_min_x + (float)x * pitch; - float min_d2 = 1e20f; - for (int s = 0; s < num_segs; ++s) { - float sx0 = layer_segs[s].x0; - float sy0 = layer_segs[s].y0; - float sx1 = layer_segs[s].x1; - float sy1 = layer_segs[s].y1; - float dx = sx1 - sx0; - float dy = sy1 - sy0; - float l2 = dx * dx + dy * dy; - float t = (l2 > 1e-12f) ? ((px - sx0) * dx + (py - sy0) * dy) / l2 : 0.0f; - t = std::max(0.0f, std::min(1.0f, t)); - float qx = sx0 + t * dx; - float qy = sy0 + t * dy; - float d2 = (px - qx) * (px - qx) + (py - qy) * (py - qy); - if (d2 < min_d2) min_d2 = d2; + + for (int s = 0; s < num_segs; ++s) { + float sx0 = layer_segs[s].x0; + float sy0 = layer_segs[s].y0; + float sx1 = layer_segs[s].x1; + float sy1 = layer_segs[s].y1; + + float min_x = std::min(sx0, sx1) - SEARCH_R; + float max_x = std::max(sx0, sx1) + SEARCH_R; + float min_y = std::min(sy0, sy1) - SEARCH_R; + float max_y = std::max(sy0, sy1) + SEARCH_R; + + int tx0 = std::max(0, (int)((min_x - bounds_min_x) * inv_pitch) / TILE_SZ); + int tx1 = std::min(num_tiles_x - 1, (int)((max_x - bounds_min_x) * inv_pitch) / TILE_SZ); + int ty0 = std::max(0, (int)((min_y - bounds_min_y) * inv_pitch) / TILE_SZ); + int ty1 = std::min(num_tiles_y - 1, (int)((max_y - bounds_min_y) * inv_pitch) / TILE_SZ); + + for (int ty = ty0; ty <= ty1; ++ty) { + for (int tx = tx0; tx <= tx1; ++tx) { + tile_segs[ty * num_tiles_x + tx].push_back(s); + } + } + } + + for (int ty = 0; ty < num_tiles_y; ++ty) { + int y_start = ty * TILE_SZ; + int y_end = std::min(rows, y_start + TILE_SZ); + for (int tx = 0; tx < num_tiles_x; ++tx) { + int x_start = tx * TILE_SZ; + int x_end = std::min(cols, x_start + TILE_SZ); + const auto& candidates = tile_segs[ty * num_tiles_x + tx]; + + if (candidates.empty()) { + for (int y = y_start; y < y_end; ++y) { + for (int x = x_start; x < x_end; ++x) { + uint8_t m = layer_mask[(size_t)y * cols + x]; + layer_sdf[(size_t)y * cols + x] = (m > 0) ? SEARCH_R : -SEARCH_R; + } + } + continue; + } + + for (int y = y_start; y < y_end; ++y) { + float py = bounds_min_y + (float)y * pitch; + for (int x = x_start; x < x_end; ++x) { + float px = bounds_min_x + (float)x * pitch; + float min_d2 = 1e20f; + for (int s_idx : candidates) { + float sx0 = layer_segs[s_idx].x0; + float sy0 = layer_segs[s_idx].y0; + float sx1 = layer_segs[s_idx].x1; + float sy1 = layer_segs[s_idx].y1; + float dx = sx1 - sx0; + float dy = sy1 - sy0; + float l2 = dx * dx + dy * dy; + float t = (l2 > 1e-12f) ? ((px - sx0) * dx + (py - sy0) * dy) / l2 : 0.0f; + t = std::max(0.0f, std::min(1.0f, t)); + float qx = sx0 + t * dx; + float qy = sy0 + t * dy; + float d2 = (px - qx) * (px - qx) + (py - qy) * (py - qy); + if (d2 < min_d2) min_d2 = d2; + } + float d = (min_d2 < 1e19f) ? sqrtf(min_d2) : SEARCH_R; + uint8_t m = layer_mask[(size_t)y * cols + x]; + layer_sdf[(size_t)y * cols + x] = (m > 0) ? d : -d; + } } - float d = sqrtf(min_d2); - uint8_t m = layer_mask[(size_t)y * cols + x]; - layer_sdf[(size_t)y * cols + x] = (m > 0) ? d : -d; } } diff --git a/src/gpu_sdf.cu b/src/gpu_sdf.cu index b507893..0036492 100644 --- a/src/gpu_sdf.cu +++ b/src/gpu_sdf.cu @@ -12,6 +12,29 @@ #include #include #include +#include + +static float* g_pinned_sdf_buf = nullptr; +static size_t g_pinned_sdf_capacity = 0; +static std::mutex g_pinned_sdf_mutex; + +static float* get_pinned_sdf_buffer(size_t required_bytes) { + std::lock_guard lock(g_pinned_sdf_mutex); + if (g_pinned_sdf_capacity < required_bytes) { + if (g_pinned_sdf_buf) { + cudaFreeHost(g_pinned_sdf_buf); + g_pinned_sdf_buf = nullptr; + g_pinned_sdf_capacity = 0; + } + size_t alloc_bytes = required_bytes + 16 * 1024 * 1024; + if (cudaMallocHost(&g_pinned_sdf_buf, alloc_bytes) == cudaSuccess) { + g_pinned_sdf_capacity = alloc_bytes; + } else { + return nullptr; + } + } + return g_pinned_sdf_buf; +} #define CUDA_CHECK(call) do { \ cudaError_t err = call; \ @@ -702,11 +725,13 @@ __global__ void gpu_compute_exact_mesh_sdf_kernel( float bounds_min_y, float pitch_mm ) { - int x = blockIdx.x * blockDim.x + threadIdx.x; - int y = blockIdx.y * blockDim.y + threadIdx.y; - int l = blockIdx.z; + int tx = threadIdx.x; + int ty = threadIdx.y; + int tid = ty * blockDim.x + tx; - if (x >= cols || y >= rows || l >= num_layers) return; + int x = blockIdx.x * blockDim.x + tx; + int y = blockIdx.y * blockDim.y + ty; + int l = blockIdx.z; size_t layer_pixel_offset = (size_t)l * rows * cols; size_t idx = layer_pixel_offset + (size_t)y * cols + x; @@ -715,37 +740,89 @@ __global__ void gpu_compute_exact_mesh_sdf_kernel( if (seg_count > max_segs_per_layer) seg_count = max_segs_per_layer; if (seg_count == 0) { - out_sdf[idx] = -1e5f; + if (x < cols && y < rows && l < num_layers) { + out_sdf[idx] = -1e5f; + } return; } + // Tile bounding box in world mm + float tile_x0 = bounds_min_x + (float)(blockIdx.x * blockDim.x) * pitch_mm; + float tile_y0 = bounds_min_y + (float)(blockIdx.y * blockDim.y) * pitch_mm; + float tile_x1 = bounds_min_x + (float)min((blockIdx.x + 1) * (int)blockDim.x - 1, cols - 1) * pitch_mm; + float tile_y1 = bounds_min_y + (float)min((blockIdx.y + 1) * (int)blockDim.y - 1, rows - 1) * pitch_mm; + + const float SEARCH_R = 2.0f; // Max perimeter search range in mm + float bx0 = tile_x0 - SEARCH_R; + float by0 = tile_y0 - SEARCH_R; + float bx1 = tile_x1 + SEARCH_R; + float by1 = tile_y1 + SEARCH_R; + float px = bounds_min_x + (float)x * pitch_mm; float py = bounds_min_y + (float)y * pitch_mm; - float min_d2 = 1e10f; + const GpuSegment2D* segs = layer_segments + (size_t)l * max_segs_per_layer; - for (int s = 0; s < seg_count; ++s) { - GpuSegment2D seg = segs[s]; - float vx = seg.x1 - seg.x0; - float vy = seg.y1 - seg.y0; - float wx = px - seg.x0; - float wy = py - seg.y0; - float len_sq = vx * vx + vy * vy; - float t = (len_sq > 1e-10f) ? (wx * vx + wy * vy) / len_sq : 0.0f; - if (t < 0.0f) t = 0.0f; - else if (t > 1.0f) t = 1.0f; - float qx = seg.x0 + t * vx; - float qy = seg.y0 + t * vy; - float dx = px - qx; - float dy = py - qy; - float d2 = dx * dx + dy * dy; - if (d2 < min_d2) min_d2 = d2; + __shared__ GpuSegment2D s_active_segs[256]; + __shared__ int s_active_count; + + for (int chunk_start = 0; chunk_start < seg_count; chunk_start += 256) { + if (tid == 0) { + s_active_count = 0; + } + __syncthreads(); + + int s_idx = chunk_start + tid; + if (s_idx < seg_count) { + GpuSegment2D seg = segs[s_idx]; + float sx_min = fminf(seg.x0, seg.x1); + float sx_max = fmaxf(seg.x0, seg.x1); + float sy_min = fminf(seg.y0, seg.y1); + float sy_max = fmaxf(seg.y0, seg.y1); + + if (sx_max >= bx0 && sx_min <= bx1 && sy_max >= by0 && sy_min <= by1) { + int pos = atomicAdd(&s_active_count, 1); + if (pos < 256) { + s_active_segs[pos] = seg; + } + } + } + __syncthreads(); + + int num_active = min(s_active_count, 256); + if (num_active > 0 && x < cols && y < rows) { + for (int k = 0; k < num_active; ++k) { + GpuSegment2D seg = s_active_segs[k]; + float vx = seg.x1 - seg.x0; + float vy = seg.y1 - seg.y0; + float wx = px - seg.x0; + float wy = py - seg.y0; + float len_sq = vx * vx + vy * vy; + float t = (len_sq > 1e-10f) ? (wx * vx + wy * vy) / len_sq : 0.0f; + if (t < 0.0f) t = 0.0f; + else if (t > 1.0f) t = 1.0f; + float qx = seg.x0 + t * vx; + float qy = seg.y0 + t * vy; + float dx = px - qx; + float dy = py - qy; + float d2 = dx * dx + dy * dy; + if (d2 < min_d2) min_d2 = d2; + } + } + __syncthreads(); } - float d = sqrtf(min_d2); - uint8_t m = masks[idx]; - out_sdf[idx] = (m > 0) ? d : -d; + if (x < cols && y < rows && l < num_layers) { + uint8_t m = masks[idx]; + float d; + if (min_d2 < 1e9f) { + d = sqrtf(min_d2); + } else { + d = SEARCH_R; + } + out_sdf[idx] = (m > 0) ? d : -d; + } } GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_mesh( @@ -805,26 +882,26 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_mesh( if (max_segs_per_layer < 256) max_segs_per_layer = 256; size_t total_segs = (size_t)num_layers * max_segs_per_layer; - GpuTriangle3D* d_triangles = nullptr; - unsigned int* d_parity = nullptr; - uint8_t* d_masks = nullptr; - GpuSegment2D* d_layer_segments = nullptr; - int* d_layer_seg_counts = nullptr; + size_t tri_bytes_aligned = (tri_bytes + 255) & ~255ULL; + size_t parity_bytes_aligned = (total_words * sizeof(unsigned int) + 255) & ~255ULL; + size_t masks_bytes_aligned = (total_volume * sizeof(uint8_t) + 255) & ~255ULL; + size_t segs_bytes_aligned = (total_segs * sizeof(GpuSegment2D) + 255) & ~255ULL; + size_t counts_bytes_aligned = (num_layers * sizeof(int) + 255) & ~255ULL; - if (cudaMalloc(&d_triangles, tri_bytes) != cudaSuccess || - cudaMalloc(&d_parity, total_words * sizeof(unsigned int)) != cudaSuccess || - cudaMalloc(&d_masks, total_volume * sizeof(uint8_t)) != cudaSuccess || - cudaMalloc(&d_layer_segments, total_segs * sizeof(GpuSegment2D)) != cudaSuccess || - cudaMalloc(&d_layer_seg_counts, num_layers * sizeof(int)) != cudaSuccess) { - if (d_triangles) cudaFree(d_triangles); - if (d_parity) cudaFree(d_parity); - if (d_masks) cudaFree(d_masks); - if (d_layer_segments) cudaFree(d_layer_segments); - if (d_layer_seg_counts) cudaFree(d_layer_seg_counts); + size_t workspace_total = tri_bytes_aligned + parity_bytes_aligned + masks_bytes_aligned + segs_bytes_aligned + counts_bytes_aligned; + + uint8_t* d_workspace = nullptr; + if (cudaMalloc(&d_workspace, workspace_total) != cudaSuccess) { gpu_sdf_destroy(stack); return nullptr; } + GpuTriangle3D* d_triangles = (GpuTriangle3D*)(d_workspace); + unsigned int* d_parity = (unsigned int*)(d_workspace + tri_bytes_aligned); + uint8_t* d_masks = (uint8_t*)(d_workspace + tri_bytes_aligned + parity_bytes_aligned); + GpuSegment2D* d_layer_segments = (GpuSegment2D*)(d_workspace + tri_bytes_aligned + parity_bytes_aligned + masks_bytes_aligned); + int* d_layer_seg_counts = (int*)(d_workspace + tri_bytes_aligned + parity_bytes_aligned + masks_bytes_aligned + segs_bytes_aligned); + cudaMemsetAsync(d_parity, 0, total_words * sizeof(unsigned int), stack->stream); cudaMemsetAsync(d_layer_seg_counts, 0, num_layers * sizeof(int), stack->stream); cudaMemcpyAsync(d_triangles, triangles, tri_bytes, cudaMemcpyHostToDevice, stack->stream); @@ -851,12 +928,7 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_mesh( ); cudaStreamSynchronize(stack->stream); - - cudaFree(d_triangles); - cudaFree(d_parity); - cudaFree(d_masks); - cudaFree(d_layer_segments); - cudaFree(d_layer_seg_counts); + cudaFree(d_workspace); return stack; } @@ -1102,6 +1174,284 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_shell( return out_stack; } +struct GpuDirectedSegment { + int src_edge; + int dst_edge; + GpuPoint2D pt; +}; + +__device__ inline GpuPoint2D interpolate_gpu_edge( + int edge_idx, int x, int y, + float c0, float c1, float c2, float c3, + float iso_val, float min_x, float min_y, float pitch_mm, int subpixel_lerp +) { + if (!subpixel_lerp) { + if (edge_idx == 0) return GpuPoint2D{ min_x + (x + 0.5f) * pitch_mm, min_y + (float)y * pitch_mm }; + if (edge_idx == 1) return GpuPoint2D{ min_x + (float)(x + 1) * pitch_mm, min_y + (y + 0.5f) * pitch_mm }; + if (edge_idx == 2) return GpuPoint2D{ min_x + (x + 0.5f) * pitch_mm, min_y + (float)(y + 1) * pitch_mm }; + /* edge_idx == 3 */return GpuPoint2D{ min_x + (float)x * pitch_mm, min_y + (y + 0.5f) * pitch_mm }; + } + + float t = 0.5f; + if (edge_idx == 0) { // Bottom (c0 -> c1) + float denom = c1 - c0; + if (fabsf(denom) > 1e-7f) { + t = (iso_val - c0) / denom; + t = fminf(fmaxf(t, 0.0f), 1.0f); + } + return GpuPoint2D{ min_x + ((float)x + t) * pitch_mm, min_y + (float)y * pitch_mm }; + } else if (edge_idx == 1) { // Right (c1 -> c2) + float denom = c2 - c1; + if (fabsf(denom) > 1e-7f) { + t = (iso_val - c1) / denom; + t = fminf(fmaxf(t, 0.0f), 1.0f); + } + return GpuPoint2D{ min_x + (float)(x + 1) * pitch_mm, min_y + ((float)y + t) * pitch_mm }; + } else if (edge_idx == 2) { // Top (c3 -> c2) + float denom = c2 - c3; + if (fabsf(denom) > 1e-7f) { + t = (iso_val - c3) / denom; + t = fminf(fmaxf(t, 0.0f), 1.0f); + } + return GpuPoint2D{ min_x + ((float)x + t) * pitch_mm, min_y + (float)(y + 1) * pitch_mm }; + } else { // Left (c0 -> c3) + float denom = c3 - c0; + if (fabsf(denom) > 1e-7f) { + t = (iso_val - c0) / denom; + t = fminf(fmaxf(t, 0.0f), 1.0f); + } + return GpuPoint2D{ min_x + (float)x * pitch_mm, min_y + ((float)y + t) * pitch_mm }; + } +} + +__global__ void gpu_ms_extract_segments_kernel( + const float* __restrict__ d_sdf, + int cols, int rows, int num_layers, + float bounds_min_x, float bounds_min_y, float pitch_mm, + const float* __restrict__ d_offsets, + int num_offsets, + int subpixel_lerp, + int max_segs_per_pass, + GpuDirectedSegment* __restrict__ d_out_segments, + int* __restrict__ d_out_segment_counts +) { + int x = blockIdx.x * blockDim.x + threadIdx.x; + int y = blockIdx.y * blockDim.y + threadIdx.y; + int pass_id = blockIdx.z; + + int total_passes = num_layers * num_offsets; + if (pass_id >= total_passes) return; + + int l = pass_id / num_offsets; + int o = pass_id % num_offsets; + float iso_val = d_offsets[o]; + + int stride_w = cols + 1; + size_t layer_offset = (size_t)l * rows * cols; + const float* layer_sdf = d_sdf + layer_offset; + + if (x >= cols - 1 || y >= rows - 1) return; + + const float* row0 = layer_sdf + (size_t)y * cols; + const float* row1 = layer_sdf + (size_t)(y + 1) * cols; + + float c0 = row0[x]; + float c1 = row0[x + 1]; + float c2 = row1[x + 1]; + float c3 = row1[x]; + + int b0 = (c0 >= iso_val) ? 1 : 0; + int b1 = (c1 >= iso_val) ? 1 : 0; + int b2 = (c2 >= iso_val) ? 1 : 0; + int b3 = (c3 >= iso_val) ? 1 : 0; + + int case_idx = b0 | (b1 << 1) | (b2 << 2) | (b3 << 3); + if (case_idx == 0 || case_idx == 15) return; + + int e_bottom = 2 * (y * stride_w + x) + 0; + int e_top = 2 * ((y + 1) * stride_w + x) + 0; + int e_left = 2 * (y * stride_w + x) + 1; + int e_right = 2 * (y * stride_w + (x + 1)) + 1; + + auto emit_segment = [&](int local_src, int global_src, int global_dst) { + GpuPoint2D pt = interpolate_gpu_edge(local_src, x, y, c0, c1, c2, c3, iso_val, bounds_min_x, bounds_min_y, pitch_mm, subpixel_lerp); + int pos = atomicAdd(&d_out_segment_counts[pass_id], 1); + if (pos < max_segs_per_pass) { + d_out_segments[(size_t)pass_id * max_segs_per_pass + pos] = GpuDirectedSegment{ global_src, global_dst, pt }; + } + }; + + switch (case_idx) { + case 1: emit_segment(3, e_left, e_bottom); break; + case 2: emit_segment(0, e_bottom, e_right); break; + case 3: emit_segment(3, e_left, e_right); break; + case 4: emit_segment(1, e_right, e_top); break; + case 5: { + float center = (c0 + c1 + c2 + c3) * 0.25f; + if (center >= iso_val) { + emit_segment(3, e_left, e_top); + emit_segment(0, e_bottom, e_right); + } else { + emit_segment(3, e_left, e_bottom); + emit_segment(1, e_right, e_top); + } + break; + } + case 6: emit_segment(0, e_bottom, e_top); break; + case 7: emit_segment(3, e_left, e_top); break; + case 8: emit_segment(2, e_top, e_left); break; + case 9: emit_segment(2, e_top, e_bottom); break; + case 10: { + float center = (c0 + c1 + c2 + c3) * 0.25f; + if (center >= iso_val) { + emit_segment(0, e_bottom, e_left); + emit_segment(2, e_top, e_right); + } else { + emit_segment(2, e_top, e_left); + emit_segment(0, e_bottom, e_right); + } + break; + } + case 11: emit_segment(2, e_top, e_right); break; + case 12: emit_segment(1, e_right, e_left); break; + case 13: emit_segment(1, e_right, e_bottom); break; + case 14: emit_segment(0, e_bottom, e_left); break; + } +} + +__global__ void gpu_prefix_scan_counts_kernel( + const int* __restrict__ d_counts, + int* __restrict__ d_offsets, + int total_passes +) { + if (threadIdx.x == 0 && blockIdx.x == 0) { + int sum = 0; + for (int i = 0; i < total_passes; ++i) { + d_offsets[i] = sum; + sum += d_counts[i]; + } + d_offsets[total_passes] = sum; + } +} + +__global__ void gpu_compact_segments_kernel( + const GpuDirectedSegment* __restrict__ d_src_segments, + const int* __restrict__ d_counts, + const int* __restrict__ d_pass_offsets, + int max_segs_per_pass, + int total_passes, + GpuDirectedSegment* __restrict__ d_dst_segments +) { + int p = blockIdx.x; + int tid = threadIdx.x; + if (p >= total_passes) return; + + int count = d_counts[p]; + int dst_offset = d_pass_offsets[p]; + const GpuDirectedSegment* src = d_src_segments + (size_t)p * max_segs_per_pass; + GpuDirectedSegment* dst = d_dst_segments + dst_offset; + + for (int i = tid; i < count; i += blockDim.x) { + dst[i] = src[i]; + } +} + +struct CompactEdgeMap { + int cap; + int mask; + int* keys; + int* dsts; + GpuPoint2D* pts; + uint8_t* visited; + + CompactEdgeMap() : cap(0), mask(0), keys(nullptr), dsts(nullptr), pts(nullptr), visited(nullptr) {} + + void ensure_capacity(int required) { + int target = 128; + while (target < required * 2) target <<= 1; + if (target <= cap) return; + + if (keys) { free(keys); free(dsts); free(pts); free(visited); } + cap = target; + mask = cap - 1; + keys = (int*)malloc(cap * sizeof(int)); + dsts = (int*)malloc(cap * sizeof(int)); + pts = (GpuPoint2D*)malloc(cap * sizeof(GpuPoint2D)); + visited = (uint8_t*)malloc(cap * sizeof(uint8_t)); + memset(keys, -1, cap * sizeof(int)); + memset(visited, 0, cap * sizeof(uint8_t)); + } + void free_mem() { + if (keys) { free(keys); free(dsts); free(pts); free(visited); } + keys = nullptr; cap = 0; + } + inline void insert(int src, int dst, GpuPoint2D pt) { + int idx = ((unsigned int)src * 2654435761u) & mask; + while (keys[idx] != -1 && keys[idx] != src) { + idx = (idx + 1) & mask; + } + keys[idx] = src; + dsts[idx] = dst; + pts[idx] = pt; + visited[idx] = 0; + } + inline int find_slot(int src) const { + int idx = ((unsigned int)src * 2654435761u) & mask; + while (keys[idx] != -1) { + if (keys[idx] == src) return idx; + idx = (idx + 1) & mask; + } + return -1; + } + inline void clear_slot(int idx) { + keys[idx] = -1; + visited[idx] = 0; + } +}; + +struct GpuVecWorkspace { + GpuDirectedSegment* d_segments = nullptr; + GpuDirectedSegment* d_compact_segments = nullptr; + int* d_counts = nullptr; + int* d_offsets = nullptr; + float* d_iso_offsets = nullptr; + int* h_offsets = nullptr; + GpuDirectedSegment* h_compact_segments = nullptr; + size_t max_total_segs = 0; + int total_passes_cap = 0; + int num_offsets_cap = 0; + std::mutex mtx; + + void ensure(int total_passes, int max_segs_per_pass, int num_offsets) { + std::lock_guard lock(mtx); + size_t needed_segs = (size_t)total_passes * max_segs_per_pass; + if (needed_segs > max_total_segs) { + if (d_segments) cudaFree(d_segments); + if (d_compact_segments) cudaFree(d_compact_segments); + if (h_compact_segments) cudaFreeHost(h_compact_segments); + cudaMalloc(&d_segments, needed_segs * sizeof(GpuDirectedSegment)); + cudaMalloc(&d_compact_segments, needed_segs * sizeof(GpuDirectedSegment)); + cudaMallocHost(&h_compact_segments, needed_segs * sizeof(GpuDirectedSegment)); + max_total_segs = needed_segs; + } + if (total_passes > total_passes_cap) { + if (d_counts) cudaFree(d_counts); + if (d_offsets) cudaFree(d_offsets); + if (h_offsets) cudaFreeHost(h_offsets); + cudaMalloc(&d_counts, total_passes * sizeof(int)); + cudaMalloc(&d_offsets, (total_passes + 1) * sizeof(int)); + cudaMallocHost(&h_offsets, (total_passes + 1) * sizeof(int)); + total_passes_cap = total_passes; + } + if (num_offsets > num_offsets_cap) { + if (d_iso_offsets) cudaFree(d_iso_offsets); + cudaMalloc(&d_iso_offsets, num_offsets * sizeof(float)); + num_offsets_cap = num_offsets; + } + } +}; +static GpuVecWorkspace g_vec_workspace; + GPU_SDF_API int gpu_sdf_vectorize( GpuSdfStackHandle stack, const GpuOffsetConfig* config, @@ -1137,117 +1487,188 @@ GPU_SDF_API int gpu_sdf_vectorize( resolved_offsets[P] = (float)P * W; } size_t num_offsets = resolved_offsets.size(); + int total_passes = (int)(num_layers * num_offsets); float simplify_tol = config->simplify_tolerance_mm; int subpixel_lerp = config->subpixel_lerp; + int max_segs_per_pass = 8192; - size_t total_voxels = (size_t)num_layers * rows * cols; - size_t sdf_bytes = total_voxels * sizeof(float); + g_vec_workspace.ensure(total_passes, max_segs_per_pass, (int)num_offsets); - float* h_pinned_sdf = nullptr; - CUDA_CHECK_INT(cudaMallocHost(&h_pinned_sdf, sdf_bytes)); - CUDA_CHECK_INT(cudaMemcpyAsync(h_pinned_sdf, stack->d_sdf, sdf_bytes, - cudaMemcpyDeviceToHost, stack->stream)); - CUDA_CHECK_INT(cudaStreamSynchronize(stack->stream)); + cudaMemcpyAsync(g_vec_workspace.d_iso_offsets, resolved_offsets.data(), num_offsets * sizeof(float), cudaMemcpyHostToDevice, stack->stream); + cudaMemsetAsync(g_vec_workspace.d_counts, 0, total_passes * sizeof(int), stack->stream); + + auto t_k_0 = std::chrono::high_resolution_clock::now(); + + dim3 block(16, 16); + dim3 grid((cols + 15) / 16, (rows + 15) / 16, total_passes); + gpu_ms_extract_segments_kernel<<stream>>>( + stack->d_sdf, cols, rows, (int)num_layers, bounds_min_x, bounds_min_y, pitch_mm, + g_vec_workspace.d_iso_offsets, (int)num_offsets, subpixel_lerp, max_segs_per_pass, + g_vec_workspace.d_segments, g_vec_workspace.d_counts + ); + gpu_prefix_scan_counts_kernel<<<1, 1, 0, stack->stream>>>( + g_vec_workspace.d_counts, g_vec_workspace.d_offsets, total_passes + ); + gpu_compact_segments_kernel<<stream>>>( + g_vec_workspace.d_segments, g_vec_workspace.d_counts, g_vec_workspace.d_offsets, + max_segs_per_pass, total_passes, g_vec_workspace.d_compact_segments + ); + cudaStreamSynchronize(stack->stream); + auto t_k_1 = std::chrono::high_resolution_clock::now(); + + // Fast PCIe download of active segments + auto t_d2h_0 = std::chrono::high_resolution_clock::now(); + cudaMemcpyAsync(g_vec_workspace.h_offsets, g_vec_workspace.d_offsets, (total_passes + 1) * sizeof(int), cudaMemcpyDeviceToHost, stack->stream); + cudaStreamSynchronize(stack->stream); + + size_t total_active_segments = (size_t)g_vec_workspace.h_offsets[total_passes]; + cudaMemcpyAsync(g_vec_workspace.h_compact_segments, g_vec_workspace.d_compact_segments, + total_active_segments * sizeof(GpuDirectedSegment), cudaMemcpyDeviceToHost, stack->stream); + cudaStreamSynchronize(stack->stream); + auto t_d2h_1 = std::chrono::high_resolution_clock::now(); std::vector z_heights_um(num_layers); gpu_sdf_get_layer_z_um(stack, z_heights_um.data()); - // Multi-threaded C++ Sub-Pixel LERP Vectorization + // Multi-threaded L1-Resident Loop Chaining & Simplification 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; + for (size_t l = 0; l < num_layers; ++l) { + res_array[l].z_height = (float)z_heights_um[l] * 0.001f; + res_array[l].num_perimeters = P; + res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter)); + memset(res_array[l].perimeters, 0, P * sizeof(GpuOffsetPerimeter)); + for (int o = 0; o < P; ++o) { + res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o]; + } + res_array[l].infill.offset_distance_mm = resolved_offsets[P]; + res_array[l].infill.width = cols; + res_array[l].infill.height = rows; + } - 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); + #pragma omp parallel + { + CompactEdgeMap map; + std::vector loop_pts; + loop_pts.reserve(2000); + std::vector active_slots; + active_slots.reserve(4000); + std::vector pass_contours; + pass_contours.reserve(32); - for (size_t l = start_layer; l < end_layer; ++l) { - res_array[l].z_height = (float)z_heights_um[l] * 0.001f; - res_array[l].num_perimeters = P; - res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter)); + #pragma omp for schedule(dynamic, 8) + for (int p = 0; p < total_passes; ++p) { + size_t l = (size_t)(p / (int)num_offsets); + int o = p % (int)num_offsets; - const float* layer_sdf_ptr = h_pinned_sdf + l * ((size_t)rows * cols); + int start_idx = g_vec_workspace.h_offsets[p]; + int end_idx = g_vec_workspace.h_offsets[p + 1]; + int n_segs = end_idx - start_idx; - for (int o = 0; o < P; ++o) { - float iso_thresh = resolved_offsets[o]; + if (n_segs < 3) { + continue; + } - std::vector contours; - vectorize_sdf_layer(layer_sdf_ptr, cols, rows, bounds_min_x, bounds_min_y, pitch_mm, - iso_thresh, simplify_tol, subpixel_lerp, - next_edge_buf, active_edges_buf, visited_buf, contours); + map.ensure_capacity(n_segs); + const GpuDirectedSegment* segs = g_vec_workspace.h_compact_segments + start_idx; + active_slots.clear(); + pass_contours.clear(); - res_array[l].perimeters[o].offset_distance_mm = iso_thresh; - res_array[l].perimeters[o].num_contours = contours.size(); - if (!contours.empty()) { - res_array[l].perimeters[o].contours = (GpuContour*)malloc(contours.size() * sizeof(GpuContour)); - memcpy(res_array[l].perimeters[o].contours, contours.data(), contours.size() * sizeof(GpuContour)); + for (int i = 0; i < n_segs; ++i) { + map.insert(segs[i].src_edge, segs[i].dst_edge, segs[i].pt); + } + + for (int i = 0; i < n_segs; ++i) { + int slot = map.find_slot(segs[i].src_edge); + if (slot == -1 || map.visited[slot]) continue; + + loop_pts.clear(); + int cur_slot = slot; + + while (cur_slot != -1 && !map.visited[cur_slot]) { + map.visited[cur_slot] = 1; + active_slots.push_back(cur_slot); + GpuPoint2D pt = map.pts[cur_slot]; + + size_t n = loop_pts.size(); + if (n >= 2) { + float dx1 = loop_pts[n-1].x - loop_pts[n-2].x; + float dy1 = loop_pts[n-1].y - loop_pts[n-2].y; + float dx2 = pt.x - loop_pts[n-1].x; + float dy2 = pt.y - loop_pts[n-1].y; + + float cross = dx1 * dy2 - dy1 * dx2; + if (std::abs(cross) < 1e-7f) { + loop_pts[n-1] = pt; // simplify collinear vertex + } else { + loop_pts.push_back(pt); + } + } else { + loop_pts.push_back(pt); + } + + int next_src = map.dsts[cur_slot]; + cur_slot = map.find_slot(next_src); + if (cur_slot == slot) break; + if (loop_pts.size() > 100000) break; + } + + if (loop_pts.size() >= 3) { + std::vector simplified = (simplify_tol > 1e-6f) + ? simplify_contour_rdp(loop_pts, simplify_tol) + : loop_pts; + if (simplified.size() >= 3) { + GpuContour c; + c.num_points = simplified.size(); + c.points = (GpuPoint2D*)malloc(c.num_points * sizeof(GpuPoint2D)); + memcpy(c.points, simplified.data(), c.num_points * sizeof(GpuPoint2D)); + c.is_hole = 0; + pass_contours.push_back(c); + } + } + } + + for (int s : active_slots) { + map.clear_slot(s); + } + + if (!pass_contours.empty()) { + if (o < P) { + res_array[l].perimeters[o].num_contours = pass_contours.size(); + res_array[l].perimeters[o].contours = (GpuContour*)malloc(pass_contours.size() * sizeof(GpuContour)); + memcpy(res_array[l].perimeters[o].contours, pass_contours.data(), pass_contours.size() * sizeof(GpuContour)); } else { - res_array[l].perimeters[o].contours = nullptr; + res_array[l].infill.num_contours = pass_contours.size(); + res_array[l].infill.contours = (GpuContour*)malloc(pass_contours.size() * sizeof(GpuContour)); + memcpy(res_array[l].infill.contours, pass_contours.data(), pass_contours.size() * sizeof(GpuContour)); } } - - // Infill Contours (Innermost offset boundary) - float infill_iso_thresh = resolved_offsets[P]; - std::vector infill_contours; - vectorize_sdf_layer(layer_sdf_ptr, cols, rows, bounds_min_x, bounds_min_y, pitch_mm, - infill_iso_thresh, simplify_tol, subpixel_lerp, - next_edge_buf, active_edges_buf, visited_buf, infill_contours); - - res_array[l].infill.offset_distance_mm = infill_iso_thresh; - res_array[l].infill.width = cols; - res_array[l].infill.height = rows; - res_array[l].infill.num_contours = infill_contours.size(); - if (!infill_contours.empty()) { - res_array[l].infill.contours = (GpuContour*)malloc(infill_contours.size() * sizeof(GpuContour)); - memcpy(res_array[l].infill.contours, infill_contours.data(), infill_contours.size() * sizeof(GpuContour)); - } else { - res_array[l].infill.contours = nullptr; - } - - if (config->export_infill_mask) { - uint8_t* infill_mask = (uint8_t*)malloc((size_t)cols * rows); - for (int i = 0; i < cols * rows; ++i) { - infill_mask[i] = (layer_sdf_ptr[i] >= infill_iso_thresh) ? 1 : 0; - } - res_array[l].infill.mask = infill_mask; - } else { - res_array[l].infill.mask = nullptr; - } } - }; - - size_t chunk = (num_layers + num_threads - 1) / num_threads; - for (unsigned int t = 0; t < num_threads; ++t) { - size_t s = t * chunk; - size_t e = std::min(s + chunk, num_layers); - if (s < e) workers.emplace_back(vec_worker, s, e); + map.free_mem(); } - 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(); + double t_d2h_ms = std::chrono::duration(t_d2h_1 - t_d2h_0).count(); + double t_kernel_ms = std::chrono::duration(t_k_1 - t_k_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_gpu_thresh_ms = t_kernel_ms; + stats->t_d2h_ms = t_d2h_ms; stats->t_vectorize_ms = t_vec_ms; stats->t_total_pipeline_ms = t_total_ms; } - cudaFreeHost(h_pinned_sdf); - *out_results = res_array; return GPU_OFFSET_SUCCESS; } } // extern "C" + diff --git a/src/vectorizer.h b/src/vectorizer.h index d579af4..312851d 100644 --- a/src/vectorizer.h +++ b/src/vectorizer.h @@ -40,15 +40,16 @@ inline GpuPoint2D get_sdf_edge_point( int cell = edge_id >> 1; int stride_w = width + 1; int y = cell / stride_w; - int x = cell % stride_w; + int x = cell - y * stride_w; if (type == 0) { // Horizontal edge: from (x, y) to (x + 1, y) float d0 = -1e5f; float d1 = -1e5f; - if (y >= 0 && y < height) { - if (x >= 0 && x < width) d0 = sdf[y * width + x]; - if (x + 1 >= 0 && x + 1 < width) d1 = sdf[y * width + (x + 1)]; + if ((unsigned int)y < (unsigned int)height) { + size_t row_idx = (size_t)y * width; + if ((unsigned int)x < (unsigned int)width) d0 = sdf[row_idx + x]; + if ((unsigned int)(x + 1) < (unsigned int)width) d1 = sdf[row_idx + (x + 1)]; } float t = 0.5f; float denom = d1 - d0; @@ -62,9 +63,9 @@ inline GpuPoint2D get_sdf_edge_point( // Vertical edge: from (x, y) to (x, y + 1) float d0 = -1e5f; float d3 = -1e5f; - if (x >= 0 && x < width) { - if (y >= 0 && y < height) d0 = sdf[y * width + x]; - if (y + 1 >= 0 && y + 1 < height) d3 = sdf[(y + 1) * width + x]; + if ((unsigned int)x < (unsigned int)width) { + if ((unsigned int)y < (unsigned int)height) d0 = sdf[(size_t)y * width + x]; + if ((unsigned int)(y + 1) < (unsigned int)height) d3 = sdf[(size_t)(y + 1) * width + x]; } float t = 0.5f; float denom = d3 - d0; @@ -162,6 +163,70 @@ static inline std::vector simplify_contour_rdp( return simplified; } +// C++ Seam Alignment: positions starting vertex at sharpest concave corner (or sharpest angle) +static inline std::vector align_contour_seam_sharpest_cpp(const std::vector& loop) { + size_t M = loop.size(); + if (M < 3) return loop; + + // 1. Calculate signed area + double signed_area = 0.0; + for (size_t i = 0; i < M; ++i) { + size_t next_i = (i + 1 == M) ? 0 : i + 1; + signed_area += (double)loop[i].x * (double)loop[next_i].y - (double)loop[next_i].x * (double)loop[i].y; + } + signed_area *= 0.5; + if (std::abs(signed_area) < 1e-9) return loop; + + // 2. Find sharpest concave angle or fallback to sharpest angle + size_t best_concave_idx = 0; + float min_concave_cos = 2.0f; + bool found_concave = false; + + size_t best_any_idx = 0; + float min_any_cos = 2.0f; + + for (size_t i = 0; i < M; ++i) { + size_t prev_i = (i == 0) ? M - 1 : i - 1; + size_t next_i = (i + 1 == M) ? 0 : i + 1; + + float v_in_x = loop[i].x - loop[prev_i].x; + float v_in_y = loop[i].y - loop[prev_i].y; + float v_out_x = loop[next_i].x - loop[i].x; + float v_out_y = loop[next_i].y - loop[i].y; + + float l_in = std::sqrt(v_in_x * v_in_x + v_in_y * v_in_y); + float l_out = std::sqrt(v_out_x * v_out_x + v_out_y * v_out_y); + + if (l_in < 1e-6f || l_out < 1e-6f) continue; + + float cross = v_in_x * v_out_y - v_in_y * v_out_x; + bool is_concave = (signed_area > 0.0) ? (cross < -1e-6f) : (cross > 1e-6f); + + float dot = v_in_x * v_out_x + v_in_y * v_out_y; + float cos_angle = dot / (l_in * l_out); + + if (cos_angle < min_any_cos) { + min_any_cos = cos_angle; + best_any_idx = i; + } + + if (is_concave && cos_angle < min_concave_cos) { + min_concave_cos = cos_angle; + best_concave_idx = i; + found_concave = true; + } + } + + size_t best_idx = found_concave ? best_concave_idx : best_any_idx; + if (best_idx == 0) return loop; + + std::vector aligned; + aligned.reserve(M); + for (size_t i = best_idx; i < M; ++i) aligned.push_back(loop[i]); + for (size_t i = 0; i < best_idx; ++i) aligned.push_back(loop[i]); + return aligned; +} + static inline void vectorize_sdf_layer( const float* sdf, int width, int height, @@ -183,64 +248,93 @@ static inline void vectorize_sdf_layer( active_edges_buf.push_back(e0); }; - for (int y = 0; y < height - 1; ++y) { - const float* row0 = sdf + y * width; - const float* row1 = sdf + (y + 1) * width; + const int BLOCK = 8; + float max_diag_mm = 1.4142f * (float)BLOCK * pitch_mm; - for (int x = 0; x < width - 1; ++x) { - float c0 = row0[x]; - float c1 = row0[x + 1]; - float c2 = row1[x + 1]; - float c3 = row1[x]; + for (int by = 0; by < height - 1; by += BLOCK) { + int ey = std::min(by + BLOCK, height - 1); + for (int bx = 0; bx < width - 1; bx += BLOCK) { + int ex = std::min(bx + BLOCK, width - 1); - int b0 = (c0 >= iso_val) ? 1 : 0; - int b1 = (c1 >= iso_val) ? 1 : 0; - int b2 = (c2 >= iso_val) ? 1 : 0; - int b3 = (c3 >= iso_val) ? 1 : 0; + // Fast Lipschitz block culling: test 4 corners + float c_tl = sdf[by * width + bx]; + float c_tr = sdf[by * width + ex]; + float c_bl = sdf[ey * width + bx]; + float c_br = sdf[ey * width + ex]; - int case_idx = b0 | (b1 << 1) | (b2 << 2) | (b3 << 3); - if (case_idx == 0 || case_idx == 15) continue; + if (c_tl > iso_val + max_diag_mm && + c_tr > iso_val + max_diag_mm && + c_bl > iso_val + max_diag_mm && + c_br > iso_val + max_diag_mm) { + continue; // 100% Guaranteed no isocontour in this block + } + if (c_tl < iso_val - max_diag_mm && + c_tr < iso_val - max_diag_mm && + c_bl < iso_val - max_diag_mm && + c_br < iso_val - max_diag_mm) { + continue; // 100% Guaranteed no isocontour in this block + } - 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); + for (int y = by; y < ey; ++y) { + const float* row0 = sdf + y * width; + const float* row1 = sdf + (y + 1) * width; - switch (case_idx) { - case 1: add_directed_edge(e_left, e_bottom); break; - case 2: add_directed_edge(e_bottom, e_right); break; - case 3: add_directed_edge(e_left, e_right); break; - case 4: add_directed_edge(e_right, e_top); break; - case 5: { - float center = (c0 + c1 + c2 + c3) * 0.25f; - if (center >= iso_val) { - add_directed_edge(e_left, e_top); - add_directed_edge(e_bottom, e_right); - } else { - add_directed_edge(e_left, e_bottom); - add_directed_edge(e_right, e_top); + for (int x = bx; x < ex; ++x) { + float c0 = row0[x]; + float c1 = row0[x + 1]; + float c2 = row1[x + 1]; + float c3 = row1[x]; + + int b0 = (c0 >= iso_val) ? 1 : 0; + int b1 = (c1 >= iso_val) ? 1 : 0; + int b2 = (c2 >= iso_val) ? 1 : 0; + int b3 = (c3 >= iso_val) ? 1 : 0; + + int case_idx = b0 | (b1 << 1) | (b2 << 2) | (b3 << 3); + if (case_idx == 0 || case_idx == 15) continue; + + int e_bottom = get_h_edge_id(x, y, stride_w); + int e_top = get_h_edge_id(x, y + 1, stride_w); + int e_left = get_v_edge_id(x, y, stride_w); + int e_right = get_v_edge_id(x + 1, y, stride_w); + + switch (case_idx) { + case 1: add_directed_edge(e_left, e_bottom); break; + case 2: add_directed_edge(e_bottom, e_right); break; + case 3: add_directed_edge(e_left, e_right); break; + case 4: add_directed_edge(e_right, e_top); break; + case 5: { + float center = (c0 + c1 + c2 + c3) * 0.25f; + if (center >= iso_val) { + add_directed_edge(e_left, e_top); + add_directed_edge(e_bottom, e_right); + } else { + add_directed_edge(e_left, e_bottom); + add_directed_edge(e_right, e_top); + } + break; + } + case 6: add_directed_edge(e_bottom, e_top); break; + case 7: add_directed_edge(e_left, e_top); break; + case 8: add_directed_edge(e_top, e_left); break; + case 9: add_directed_edge(e_top, e_bottom); break; + case 10: { + float center = (c0 + c1 + c2 + c3) * 0.25f; + if (center >= iso_val) { + add_directed_edge(e_bottom, e_left); + add_directed_edge(e_top, e_right); + } else { + add_directed_edge(e_top, e_left); + add_directed_edge(e_bottom, e_right); + } + break; + } + case 11: add_directed_edge(e_top, e_right); break; + case 12: add_directed_edge(e_right, e_left); break; + case 13: add_directed_edge(e_right, e_bottom); break; + case 14: add_directed_edge(e_bottom, e_left); break; } - break; } - case 6: add_directed_edge(e_bottom, e_top); break; - case 7: add_directed_edge(e_left, e_top); break; - case 8: add_directed_edge(e_top, e_left); break; - case 9: add_directed_edge(e_top, e_bottom); break; - case 10: { - float center = (c0 + c1 + c2 + c3) * 0.25f; - if (center >= iso_val) { - add_directed_edge(e_bottom, e_left); - add_directed_edge(e_top, e_right); - } else { - add_directed_edge(e_top, e_left); - add_directed_edge(e_bottom, e_right); - } - break; - } - case 11: add_directed_edge(e_top, e_right); break; - case 12: add_directed_edge(e_right, e_left); break; - case 13: add_directed_edge(e_right, e_bottom); break; - case 14: add_directed_edge(e_bottom, e_left); break; } } } @@ -292,10 +386,11 @@ static inline void vectorize_sdf_layer( ? simplify_contour_rdp(loop_pts, simplify_tolerance_mm) : loop_pts; if (simplified.size() >= 3) { + std::vector aligned = align_contour_seam_sharpest_cpp(simplified); GpuContour c; - c.num_points = simplified.size(); + c.num_points = aligned.size(); c.points = (GpuPoint2D*)malloc(c.num_points * sizeof(GpuPoint2D)); - memcpy(c.points, simplified.data(), c.num_points * sizeof(GpuPoint2D)); + memcpy(c.points, aligned.data(), c.num_points * sizeof(GpuPoint2D)); c.is_hole = 0; out_contours.push_back(c); } @@ -412,10 +507,11 @@ static inline void vectorize_mask_layer( } if (loop_pts.size() >= 3) { + std::vector aligned = align_contour_seam_sharpest_cpp(loop_pts); GpuContour c; - c.num_points = loop_pts.size(); + c.num_points = aligned.size(); c.points = (GpuPoint2D*)malloc(c.num_points * sizeof(GpuPoint2D)); - memcpy(c.points, loop_pts.data(), c.num_points * sizeof(GpuPoint2D)); + memcpy(c.points, aligned.data(), c.num_points * sizeof(GpuPoint2D)); c.is_hole = 0; out_contours.push_back(c); } @@ -550,10 +646,11 @@ static inline void vectorize_bitpack_layer( } if (loop_pts.size() >= 3) { + std::vector aligned = align_contour_seam_sharpest_cpp(loop_pts); GpuContour c; - c.num_points = loop_pts.size(); + c.num_points = aligned.size(); c.points = (GpuPoint2D*)malloc(c.num_points * sizeof(GpuPoint2D)); - memcpy(c.points, loop_pts.data(), c.num_points * sizeof(GpuPoint2D)); + memcpy(c.points, aligned.data(), c.num_points * sizeof(GpuPoint2D)); c.is_hole = 0; out_contours.push_back(c); }