Merge branch 'main' of ssh://git.c2d.lv:2222/nicolas/cudaslicer

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_CXX_FLAGS "${CMAKE_CXX_FLAGS} /O2 /W3 /MP /openmp")
-    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler=/MD,/W3")
+    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_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -pthread -fopenmp -fPIC")
-    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler -fPIC,-Wall")
+    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

Makefile

@@ -2,8 +2,8 @@ NVCC ?= nvcc
 CXX ?= g++
 ARCH ?= sm_86
 
-CXXFLAGS = -O3 -std=c++17 -Wall -pthread -fPIC -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 -Iinclude -Isrc
+NVCCFLAGS = -O3 -std=c++14 -arch=$(ARCH) --use_fast_math -Xcompiler -fPIC,-Wall,-fopenmp -Iinclude -Isrc
 
 LIB_DIR = lib
 BIN_DIR = bin

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]

python/cudaslicer/bindings.cpp

@@ -107,7 +107,8 @@ py::list py_offset_polygon_layers(
     float pad_mm,
     bool export_infill_mask,
     const std::vector<float>& 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,6 +195,20 @@ 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];
+                if (output_3d && cnt.num_points > 0) {
+                    py::array_t<float> 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<float> 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) {
@@ -202,6 +217,7 @@ py::list py_offset_polygon_layers(
                     }
                     py_cnt_list.append(pts_arr);
                 }
+            }
             py_perim_dict["contours"] = py_cnt_list;
             py_perimeters.append(py_perim_dict);
         }
@@ -211,6 +227,20 @@ 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];
+            if (output_3d && cnt.num_points > 0) {
+                py::array_t<float> 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<float> 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) {
@@ -219,6 +249,7 @@ py::list py_offset_polygon_layers(
                 }
                 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<float>& 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,6 +353,20 @@ 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];
+                if (output_3d && cnt.num_points > 0) {
+                    py::array_t<float> 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<float> 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) {
@@ -329,6 +375,7 @@ py::list py_slice_mesh_and_offset(
                     }
                     py_cnt_list.append(pts_arr);
                 }
+            }
             py_perim_dict["contours"] = py_cnt_list;
             py_perimeters.append(py_perim_dict);
         }
@@ -338,6 +385,20 @@ 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];
+            if (output_3d && cnt.num_points > 0) {
+                py::array_t<float> 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<float> 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) {
@@ -346,6 +407,7 @@ py::list py_slice_mesh_and_offset(
                 }
                 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<float>& offset_distances = std::vector<float>()
+        const std::vector<float>& offset_distances = std::vector<float>(),
+        bool output_3d = false
     ) const {
         GpuOffsetConfig config = {0};
         config.perimeters_count = perimeters_count;
@@ -555,6 +618,20 @@ 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];
+                    if (output_3d && cnt.num_points > 0) {
+                        py::array_t<float> 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<float> 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) {
@@ -563,6 +640,7 @@ public:
                         }
                         py_cnt_list.append(pts_arr);
                     }
+                }
                 py_perim_dict["contours"] = py_cnt_list;
                 py_perimeters.append(py_perim_dict);
             }
@@ -572,6 +650,20 @@ 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];
+                if (output_3d && cnt.num_points > 0) {
+                    py::array_t<float> 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<float> 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) {
@@ -580,6 +672,7 @@ public:
                     }
                     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<float>(),
           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<float>(),
+             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.");

src/cpu_offset.cpp

@@ -641,18 +641,67 @@ 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<std::vector<int>> 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 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;
@@ -663,11 +712,13 @@ int cpu_offset_slice_mesh_and_offset(
                                 float d2 = (px - qx) * (px - qx) + (py - qy) * (py - qy);
                                 if (d2 < min_d2) min_d2 = d2;
                             }
-                    float d = sqrtf(min_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;
                         }
                     }
+                }
+            }
 
             res_array[l].z_height = z0;
             res_array[l].num_perimeters = P;

src/gpu_sdf.cu

@@ -12,6 +12,29 @@
 #include <thread>
 #include <map>
 #include <set>
+#include <mutex>
+
+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<std::mutex> 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 tx = threadIdx.x;
-    int y = blockIdx.y * blockDim.y + threadIdx.y;
+    int ty = threadIdx.y;
-    int l = blockIdx.z;
+    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,18 +740,60 @@ __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) {
+        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) {
+    __shared__ GpuSegment2D s_active_segs[256];
-        GpuSegment2D seg = segs[s];
+    __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;
@@ -742,11 +809,21 @@ __global__ void gpu_compute_exact_mesh_sdf_kernel(
                 float d2 = dx * dx + dy * dy;
                 if (d2 < min_d2) min_d2 = d2;
             }
+        }
+        __syncthreads();
+    }
 
