uo

CMakeLists.txt

@@ -40,7 +40,7 @@ set_target_properties(cudaslicer_offset PROPERTIES
     CUDA_SEPARABLE_COMPILATION OFF
 )
 target_compile_definitions(cudaslicer_offset PRIVATE GPU_OFFSET_EXPORTS)
-target_link_libraries(cudaslicer_offset PUBLIC cudart)
+target_link_libraries(cudaslicer_offset PUBLIC cudart OpenMP::OpenMP_CXX)
 
 # -----------------------------------------------------------------------------
 # 2. Standalone C/C++ Static Library
@@ -54,7 +54,7 @@ set_target_properties(cudaslicer_offset_static PROPERTIES
     OUTPUT_NAME "cudaslicer_offset"
     CUDA_SEPARABLE_COMPILATION OFF
 )
-target_link_libraries(cudaslicer_offset_static PUBLIC cudart)
+target_link_libraries(cudaslicer_offset_static PUBLIC cudart OpenMP::OpenMP_CXX)
 
 # -----------------------------------------------------------------------------
 # 3. C++ Benchmark Binary
@@ -108,7 +108,7 @@ set_target_properties(_cudaslicer_backend PROPERTIES
     SUFFIX "${_ext_suffix}"
 )
 target_include_directories(_cudaslicer_backend PRIVATE "${_python3_include_dir}" "${_pybind11_include_dir}")
-target_link_libraries(_cudaslicer_backend PRIVATE cudaslicer_offset_static cudart)
+target_link_libraries(_cudaslicer_backend PRIVATE cudaslicer_offset_static cudart OpenMP::OpenMP_CXX)
 
 # -----------------------------------------------------------------------------
 # Output Directories

python/cudaslicer/bindings.cpp

@@ -79,10 +79,14 @@ py::array_t<uint8_t> py_offset_masks(
     }
 
     GpuOffsetStats stats = {0};