-    float d = sqrtf(min_d2);
+    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(
     const GpuTriangle3D* triangles,
@@ -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;
+    size_t tri_bytes_aligned = (tri_bytes + 255) & ~255ULL;
-    unsigned int* d_parity = nullptr;
+    size_t parity_bytes_aligned = (total_words * sizeof(unsigned int) + 255) & ~255ULL;
-    uint8_t* d_masks = nullptr;
+    size_t masks_bytes_aligned = (total_volume * sizeof(uint8_t) + 255) & ~255ULL;
-    GpuSegment2D* d_layer_segments = nullptr;
+    size_t segs_bytes_aligned = (total_segs * sizeof(GpuSegment2D) + 255) & ~255ULL;
-    int* d_layer_seg_counts = nullptr;
+    size_t counts_bytes_aligned = (num_layers * sizeof(int) + 255) & ~255ULL;
 
-    if (cudaMalloc(&d_triangles, tri_bytes) != cudaSuccess ||
+    size_t workspace_total = tri_bytes_aligned + parity_bytes_aligned + masks_bytes_aligned + segs_bytes_aligned + counts_bytes_aligned;
-        cudaMalloc(&d_parity, total_words * sizeof(unsigned int)) != cudaSuccess ||
+
-        cudaMalloc(&d_masks, total_volume * sizeof(uint8_t)) != cudaSuccess ||
+    uint8_t* d_workspace = nullptr;
-        cudaMalloc(&d_layer_segments, total_segs * sizeof(GpuSegment2D)) != cudaSuccess ||
+    if (cudaMalloc(&d_workspace, workspace_total) != 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);
         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_workspace);
-    cudaFree(d_triangles);
-    cudaFree(d_parity);
-    cudaFree(d_masks);
-    cudaFree(d_layer_segments);
-    cudaFree(d_layer_seg_counts);
 
     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<std::mutex> 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;
+    g_vec_workspace.ensure(total_passes, max_segs_per_pass, (int)num_offsets);
-    size_t sdf_bytes = total_voxels * sizeof(float);
 
-    float* h_pinned_sdf = nullptr;
+    cudaMemcpyAsync(g_vec_workspace.d_iso_offsets, resolved_offsets.data(), num_offsets * sizeof(float), cudaMemcpyHostToDevice, stack->stream);
-    CUDA_CHECK_INT(cudaMallocHost(&h_pinned_sdf, sdf_bytes));
+    cudaMemsetAsync(g_vec_workspace.d_counts, 0, total_passes * sizeof(int), stack->stream);
-    CUDA_CHECK_INT(cudaMemcpyAsync(h_pinned_sdf, stack->d_sdf, sdf_bytes,
+
-                                   cudaMemcpyDeviceToHost, stack->stream));
+    auto t_k_0 = std::chrono::high_resolution_clock::now();
-    CUDA_CHECK_INT(cudaStreamSynchronize(stack->stream));
+
+    dim3 block(16, 16);
+    dim3 grid((cols + 15) / 16, (rows + 15) / 16, total_passes);
+    gpu_ms_extract_segments_kernel<<<grid, block, 0, stack->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<<<total_passes, 128, 0, stack->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<int32_t> 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());
+    for (size_t l = 0; l < num_layers; ++l) {
-    std::vector<std::thread> 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<int> next_edge_buf(max_edges, -1);
-        std::vector<int> active_edges_buf;
-        active_edges_buf.reserve(cols * rows / 16);
-        std::vector<uint8_t> visited_buf(max_edges, 0);
-
-        for (size_t l = start_layer; l < end_layer; ++l) {
         res_array[l].z_height = (float)z_heights_um[l] * 0.001f;
         res_array[l].num_perimeters = P;
         res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter));
-
+        memset(res_array[l].perimeters, 0, P * sizeof(GpuOffsetPerimeter));
-            const float* layer_sdf_ptr = h_pinned_sdf + l * ((size_t)rows * cols);
-
         for (int o = 0; o < P; ++o) {
-                float iso_thresh = resolved_offsets[o];
+            res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o];
-
-                std::vector<GpuContour> 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);
-
-                res_array[l].perimeters[o].offset_distance_mm = iso_thresh;
-                res_array[l].perimeters[o].num_contours = contours.size();
-                if (!contours.empty()) {
-                    res_array[l].perimeters[o].contours = (GpuContour*)malloc(contours.size() * sizeof(GpuContour));
-                    memcpy(res_array[l].perimeters[o].contours, contours.data(), contours.size() * sizeof(GpuContour));
-                } else {
-                    res_array[l].perimeters[o].contours = nullptr;
         }
-            }
+        res_array[l].infill.offset_distance_mm = resolved_offsets[P];
-
-            // Infill Contours (Innermost offset boundary)
-            float infill_iso_thresh = resolved_offsets[P];
-            std::vector<GpuContour> 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) {
+    #pragma omp parallel
-                uint8_t* infill_mask = (uint8_t*)malloc((size_t)cols * rows);
+    {
-                for (int i = 0; i < cols * rows; ++i) {
+        CompactEdgeMap map;
-                    infill_mask[i] = (layer_sdf_ptr[i] >= infill_iso_thresh) ? 1 : 0;
+        std::vector<GpuPoint2D> loop_pts;
-                }
+        loop_pts.reserve(2000);
-                res_array[l].infill.mask = infill_mask;
+        std::vector<int> active_slots;
-            } else {
+        active_slots.reserve(4000);
-                res_array[l].infill.mask = nullptr;
+        std::vector<GpuContour> pass_contours;
-            }
+        pass_contours.reserve(32);
-        }
-    };
 
-    size_t chunk = (num_layers + num_threads - 1) / num_threads;
+        #pragma omp for schedule(dynamic, 8)
-    for (unsigned int t = 0; t < num_threads; ++t) {
+        for (int p = 0; p < total_passes; ++p) {
-        size_t s = t * chunk;
+            size_t l = (size_t)(p / (int)num_offsets);
-        size_t e = std::min(s + chunk, num_layers);
+            int o = p % (int)num_offsets;
-        if (s < e) workers.emplace_back(vec_worker, s, e);
+
+            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;
+
+            if (n_segs < 3) {
+                continue;
+            }
+
+            map.ensure_capacity(n_segs);
+            const GpuDirectedSegment* segs = g_vec_workspace.h_compact_segments + start_idx;
+            active_slots.clear();
+            pass_contours.clear();
+
+            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<GpuPoint2D> 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].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));
+                }
+            }
+        }
+        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<double, std::milli>(t_vec_1 - t_vec_0).count();
+    double t_d2h_ms = std::chrono::duration<double, std::milli>(t_d2h_1 - t_d2h_0).count();
+    double t_kernel_ms = std::chrono::duration<double, std::milli>(t_k_1 - t_k_0).count();
 
     auto t_total_end = std::chrono::high_resolution_clock::now();
     double t_total_ms = std::chrono::duration<double, std::milli>(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"
+

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 ((unsigned int)y < (unsigned int)height) {
-            if (x >= 0 && x < width) d0 = sdf[y * width + x];
+            size_t row_idx = (size_t)y * width;
-            if (x + 1 >= 0 && x + 1 < width) d1 = sdf[y * width + (x + 1)];
+            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 ((unsigned int)x < (unsigned int)width) {
-            if (y >= 0 && y < height) d0 = sdf[y * width + x];
+            if ((unsigned int)y < (unsigned int)height) d0 = sdf[(size_t)y * width + x];
-            if (y + 1 >= 0 && y + 1 < height) d3 = sdf[(y + 1) * 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<GpuPoint2D> simplify_contour_rdp(
     return simplified;
 }
 
+// C++ Seam Alignment: positions starting vertex at sharpest concave corner (or sharpest angle)
+static inline std::vector<GpuPoint2D> align_contour_seam_sharpest_cpp(const std::vector<GpuPoint2D>& 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<GpuPoint2D> 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,11 +248,38 @@ static inline void vectorize_sdf_layer(
         active_edges_buf.push_back(e0);
     };
 
-    for (int y = 0; y < height - 1; ++y) {
+    const int BLOCK = 8;
+    float max_diag_mm = 1.4142f * (float)BLOCK * pitch_mm;
+
+    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);
+
+            // 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];
+
+            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
+            }
+
+            for (int y = by; y < ey; ++y) {
                 const float* row0 = sdf + y * width;
                 const float* row1 = sdf + (y + 1) * width;
 
-        for (int x = 0; x < width - 1; ++x) {
+                for (int x = bx; x < ex; ++x) {
                     float c0 = row0[x];
                     float c1 = row0[x + 1];
                     float c2 = row1[x + 1];
@@ -244,6 +336,8 @@ static inline void vectorize_sdf_layer(
                     }
                 }
             }
+        }
+    }
 
     if (active_edges_buf.empty()) return;
 
@@ -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<GpuPoint2D> 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<GpuPoint2D> 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<GpuPoint2D> 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);
             }