-    int ret = gpu_offset_compute_masks(
+    int ret;
+    {
+        py::gil_scoped_release release;
+        ret = gpu_offset_compute_masks(
             ptr_in, width, height, num_layers, pitch_mm,
             resolved_offsets.data(), num_offsets, ptr_out, &stats
         );
+    }
 
     if (!backend.empty()) {
         gpu_offset_set_backend(prev_backend);
@@ -169,11 +173,15 @@ py::list py_offset_polygon_layers(
     GpuLayerResult* results = nullptr;
     GpuOffsetStats stats = {0};
 
-    int ret = gpu_offset_polygon_layers(
+    int ret;
+    {
+        py::gil_scoped_release release;
+        ret = gpu_offset_polygon_layers(
             c_layers.data(), num_layers,
             min_x, min_y, max_x, max_y,
             &config, &results, &stats
         );
+    }
 
     if (ret == GPU_OFFSET_ERR_NO_CUDA) {
         throw std::runtime_error("CUDA backend was requested, but no functional CUDA-capable device was found.");
@@ -326,12 +334,16 @@ py::list py_slice_mesh_and_offset(
     size_t out_num_layers = 0;
     GpuOffsetStats stats = {0};
 
-    int ret = gpu_offset_slice_mesh_and_offset(
+    int ret;
+    {
+        py::gil_scoped_release release;
+        ret = gpu_offset_slice_mesh_and_offset(
             tris_ptr, num_tris,
             min_x, min_y, min_z, max_x, max_y, max_z,
             layer_height, first_layer_height, &config,
             &results, &out_num_layers, &stats
         );
+    }
 
     if (ret == GPU_OFFSET_ERR_NO_CUDA) {
         throw std::runtime_error("CUDA backend was requested, but no functional CUDA-capable device was found.");

src/cpu_offset.cpp

@@ -6,7 +6,7 @@
 #include <chrono>
 #include <cstring>
 #include <algorithm>
-#include <thread>
+#include <omp.h>
 
 // -----------------------------------------------------------------------------
 // Fast Separable 2D Exact Euclidean Distance Transform on Binary Grid
@@ -21,8 +21,8 @@ static inline void edt_1d_binary_col(const uint8_t* col_mask, int height, int st
         if (col_mask[(size_t)y * stride] == 0) {
             last = y;
         }
-        float dy = (float)(y - last);
+        int dy = y - last;
-        col_g[y] = (last >= 0) ? (dy * dy) : INF;
+        col_g[y] = (last >= 0) ? (float)(dy * dy) : INF;
     }
     // Backward pass
     last = 1000000;
@@ -30,8 +30,34 @@ static inline void edt_1d_binary_col(const uint8_t* col_mask, int height, int st
         if (col_mask[(size_t)y * stride] == 0) {
             last = y;
         }
-        float dy = (float)(last - y);
+        int dy = last - y;
-        float d2 = (last < height) ? (dy * dy) : INF;
+        float d2 = (last < height) ? (float)(dy * dy) : INF;
+        if (d2 < col_g[y]) {
+            col_g[y] = d2;
+        }
+    }
+}
+
+// Pass 1 on inverted binary column (mask > 0 are seeds)
+static inline void edt_1d_binary_col_inverted(const uint8_t* col_mask, int height, int stride, float* col_g) {
+    const float INF = 1e12f;
+    // Forward pass
+    int last = -1000000;
+    for (int y = 0; y < height; ++y) {
+        if (col_mask[(size_t)y * stride] != 0) {
+            last = y;
+        }
+        int dy = y - last;
+        col_g[y] = (last >= 0) ? (float)(dy * dy) : INF;
+    }
+    // Backward pass
+    last = 1000000;
+    for (int y = height - 1; y >= 0; --y) {
+        if (col_mask[(size_t)y * stride] != 0) {
+            last = y;
+        }
+        int dy = last - y;
+        float d2 = (last < height) ? (float)(dy * dy) : INF;
         if (d2 < col_g[y]) {
             col_g[y] = d2;
         }
@@ -78,12 +104,19 @@ static inline void edt_2d_layer(
     std::vector<float>& row_f,
     std::vector<float>& row_d,
     std::vector<int>& v_buf,
-    std::vector<float>& z_buf
+    std::vector<float>& z_buf,
+    bool invert = false
 ) {
     // Pass 1: Columns
+    if (!invert) {
         for (int x = 0; x < width; ++x) {
             edt_1d_binary_col(mask + x, height, width, g_col_buf.data() + (size_t)x * height);
         }
+    } else {
+        for (int x = 0; x < width; ++x) {
+            edt_1d_binary_col_inverted(mask + x, height, width, g_col_buf.data() + (size_t)x * height);
+        }
+    }
 
     // Pass 2: Rows
     for (int y = 0; y < height; ++y) {
@@ -98,13 +131,23 @@ static inline void edt_2d_layer(
 }
 
 // -----------------------------------------------------------------------------
-// CPU Multi-Threaded Polygon Rasterizer
+// Structures for Fast Exact Analytical Segment Distance on CPU
 // -----------------------------------------------------------------------------
-struct CpuEdge2D {
+struct CpuSegment2D {
-    float x1, y1;
+    float x0, y0, x1, y1;
-    float x2, y2;
+    float vx, vy;
+    float inv_len_sq;
+    float min_x, max_x, min_y, max_y;
 };
 
+struct TileSegNode {
+    int seg_idx;
+    int next;
+};
+
+// -----------------------------------------------------------------------------
+// CPU Multi-Threaded Polygon Rasterizer (Scanline Parity)
+// -----------------------------------------------------------------------------
 void cpu_rasterize_polygon_layers_internal(
     const GpuSliceLayer* in_layers,
     size_t num_layers,
@@ -112,60 +155,80 @@ void cpu_rasterize_polygon_layers_internal(
     int width, int height,
     float min_x, float min_y, float pitch_mm
 ) {
-    unsigned int num_threads = std::max(1u, std::thread::hardware_concurrency());
+    int row_stride_words = (width + 31) / 32;
-    std::vector<std::thread> workers;
+    float inv_pitch = 1.0f / pitch_mm;
 
-    auto worker_fn = [&](size_t start_layer, size_t end_layer) {
+    #pragma omp parallel
-        for (size_t l = start_layer; l < end_layer; ++l) {
+    {
-            uint8_t* mask = out_masks + l * ((size_t)width * height);
+        std::vector<unsigned int> parity_words((size_t)height * row_stride_words);
+
+        #pragma omp for schedule(dynamic, 4)
+        for (int l = 0; l < (int)num_layers; ++l) {
+            uint8_t* mask = out_masks + (size_t)l * ((size_t)width * height);
             const GpuSliceLayer& layer = in_layers[l];
+            std::fill(parity_words.begin(), parity_words.end(), 0u);
 
-            std::vector<CpuEdge2D> edges;
             for (size_t c = 0; c < layer.num_contours; ++c) {
                 const GpuContour& cnt = layer.contours[c];
                 if (cnt.num_points < 3) continue;
+
                 for (size_t p = 0; p < cnt.num_points; ++p) {
-                    size_t p_next = (p + 1) % cnt.num_points;
+                    size_t p_next = (p + 1 == cnt.num_points) ? 0 : p + 1;
-                    edges.push_back({ cnt.points[p].x, cnt.points[p].y,
+                    float px0 = cnt.points[p].x;
-                                      cnt.points[p_next].x, cnt.points[p_next].y });
+                    float py0 = cnt.points[p].y;
+                    float px1 = cnt.points[p_next].x;
+                    float py1 = cnt.points[p_next].y;
+
+                    float gx0 = (px0 - min_x) * inv_pitch;
+                    float gy0 = (py0 - min_y) * inv_pitch;
+                    float gx1 = (px1 - min_x) * inv_pitch;
+                    float gy1 = (py1 - min_y) * inv_pitch;
+
+                    if (gy0 == gy1) continue;
+                    if (gy0 > gy1) {
+                        std::swap(gx0, gx1);
+                        std::swap(gy0, gy1);
+                    }
+
+                    int y_start = (int)std::ceil(gy0);
+                    int y_end   = (int)std::ceil(gy1) - 1;
+                    if (y_start < 0) y_start = 0;
+                    if (y_end >= height) y_end = height - 1;
+
+                    float dy = gy1 - gy0;
+                    if (dy < 1e-7f) continue;
+
+                    for (int y = y_start; y <= y_end; ++y) {
+                        float t = ((float)y - gy0) / dy;
+                        float gx = gx0 + t * (gx1 - gx0);
+                        int x_int = (int)std::ceil(gx);
+                        if (x_int < 0) x_int = 0;
+                        if (x_int >= width) x_int = width - 1;
+
+                        int w = x_int / 32;
+                        int b = x_int % 32;
+                        parity_words[(size_t)y * row_stride_words + w] ^= (1u << b);
+                    }
                 }
             }
 
+            // Scan parity into mask
             for (int y = 0; y < height; ++y) {
-                float py = min_y + (y + 0.5f) * pitch_mm;
+                const unsigned int* row_p = parity_words.data() + (size_t)y * row_stride_words;
-                std::vector<float> x_inters;
+                uint8_t* row_m = mask + (size_t)y * width;
-
+                uint8_t inside = 0;
-                for (const auto& e : edges) {
+                int x_out = 0;
-                    if ((e.y1 <= py && e.y2 > py) || (e.y2 <= py && e.y1 > py)) {
+                for (int w = 0; w < row_stride_words && x_out < width; ++w) {
-                        float x_int = e.x1 + (py - e.y1) * (e.x2 - e.x1) / (e.y2 - e.y1);
+                    unsigned int word = row_p[w];
-                        x_inters.push_back(x_int);
+                    for (int b = 0; b < 32 && x_out < width; ++b, ++x_out) {
-                    }
+                        inside ^= (uint8_t)((word >> b) & 1u);
-                }
+                        row_m[x_out] = inside;
-
-                std::sort(x_inters.begin(), x_inters.end());
-
-                for (size_t i = 0; i + 1 < x_inters.size(); i += 2) {
-                    int x_start = std::max(0, (int)std::ceil((x_inters[i] - min_x) / pitch_mm - 0.5f));
-                    int x_end = std::min(width, (int)std::floor((x_inters[i + 1] - min_x) / pitch_mm - 0.5f) + 1);
-
-                    for (int x = x_start; x < x_end; ++x) {
-                        mask[(size_t)y * width + x] = 1;
                     }
                 }
             }
         }
-    };
-
-    size_t chunk = (num_layers + num_threads - 1) / num_threads;
-    for (unsigned int t = 0; t < num_threads; ++t) {
-        size_t s = t * chunk;
-        size_t e = std::min(s + chunk, num_layers);
-        if (s < e) {
-            workers.emplace_back(worker_fn, s, e);
     }
 }
-    for (auto& w : workers) w.join();
-}
 
 // -----------------------------------------------------------------------------
 // CPU Implementation: Core Mask Offsetting
@@ -189,43 +252,34 @@ int cpu_offset_compute_masks(
     size_t layer_pixels = (size_t)width * height;
     size_t total_volume = (size_t)num_layers * layer_pixels;
 
-    unsigned int num_threads = std::max(1u, std::thread::hardware_concurrency());
+    #pragma omp parallel
-    std::vector<std::thread> workers;
+    {
-
-    auto worker_fn = [&](size_t start_layer, size_t end_layer) {
         std::vector<float> dist_sq(layer_pixels);
         std::vector<float> g_col_buf(layer_pixels);
         std::vector<float> row_f(width), row_d(width);
         std::vector<int> v_buf(width + 1);
         std::vector<float> z_buf(width + 2);
 
-        for (size_t l = start_layer; l < end_layer; ++l) {
+        #pragma omp for schedule(dynamic, 4)
-            const uint8_t* mask_layer = in_masks + l * layer_pixels;
+        for (int l = 0; l < num_layers; ++l) {
+            const uint8_t* mask_layer = in_masks + (size_t)l * layer_pixels;
 
             // 1. Compute 2D EDT
             edt_2d_layer(mask_layer, width, height, dist_sq.data(),
-                         g_col_buf, row_f, row_d, v_buf, z_buf);
+                         g_col_buf, row_f, row_d, v_buf, z_buf, false);
 
             // 2. Thresholding for each offset distance
             for (int o = 0; o < num_offsets; ++o) {
                 float thresh = offset_distances[o];
                 float thresh_sq_px = (thresh / pitch_mm) * (thresh / pitch_mm);
-                uint8_t* out_layer = out_masks + (size_t)o * total_volume + l * layer_pixels;
+                uint8_t* out_layer = out_masks + (size_t)o * total_volume + (size_t)l * layer_pixels;
 
                 for (size_t i = 0; i < layer_pixels; ++i) {
                     out_layer[i] = (mask_layer[i] > 0 && dist_sq[i] >= thresh_sq_px) ? 1 : 0;
                 }
             }
         }
-    };
-
-    size_t chunk = (num_layers + num_threads - 1) / num_threads;
-    for (unsigned int t = 0; t < num_threads; ++t) {
-        size_t s = t * chunk;
-        size_t e = std::min(s + chunk, (size_t)num_layers);
-        if (s < e) workers.emplace_back(worker_fn, s, e);
     }
-    for (auto& w : workers) w.join();
 
     auto t_end = std::chrono::high_resolution_clock::now();
     double total_ms = std::chrono::duration<double, std::milli>(t_end - t_start).count();
@@ -294,12 +348,11 @@ int cpu_offset_polygon_layers(
     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<std::thread> workers;
-
     int row_stride_words = (cols + 31) / 32;
+    bool subpixel_lerp = config->subpixel_lerp != 0;
 
-    auto worker_fn = [&](size_t start_layer, size_t end_layer) {
+    #pragma omp parallel
+    {
         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);
@@ -307,7 +360,14 @@ int cpu_offset_polygon_layers(
         active_edges_buf.reserve(cols * rows / 16);
         std::vector<uint8_t> visited_buf(max_edges, 0);
 
-        std::vector<float> dist_sq(layer_pixels);
+        std::vector<float> dist_sq_ext(layer_pixels);
+        std::vector<float> dist_sq_int;
+        std::vector<float> layer_sdf;
+        if (subpixel_lerp) {
+            dist_sq_int.resize(layer_pixels);
+            layer_sdf.resize(layer_pixels);
+        }
+
         std::vector<float> g_col_buf(layer_pixels);
         std::vector<float> row_f(cols), row_d(cols);
         std::vector<int> v_buf(cols + 1);
@@ -315,12 +375,26 @@ int cpu_offset_polygon_layers(
 
         std::vector<unsigned int> bit_words((size_t)rows * row_stride_words);
 
-        for (size_t l = start_layer; l < end_layer; ++l) {
+        #pragma omp for schedule(dynamic, 2)
-            const uint8_t* mask_ptr = h_masks.data() + l * layer_pixels;
+        for (int l = 0; l < (int)num_layers; ++l) {
+            const uint8_t* mask_ptr = h_masks.data() + (size_t)l * layer_pixels;
 
-            // Compute exact 2D EDT
+            // Compute exact 2D EDT (exterior)
-            edt_2d_layer(mask_ptr, cols, rows, dist_sq.data(),
+            edt_2d_layer(mask_ptr, cols, rows, dist_sq_ext.data(),
-                         g_col_buf, row_f, row_d, v_buf, z_buf);
+                         g_col_buf, row_f, row_d, v_buf, z_buf, false);
+
+            if (subpixel_lerp) {
+                edt_2d_layer(mask_ptr, cols, rows, dist_sq_int.data(),
+                             g_col_buf, row_f, row_d, v_buf, z_buf, true);
+
+                for (size_t i = 0; i < layer_pixels; ++i) {
+                    if (mask_ptr[i] > 0) {
+                        layer_sdf[i] = sqrtf(dist_sq_ext[i]) * pitch;
+                    } else {
+                        layer_sdf[i] = -sqrtf(dist_sq_int[i]) * pitch;
+                    }
+                }
+            }
 
             res_array[l].z_height = in_layers[l].z_height;
             res_array[l].num_perimeters = P;
@@ -328,13 +402,20 @@ int cpu_offset_polygon_layers(
 
             // Perimeter passes
             for (int o = 0; o < P; ++o) {
+                std::vector<GpuContour> contours;
+
+                if (subpixel_lerp) {
+                    float iso_thresh = resolved_offsets[o];
+                    vectorize_sdf_layer(layer_sdf.data(), cols, rows, bounds_min_x, bounds_min_y, pitch,
+                                        iso_thresh, config->simplify_tolerance_mm, 1,
+                                        next_edge_buf, active_edges_buf, visited_buf, contours);
+                } else {
                     float thresh = resolved_offsets[o];
                     float thresh_sq_px = (thresh / pitch) * (thresh / pitch);
 
-                // Bit-pack thresholded mask
                     for (int y = 0; y < rows; ++y) {
                         const uint8_t* row_m = mask_ptr + (size_t)y * cols;
-                    const float* row_d2 = dist_sq.data() + (size_t)y * cols;
+                        const float* row_d2 = dist_sq_ext.data() + (size_t)y * cols;
                         unsigned int* row_w = bit_words.data() + (size_t)y * row_stride_words;
 
                         for (int w = 0; w < row_stride_words; ++w) {
@@ -352,11 +433,10 @@ int cpu_offset_polygon_layers(
                         }
                     }
 
-                // Marching squares vectorization
-                std::vector<GpuContour> contours;
                     vectorize_bitpack_layer(bit_words.data(), cols, rows, row_stride_words,
                                             bounds_min_x, bounds_min_y, pitch,
                                             next_edge_buf, active_edges_buf, visited_buf, contours);
+                }
 
                 res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o];
                 res_array[l].perimeters[o].num_contours = contours.size();
@@ -369,12 +449,18 @@ int cpu_offset_polygon_layers(
             }
 
             // Infill pass (index P)
+            std::vector<GpuContour> infill_contours;
             float infill_thresh = resolved_offsets[P];
-            float infill_thresh_sq_px = (infill_thresh / pitch) * (infill_thresh / pitch);
 
+            if (subpixel_lerp) {
+                vectorize_sdf_layer(layer_sdf.data(), cols, rows, bounds_min_x, bounds_min_y, pitch,
+                                    infill_thresh, config->simplify_tolerance_mm, 1,
+                                    next_edge_buf, active_edges_buf, visited_buf, infill_contours);
+            } else {
+                float infill_thresh_sq_px = (infill_thresh / pitch) * (infill_thresh / pitch);
                 for (int y = 0; y < rows; ++y) {
                     const uint8_t* row_m = mask_ptr + (size_t)y * cols;
-                const float* row_d2 = dist_sq.data() + (size_t)y * cols;
+                    const float* row_d2 = dist_sq_ext.data() + (size_t)y * cols;
                     unsigned int* row_w = bit_words.data() + (size_t)y * row_stride_words;
 
                     for (int w = 0; w < row_stride_words; ++w) {
@@ -392,12 +478,12 @@ int cpu_offset_polygon_layers(
                     }
                 }
 
-            std::vector<GpuContour> infill_contours;
                 vectorize_bitpack_layer(bit_words.data(), cols, rows, row_stride_words,
                                         bounds_min_x, bounds_min_y, pitch,
                                         next_edge_buf, active_edges_buf, visited_buf, infill_contours);
+            }
 
-            res_array[l].infill.offset_distance_mm = resolved_offsets[P];
+            res_array[l].infill.offset_distance_mm = infill_thresh;
             res_array[l].infill.width = cols;
             res_array[l].infill.height = rows;
             res_array[l].infill.num_contours = infill_contours.size();
@@ -409,16 +495,15 @@ int cpu_offset_polygon_layers(
             }
 
             if (config->export_infill_mask) {
-                uint8_t* infill_mask = (uint8_t*)malloc((size_t)cols * rows);
+                uint8_t* infill_mask = (uint8_t*)malloc(layer_pixels);
-                for (int y = 0; y < rows; ++y) {
+                if (subpixel_lerp) {
-                    const unsigned int* row_words = bit_words.data() + (size_t)y * row_stride_words;
+                    for (size_t i = 0; i < layer_pixels; ++i) {
-                    uint8_t* row_out = infill_mask + (size_t)y * cols;
+                        infill_mask[i] = (layer_sdf[i] >= infill_thresh) ? 1 : 0;
-                    int x_out = 0;
-                    for (int w = 0; w < row_stride_words && x_out < cols; ++w) {
-                        unsigned int word = row_words[w];
-                        for (int b = 0; b < 32 && x_out < cols; ++b, ++x_out) {
-                            row_out[x_out] = (uint8_t)((word >> b) & 1u);
                     }
+                } else {
+                    float infill_thresh_sq_px = (infill_thresh / pitch) * (infill_thresh / pitch);
+                    for (size_t i = 0; i < layer_pixels; ++i) {
+                        infill_mask[i] = (mask_ptr[i] > 0 && dist_sq_ext[i] >= infill_thresh_sq_px) ? 1 : 0;
                     }
                 }
                 res_array[l].infill.mask = infill_mask;
@@ -426,15 +511,7 @@ int cpu_offset_polygon_layers(
                 res_array[l].infill.mask = nullptr;
             }
         }
-    };
-
-    size_t chunk = (num_layers + num_threads - 1) / num_threads;
-    for (unsigned int t = 0; t < num_threads; ++t) {
-        size_t s = t * chunk;
-        size_t e = std::min(s + chunk, num_layers);
-        if (s < e) workers.emplace_back(worker_fn, s, e);
     }
-    for (auto& w : workers) w.join();
 
     auto t_total_end = std::chrono::high_resolution_clock::now();
     double t_total_ms = std::chrono::duration<double, std::milli>(t_total_end - t_total_start).count();
@@ -442,7 +519,7 @@ int cpu_offset_polygon_layers(
     if (stats) {
         memset(stats, 0, sizeof(GpuOffsetStats));
         stats->t_raster_ms = t_raster_ms;
-        stats->t_gpu_jfa_ms = t_total_ms - t_raster_ms; // EDT + Vectorize
+        stats->t_gpu_jfa_ms = t_total_ms - t_raster_ms;
         stats->t_total_pipeline_ms = t_total_ms;
     }
 
@@ -510,6 +587,8 @@ int cpu_offset_slice_mesh_and_offset(
 
     int row_stride_words = (cols + 31) / 32;
     size_t layer_pixels = (size_t)cols * rows;
+    float inv_pitch = 1.0f / pitch;
+    bool subpixel_lerp = config->subpixel_lerp != 0;
 
     // 1. Spatial Triangle Indexing into layers
     auto t_idx_0 = std::chrono::high_resolution_clock::now();
@@ -532,14 +611,19 @@ int cpu_offset_slice_mesh_and_offset(
     auto t_idx_1 = std::chrono::high_resolution_clock::now();
     double t_idx_ms = std::chrono::duration<double, std::milli>(t_idx_1 - t_idx_0).count();
 
-    // 2. Parallel Layer Slicing, Parity Scan, Exact Segment SDF, Thresholding, Vectorization
+    // 2. Parallel Layer Slicing, Parity Scan, Exact SDF / EDT, Thresholding, Vectorization
     GpuLayerResult* res_array = (GpuLayerResult*)malloc(num_layers * sizeof(GpuLayerResult));
     memset(res_array, 0, num_layers * sizeof(GpuLayerResult));
 
-    unsigned int num_threads = std::max(1u, std::thread::hardware_concurrency());
+    const int TILE_SZ = 16;
-    std::vector<std::thread> workers;
+    int num_tiles_x = (cols + TILE_SZ - 1) / TILE_SZ;
+    int num_tiles_y = (rows + TILE_SZ - 1) / TILE_SZ;
+    int total_tiles = num_tiles_x * num_tiles_y;
+    float max_offset = resolved_offsets.empty() ? 2.0f : resolved_offsets.back();
+    float SEARCH_R = std::max(1.5f, max_offset + 2.0f * pitch);
 
-    auto worker_fn = [&](size_t start_layer, size_t end_layer) {
+    #pragma omp parallel
+    {
         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);
@@ -548,19 +632,40 @@ int cpu_offset_slice_mesh_and_offset(
         std::vector<uint8_t> visited_buf(max_edges, 0);
 
         std::vector<unsigned int> parity_words((size_t)rows * row_stride_words, 0);
+        std::vector<unsigned int> bit_words((size_t)rows * row_stride_words, 0);
         std::vector<uint8_t> layer_mask(layer_pixels, 0);
-        std::vector<float> layer_sdf(layer_pixels);
 
-        struct CpuSegment2D {
+        std::vector<float> layer_sdf;
-            float x0, y0, x1, y1;
-        };
         std::vector<CpuSegment2D> layer_segs;
-        layer_segs.reserve(256);
+        std::vector<int> tile_head;
+        std::vector<TileSegNode> seg_nodes;
+        int candidate_buf[1024];
 
-        for (size_t l = start_layer; l < end_layer; ++l) {
+        std::vector<float> dist_sq_ext;
+        std::vector<float> g_col_buf;
+        std::vector<float> row_f, row_d;
+        std::vector<int> v_buf;
+        std::vector<float> z_buf;
+
+        if (subpixel_lerp) {
+            layer_sdf.resize(layer_pixels);
+            layer_segs.reserve(2048);
+            tile_head.assign(total_tiles, -1);
+            seg_nodes.reserve(16384);
+        } else {
+            dist_sq_ext.resize(layer_pixels);
+            g_col_buf.resize(layer_pixels);
+            row_f.resize(cols);
+            row_d.resize(cols);
+            v_buf.resize(cols + 1);
+            z_buf.resize(cols + 2);
+        }
+
+        #pragma omp for schedule(dynamic, 2)
+        for (int l = 0; l < num_layers; ++l) {
             float z0 = z_first + l * layer_height;
             std::fill(parity_words.begin(), parity_words.end(), 0u);
-            layer_segs.clear();
+            if (subpixel_lerp) layer_segs.clear();
 
             // Slicing & Parity Rasterization for Layer l
             const auto& tri_indices = layer_tri_indices[l];
@@ -590,12 +695,25 @@ int cpu_offset_slice_mesh_and_offset(
 
                 if (count != 2) continue;
 
-                layer_segs.push_back({ pts_x[0], pts_y[0], pts_x[1], pts_y[1] });
+                float x0 = pts_x[0], y0 = pts_y[0], x1 = pts_x[1], y1 = pts_y[1];
+                if (subpixel_lerp) {
+                    float vx = x1 - x0;
+                    float vy = y1 - y0;
+                    float len_sq = vx * vx + vy * vy;
+                    float inv_len_sq = (len_sq > 1e-10f) ? (1.0f / len_sq) : 0.0f;
+                    layer_segs.push_back({
+                        x0, y0, x1, y1,
+                        vx, vy,
+                        inv_len_sq,
+                        std::min(x0, x1), std::max(x0, x1),
+                        std::min(y0, y1), std::max(y0, y1)
+                    });
+                }
 
-                float gx0 = (pts_x[0] - bounds_min_x) / pitch;
+                float gx0 = (x0 - bounds_min_x) * inv_pitch;
-                float gy0 = (pts_y[0] - bounds_min_y) / pitch;
+                float gy0 = (y0 - bounds_min_y) * inv_pitch;
-                float gx1 = (pts_x[1] - bounds_min_x) / pitch;
+                float gx1 = (x1 - bounds_min_x) * inv_pitch;
-                float gy1 = (pts_y[1] - bounds_min_y) / pitch;
+                float gy1 = (y1 - bounds_min_y) * inv_pitch;
 
                 if (gy0 == gy1) continue;
 
@@ -627,6 +745,7 @@ int cpu_offset_slice_mesh_and_offset(
             }
 
             // Scan Parity Words into Binary Mask
+            bool has_inside = false;
             for (int y = 0; y < rows; ++y) {
                 const unsigned int* row_p = parity_words.data() + (size_t)y * row_stride_words;
                 uint8_t* row_m = layer_mask.data() + (size_t)y * cols;
@@ -637,39 +756,50 @@ int cpu_offset_slice_mesh_and_offset(
                     for (int b = 0; b < 32 && x_out < cols; ++b, ++x_out) {
                         inside ^= (uint8_t)((word >> b) & 1u);
                         row_m[x_out] = inside;
+                        if (inside) has_inside = true;
                     }
                 }
             }
 
-            // Fast 2D Tile Spatial Binning for analytical segment distance on CPU
+            res_array[l].z_height = z0;
-            const int TILE_SZ = 16;
+            res_array[l].num_perimeters = P;
-            int num_tiles_x = (cols + TILE_SZ - 1) / TILE_SZ;
+            res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter));
-            int num_tiles_y = (rows + TILE_SZ - 1) / TILE_SZ;
+
-            std::vector<std::vector<int>> tile_segs(num_tiles_x * num_tiles_y);
+            if (!has_inside || (subpixel_lerp && layer_segs.empty())) {
+                for (int o = 0; o < P; ++o) {
+                    res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o];
+                    res_array[l].perimeters[o].num_contours = 0;
+                    res_array[l].perimeters[o].contours = nullptr;
+                }
+                res_array[l].infill.offset_distance_mm = resolved_offsets[P];
+                res_array[l].infill.width = cols;
+                res_array[l].infill.height = rows;
+                res_array[l].infill.num_contours = 0;
+                res_array[l].infill.contours = nullptr;
+                res_array[l].infill.mask = nullptr;
+                continue;
+            }
+
+            if (subpixel_lerp) {
+                // Exact analytical segment distance field with flat linked-list tile binning
+                std::fill(tile_head.begin(), tile_head.end(), -1);
+                seg_nodes.clear();
 
-            const float SEARCH_R = 4.0f;
-            float inv_pitch = 1.0f / pitch;
                 int num_segs = (int)layer_segs.size();
-
                 for (int s = 0; s < num_segs; ++s) {
-                float sx0 = layer_segs[s].x0;
+                    const auto& seg = layer_segs[s];
-                float sy0 = layer_segs[s].y0;
+                    int tx0 = std::max(0, (int)((seg.min_x - SEARCH_R - bounds_min_x) * inv_pitch) / TILE_SZ);
-                float sx1 = layer_segs[s].x1;
+                    int tx1 = std::min(num_tiles_x - 1, (int)((seg.max_x + SEARCH_R - bounds_min_x) * inv_pitch) / TILE_SZ);
-                float sy1 = layer_segs[s].y1;
+                    int ty0 = std::max(0, (int)((seg.min_y - SEARCH_R - bounds_min_y) * inv_pitch) / TILE_SZ);
-
+                    int ty1 = std::min(num_tiles_y - 1, (int)((seg.max_y + SEARCH_R - bounds_min_y) * inv_pitch) / TILE_SZ);
-                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) {
+                        int row_tile_idx = ty * num_tiles_x;
                         for (int tx = tx0; tx <= tx1; ++tx) {
-                        tile_segs[ty * num_tiles_x + tx].push_back(s);
+                            int t_idx = row_tile_idx + tx;
+                            int node = (int)seg_nodes.size();
+                            seg_nodes.push_back({ s, tile_head[t_idx] });
+                            tile_head[t_idx] = node;
                         }
                     }
                 }
@@ -677,12 +807,15 @@ int cpu_offset_slice_mesh_and_offset(
                 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);
+                    int row_tile_idx = ty * num_tiles_x;
+
                     for (int tx = 0; tx < num_tiles_x; ++tx) {
+                        int t_idx = row_tile_idx + 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()) {
+                        int head = tile_head[t_idx];
+                        if (head == -1) {
                             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];
@@ -692,27 +825,35 @@ int cpu_offset_slice_mesh_and_offset(
                             continue;
                         }
 
+                        int num_cand = 0;
+                        int cur = head;
+                        while (cur != -1 && num_cand < 1024) {
+                            candidate_buf[num_cand++] = seg_nodes[cur].seg_idx;
+                            cur = seg_nodes[cur].next;
+                        }
+
                         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;
+                                float min_d2 = 1e10f;
-                            for (int s_idx : candidates) {
+
-                                float sx0 = layer_segs[s_idx].x0;
+                                for (int k = 0; k < num_cand; ++k) {
-                                float sy0 = layer_segs[s_idx].y0;
+                                    int s_idx = candidate_buf[k];
-                                float sx1 = layer_segs[s_idx].x1;
+                                    const auto& seg = layer_segs[s_idx];
-                                float sy1 = layer_segs[s_idx].y1;
+                                    float wx = px - seg.x0;
-                                float dx = sx1 - sx0;
+                                    float wy = py - seg.y0;
-                                float dy = sy1 - sy0;
+                                    float t = (wx * seg.vx + wy * seg.vy) * seg.inv_len_sq;
-                                float l2 = dx * dx + dy * dy;
+                                    t = (t < 0.0f) ? 0.0f : ((t > 1.0f) ? 1.0f : t);
-                                float t = (l2 > 1e-12f) ? ((px - sx0) * dx + (py - sy0) * dy) / l2 : 0.0f;
+                                    float qx = seg.x0 + t * seg.vx;
-                                t = std::max(0.0f, std::min(1.0f, t));
+                                    float qy = seg.y0 + t * seg.vy;
-                                float qx = sx0 + t * dx;
+                                    float dx = px - qx;
-                                float qy = sy0 + t * dy;
+                                    float dy = py - qy;
-                                float d2 = (px - qx) * (px - qx) + (py - qy) * (py - qy);
+                                    float d2 = dx * dx + dy * dy;
                                     if (d2 < min_d2) min_d2 = d2;
                                 }
-                            float d = (min_d2 < 1e19f) ? sqrtf(min_d2) : SEARCH_R;
+
+                                float d = (min_d2 < 1e9f) ? std::sqrt(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;
                             }
@@ -720,17 +861,12 @@ int cpu_offset_slice_mesh_and_offset(
                     }
                 }
 
-            res_array[l].z_height = z0;
+                // Perimeters with continuous SDF
-            res_array[l].num_perimeters = P;
-            res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter));
-
-            // Perimeters with continuous SDF and sub-pixel LERP
                 for (int o = 0; o < P; ++o) {
                     float iso_thresh = resolved_offsets[o];
-
                     std::vector<GpuContour> contours;
                     vectorize_sdf_layer(layer_sdf.data(), cols, rows, bounds_min_x, bounds_min_y, pitch,
-                                    iso_thresh, config->simplify_tolerance_mm, config->subpixel_lerp,
+                                        iso_thresh, config->simplify_tolerance_mm, 1,
                                         next_edge_buf, active_edges_buf, visited_buf, contours);
 
                     res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o];
@@ -744,13 +880,13 @@ int cpu_offset_slice_mesh_and_offset(
                 }
 
                 // Infill pass
-            float infill_iso_thresh = resolved_offsets[P];
+                float infill_thresh = resolved_offsets[P];
                 std::vector<GpuContour> infill_contours;
                 vectorize_sdf_layer(layer_sdf.data(), cols, rows, bounds_min_x, bounds_min_y, pitch,
-                                infill_iso_thresh, config->simplify_tolerance_mm, config->subpixel_lerp,
+                                    infill_thresh, config->simplify_tolerance_mm, 1,
                                     next_edge_buf, active_edges_buf, visited_buf, infill_contours);
 
-            res_array[l].infill.offset_distance_mm = infill_iso_thresh;
+                res_array[l].infill.offset_distance_mm = infill_thresh;
                 res_array[l].infill.width = cols;
                 res_array[l].infill.height = rows;
                 res_array[l].infill.num_contours = infill_contours.size();
@@ -764,22 +900,109 @@ int cpu_offset_slice_mesh_and_offset(
                 if (config->export_infill_mask) {
                     uint8_t* infill_mask = (uint8_t*)malloc(layer_pixels);
                     for (size_t i = 0; i < layer_pixels; ++i) {
-                    infill_mask[i] = (layer_sdf[i] >= infill_iso_thresh) ? 1 : 0;
+                        infill_mask[i] = (layer_sdf[i] >= infill_thresh) ? 1 : 0;
+                    }
+                    res_array[l].infill.mask = infill_mask;
+                } else {
+                    res_array[l].infill.mask = nullptr;
+                }
+
+            } else {
+                // Discrete Marching Squares (No LERP) via Exact 2D EDT & Bitpacking
+                edt_2d_layer(layer_mask.data(), cols, rows, dist_sq_ext.data(),
+                             g_col_buf, row_f, row_d, v_buf, z_buf, false);
+
+                for (int o = 0; o < P; ++o) {
+                    float thresh = resolved_offsets[o];
+                    float thresh_sq_px = (thresh / pitch) * (thresh / pitch);
+
+                    for (int y = 0; y < rows; ++y) {
+                        const uint8_t* row_m = layer_mask.data() + (size_t)y * cols;
+                        const float* row_d2 = dist_sq_ext.data() + (size_t)y * cols;
+                        unsigned int* row_w = bit_words.data() + (size_t)y * row_stride_words;
+
+                        for (int w = 0; w < row_stride_words; ++w) {
+                            unsigned int packed = 0;
+                            int x_base = w * 32;
+                            for (int b = 0; b < 32; ++b) {
+                                int x = x_base + b;
+                                if (x < cols) {
+                                    if (row_m[x] > 0 && row_d2[x] >= thresh_sq_px) {
+                                        packed |= (1u << b);
+                                    }
+                                }
+                            }
+                            row_w[w] = packed;
+                        }
+                    }
+
+                    std::vector<GpuContour> contours;
+                    vectorize_bitpack_layer(bit_words.data(), cols, rows, row_stride_words,
+                                            bounds_min_x, bounds_min_y, pitch,
+                                            next_edge_buf, active_edges_buf, visited_buf, contours);
+
+                    res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o];
+                    res_array[l].perimeters[o].num_contours = contours.size();
+                    if (!contours.empty()) {
+                        res_array[l].perimeters[o].contours = (GpuContour*)malloc(contours.size() * sizeof(GpuContour));
+                        memcpy(res_array[l].perimeters[o].contours, contours.data(), contours.size() * sizeof(GpuContour));
+                    } else {
+                        res_array[l].perimeters[o].contours = nullptr;
+                    }
+                }
+
+                // Infill pass (index P)
+                float infill_thresh = resolved_offsets[P];
+                float infill_thresh_sq_px = (infill_thresh / pitch) * (infill_thresh / pitch);
+
+                for (int y = 0; y < rows; ++y) {
+                    const uint8_t* row_m = layer_mask.data() + (size_t)y * cols;
+                    const float* row_d2 = dist_sq_ext.data() + (size_t)y * cols;
+                    unsigned int* row_w = bit_words.data() + (size_t)y * row_stride_words;
+
+                    for (int w = 0; w < row_stride_words; ++w) {
+                        unsigned int packed = 0;
+                        int x_base = w * 32;
+                        for (int b = 0; b < 32; ++b) {
+                            int x = x_base + b;
+                            if (x < cols) {
+                                if (row_m[x] > 0 && row_d2[x] >= infill_thresh_sq_px) {
+                                    packed |= (1u << b);
+                                }
+                            }
+                        }
+                        row_w[w] = packed;
+                    }
+                }
+
+                std::vector<GpuContour> infill_contours;
+                vectorize_bitpack_layer(bit_words.data(), cols, rows, row_stride_words,
+                                        bounds_min_x, bounds_min_y, pitch,
+                                        next_edge_buf, active_edges_buf, visited_buf, infill_contours);
+
+                res_array[l].infill.offset_distance_mm = infill_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(layer_pixels);
+                    for (size_t i = 0; i < layer_pixels; ++i) {
+                        infill_mask[i] = (layer_mask[i] > 0 && dist_sq_ext[i] >= infill_thresh_sq_px) ? 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, (size_t)num_layers);
-        if (s < e) workers.emplace_back(worker_fn, s, e);
         }
-    for (auto& w : workers) w.join();
+    }
 
     auto t_total_end = std::chrono::high_resolution_clock::now();
     double t_total_ms = std::chrono::duration<double, std::milli>(t_total_end - t_total_start).count();

src/vectorizer.h

@@ -7,6 +7,7 @@
 #include <cstring>
 #include <cstdlib>
 #include <cstdint>
+#include <algorithm>
 
 inline int get_h_edge_id(int x, int y, int stride_w) {
     return 2 * (y * stride_w + x) + 0;
@@ -82,13 +83,16 @@ inline GpuPoint2D get_sdf_edge_point(
 static inline void rdp_simplify_recursive(
     const std::vector<GpuPoint2D>& pts,
     size_t i, size_t j,
+    size_t n,
     float eps_sq,
-    std::vector<bool>& keep
+    std::vector<uint8_t>& keep
 ) {
     if (j <= i + 1) return;
 
-    float p1x = pts[i].x, p1y = pts[i].y;
+    const GpuPoint2D& p1 = pts[i % n];
-    float p2x = pts[j].x, p2y = pts[j].y;
+    const GpuPoint2D& p2 = pts[j % n];
+    float p1x = p1.x, p1y = p1.y;
+    float p2x = p2.x, p2y = p2.y;
     float dx = p2x - p1x;
     float dy = p2y - p1y;
     float len_sq = dx * dx + dy * dy;
@@ -97,7 +101,8 @@ static inline void rdp_simplify_recursive(
     size_t max_idx = i;
 
     for (size_t k = i + 1; k < j; ++k) {
-        float px = pts[k].x, py = pts[k].y;
+        const GpuPoint2D& pk = pts[k % n];
+        float px = pk.x, py = pk.y;
         float dist_sq = 0.0f;
 
         if (len_sq < 1e-10f) {
@@ -116,9 +121,9 @@ static inline void rdp_simplify_recursive(
     }
 
     if (max_dist_sq > eps_sq) {
-        keep[max_idx] = true;
+        keep[max_idx] = 1;
-        rdp_simplify_recursive(pts, i, max_idx, eps_sq, keep);
+        rdp_simplify_recursive(pts, i, max_idx, n, eps_sq, keep);
-        rdp_simplify_recursive(pts, max_idx, j, eps_sq, keep);
+        rdp_simplify_recursive(pts, max_idx, j, n, eps_sq, keep);
     }
 }
 
@@ -144,16 +149,13 @@ static inline std::vector<GpuPoint2D> simplify_contour_rdp(
         }
     }
 
-    std::vector<GpuPoint2D> closed_pts = loop_pts;
+    std::vector<uint8_t> keep(n + 1, 0);
-    closed_pts.push_back(loop_pts[0]);
+    keep[0] = 1;
+    keep[split_idx] = 1;
+    keep[n] = 1;
 
-    std::vector<bool> keep(n + 1, false);
+    rdp_simplify_recursive(loop_pts, 0, split_idx, n, eps_sq, keep);
-    keep[0] = true;
+    rdp_simplify_recursive(loop_pts, split_idx, n, n, eps_sq, keep);
-    keep[split_idx] = true;
-    keep[n] = true;
-
-    rdp_simplify_recursive(closed_pts, 0, split_idx, eps_sq, keep);
-    rdp_simplify_recursive(closed_pts, split_idx, n, eps_sq, keep);
 
     std::vector<GpuPoint2D> simplified;
     simplified.reserve(n);
@@ -164,9 +166,9 @@ static inline std::vector<GpuPoint2D> simplify_contour_rdp(
 }
 
 // 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) {
+static inline void align_contour_seam_sharpest_cpp_inplace(std::vector<GpuPoint2D>& loop) {
     size_t M = loop.size();
-    if (M < 3) return loop;
+    if (M < 3) return;
 
     // 1. Calculate signed area
     double signed_area = 0.0;
@@ -175,7 +177,7 @@ static inline std::vector<GpuPoint2D> align_contour_seam_sharpest_cpp(const std:
         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;
+    if (std::abs(signed_area) < 1e-9) return;
 
     // 2. Find sharpest concave angle or fallback to sharpest angle
     size_t best_concave_idx = 0;
@@ -218,12 +220,14 @@ static inline std::vector<GpuPoint2D> align_contour_seam_sharpest_cpp(const std:
     }
 
     size_t best_idx = found_concave ? best_concave_idx : best_any_idx;
-    if (best_idx == 0) return loop;
+    if (best_idx != 0) {
+        std::rotate(loop.begin(), loop.begin() + best_idx, loop.end());
+    }
+}
 
-    std::vector<GpuPoint2D> aligned;
+static inline std::vector<GpuPoint2D> align_contour_seam_sharpest_cpp(const std::vector<GpuPoint2D>& loop) {
-    aligned.reserve(M);
+    std::vector<GpuPoint2D> aligned = loop;
-    for (size_t i = best_idx; i < M; ++i) aligned.push_back(loop[i]);
+    align_contour_seam_sharpest_cpp_inplace(aligned);
-    for (size_t i = 0; i < best_idx; ++i) aligned.push_back(loop[i]);
     return aligned;
 }