diff --git a/CMakeLists.txt b/CMakeLists.txt index 72b4664..e080eb3 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/python/cudaslicer/bindings.cpp b/python/cudaslicer/bindings.cpp index 12e17cf..dfaa363 100644 --- a/python/cudaslicer/bindings.cpp +++ b/python/cudaslicer/bindings.cpp @@ -79,10 +79,14 @@ py::array_t py_offset_masks( } GpuOffsetStats stats = {0}; - int ret = gpu_offset_compute_masks( - ptr_in, width, height, num_layers, pitch_mm, - resolved_offsets.data(), num_offsets, ptr_out, &stats - ); + 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( - c_layers.data(), num_layers, - min_x, min_y, max_x, max_y, - &config, &results, &stats - ); + 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( - 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 - ); + 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."); diff --git a/src/cpu_offset.cpp b/src/cpu_offset.cpp index ff6fd82..1eb66ee 100644 --- a/src/cpu_offset.cpp +++ b/src/cpu_offset.cpp @@ -6,7 +6,7 @@ #include #include #include -#include +#include // ----------------------------------------------------------------------------- // 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); - col_g[y] = (last >= 0) ? (dy * dy) : INF; + int dy = y - last; + 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); - float d2 = (last < height) ? (dy * dy) : INF; + int dy = last - y; + 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,11 +104,18 @@ static inline void edt_2d_layer( std::vector& row_f, std::vector& row_d, std::vector& v_buf, - std::vector& z_buf + std::vector& z_buf, + bool invert = false ) { // Pass 1: Columns - for (int x = 0; x < width; ++x) { - edt_1d_binary_col(mask + x, height, width, g_col_buf.data() + (size_t)x * height); + 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 @@ -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 { - float x1, y1; - float x2, y2; +struct CpuSegment2D { + float x0, y0, x1, y1; + 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,59 +155,79 @@ 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()); - std::vector workers; + int row_stride_words = (width + 31) / 32; + float inv_pitch = 1.0f / pitch_mm; - auto worker_fn = [&](size_t start_layer, size_t end_layer) { - for (size_t l = start_layer; l < end_layer; ++l) { - uint8_t* mask = out_masks + l * ((size_t)width * height); + #pragma omp parallel + { + std::vector 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 edges; for (size_t c = 0; c < layer.num_contours; ++c) { const GpuContour& cnt = layer.contours[c]; if (cnt.num_points < 3) continue; + for (size_t p = 0; p < cnt.num_points; ++p) { - size_t p_next = (p + 1) % cnt.num_points; - edges.push_back({ cnt.points[p].x, cnt.points[p].y, - cnt.points[p_next].x, cnt.points[p_next].y }); + size_t p_next = (p + 1 == cnt.num_points) ? 0 : p + 1; + float px0 = cnt.points[p].x; + 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; - std::vector x_inters; - - for (const auto& e : edges) { - if ((e.y1 <= py && e.y2 > py) || (e.y2 <= py && e.y1 > py)) { - float x_int = e.x1 + (py - e.y1) * (e.x2 - e.x1) / (e.y2 - e.y1); - x_inters.push_back(x_int); - } - } - - std::sort(x_inters.begin(), x_inters.end()); - - for (size_t i = 0; i + 1 < x_inters.size(); i += 2) { - int x_start = std::max(0, (int)std::ceil((x_inters[i] - min_x) / pitch_mm - 0.5f)); - int x_end = std::min(width, (int)std::floor((x_inters[i + 1] - min_x) / pitch_mm - 0.5f) + 1); - - for (int x = x_start; x < x_end; ++x) { - mask[(size_t)y * width + x] = 1; + const unsigned int* row_p = parity_words.data() + (size_t)y * row_stride_words; + uint8_t* row_m = mask + (size_t)y * width; + uint8_t inside = 0; + int x_out = 0; + for (int w = 0; w < row_stride_words && x_out < width; ++w) { + unsigned int word = row_p[w]; + for (int b = 0; b < 32 && x_out < width; ++b, ++x_out) { + inside ^= (uint8_t)((word >> b) & 1u); + row_m[x_out] = inside; } } } } - }; - - 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(); } // ----------------------------------------------------------------------------- @@ -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()); - std::vector workers; - - auto worker_fn = [&](size_t start_layer, size_t end_layer) { + #pragma omp parallel + { std::vector dist_sq(layer_pixels); std::vector g_col_buf(layer_pixels); std::vector row_f(width), row_d(width); std::vector v_buf(width + 1); std::vector z_buf(width + 2); - for (size_t l = start_layer; l < end_layer; ++l) { - const uint8_t* mask_layer = in_masks + l * layer_pixels; + #pragma omp for schedule(dynamic, 4) + 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(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 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 next_edge_buf(max_edges, -1); @@ -307,7 +360,14 @@ int cpu_offset_polygon_layers( active_edges_buf.reserve(cols * rows / 16); std::vector visited_buf(max_edges, 0); - std::vector dist_sq(layer_pixels); + std::vector dist_sq_ext(layer_pixels); + std::vector dist_sq_int; + std::vector layer_sdf; + if (subpixel_lerp) { + dist_sq_int.resize(layer_pixels); + layer_sdf.resize(layer_pixels); + } + std::vector g_col_buf(layer_pixels); std::vector row_f(cols), row_d(cols); std::vector v_buf(cols + 1); @@ -315,12 +375,26 @@ int cpu_offset_polygon_layers( std::vector bit_words((size_t)rows * row_stride_words); - for (size_t l = start_layer; l < end_layer; ++l) { - const uint8_t* mask_ptr = h_masks.data() + l * layer_pixels; + #pragma omp for schedule(dynamic, 2) + 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 - edt_2d_layer(mask_ptr, cols, rows, dist_sq.data(), - g_col_buf, row_f, row_d, v_buf, z_buf); + // Compute exact 2D EDT (exterior) + edt_2d_layer(mask_ptr, cols, rows, dist_sq_ext.data(), + 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,35 +402,41 @@ int cpu_offset_polygon_layers( // Perimeter passes for (int o = 0; o < P; ++o) { - float thresh = resolved_offsets[o]; - float thresh_sq_px = (thresh / pitch) * (thresh / pitch); + std::vector contours; - // Bit-pack thresholded mask - for (int y = 0; y < rows; ++y) { - const uint8_t* row_m = mask_ptr + (size_t)y * cols; - const float* row_d2 = dist_sq.data() + (size_t)y * cols; - unsigned int* row_w = bit_words.data() + (size_t)y * row_stride_words; + 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); - 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); + for (int y = 0; y < rows; ++y) { + const uint8_t* row_m = mask_ptr + (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; } - row_w[w] = packed; } - } - // Marching squares vectorization - std::vector contours; - vectorize_bitpack_layer(bit_words.data(), cols, rows, row_stride_words, - bounds_min_x, bounds_min_y, pitch, - next_edge_buf, active_edges_buf, visited_buf, contours); + 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,35 +449,41 @@ int cpu_offset_polygon_layers( } // Infill pass (index P) + std::vector infill_contours; float infill_thresh = resolved_offsets[P]; - float infill_thresh_sq_px = (infill_thresh / pitch) * (infill_thresh / pitch); - for (int y = 0; y < rows; ++y) { - const uint8_t* row_m = mask_ptr + (size_t)y * cols; - const float* row_d2 = dist_sq.data() + (size_t)y * cols; - unsigned int* row_w = bit_words.data() + (size_t)y * row_stride_words; + 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_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); + 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; } - row_w[w] = packed; } + + 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); } - std::vector infill_contours; - vectorize_bitpack_layer(bit_words.data(), cols, rows, row_stride_words, - bounds_min_x, bounds_min_y, pitch, - next_edge_buf, active_edges_buf, visited_buf, infill_contours); - - res_array[l].infill.offset_distance_mm = resolved_offsets[P]; + res_array[l].infill.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); - for (int y = 0; y < rows; ++y) { - const unsigned int* row_words = bit_words.data() + (size_t)y * row_stride_words; - uint8_t* row_out = infill_mask + (size_t)y * cols; - int x_out = 0; - for (int w = 0; w < row_stride_words && x_out < cols; ++w) { - unsigned int word = row_words[w]; - for (int b = 0; b < 32 && x_out < cols; ++b, ++x_out) { - row_out[x_out] = (uint8_t)((word >> b) & 1u); - } + uint8_t* infill_mask = (uint8_t*)malloc(layer_pixels); + if (subpixel_lerp) { + for (size_t i = 0; i < layer_pixels; ++i) { + infill_mask[i] = (layer_sdf[i] >= infill_thresh) ? 1 : 0; + } + } 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(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(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()); - std::vector workers; + 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; + 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 next_edge_buf(max_edges, -1); @@ -548,19 +632,40 @@ int cpu_offset_slice_mesh_and_offset( std::vector visited_buf(max_edges, 0); std::vector parity_words((size_t)rows * row_stride_words, 0); + std::vector bit_words((size_t)rows * row_stride_words, 0); std::vector layer_mask(layer_pixels, 0); - std::vector layer_sdf(layer_pixels); - struct CpuSegment2D { - float x0, y0, x1, y1; - }; + std::vector layer_sdf; std::vector layer_segs; - layer_segs.reserve(256); + std::vector tile_head; + std::vector seg_nodes; + int candidate_buf[1024]; - for (size_t l = start_layer; l < end_layer; ++l) { + std::vector dist_sq_ext; + std::vector g_col_buf; + std::vector row_f, row_d; + std::vector v_buf; + std::vector 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 gy0 = (pts_y[0] - bounds_min_y) / pitch; - float gx1 = (pts_x[1] - bounds_min_x) / pitch; - float gy1 = (pts_y[1] - bounds_min_y) / pitch; + float gx0 = (x0 - bounds_min_x) * inv_pitch; + float gy0 = (y0 - bounds_min_y) * inv_pitch; + float gx1 = (x1 - bounds_min_x) * inv_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,85 +756,7 @@ 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; - } - } - } - - // Fast 2D Tile Spatial Binning for analytical segment distance on CPU - const int TILE_SZ = 16; - int num_tiles_x = (cols + TILE_SZ - 1) / TILE_SZ; - int num_tiles_y = (rows + TILE_SZ - 1) / TILE_SZ; - std::vector> tile_segs(num_tiles_x * num_tiles_y); - - const float SEARCH_R = 4.0f; - float inv_pitch = 1.0f / pitch; - int num_segs = (int)layer_segs.size(); - - for (int s = 0; s < num_segs; ++s) { - float sx0 = layer_segs[s].x0; - float sy0 = layer_segs[s].y0; - float sx1 = layer_segs[s].x1; - float sy1 = layer_segs[s].y1; - - float min_x = std::min(sx0, sx1) - SEARCH_R; - float max_x = std::max(sx0, sx1) + SEARCH_R; - float min_y = std::min(sy0, sy1) - SEARCH_R; - float max_y = std::max(sy0, sy1) + SEARCH_R; - - int tx0 = std::max(0, (int)((min_x - bounds_min_x) * inv_pitch) / TILE_SZ); - int tx1 = std::min(num_tiles_x - 1, (int)((max_x - bounds_min_x) * inv_pitch) / TILE_SZ); - int ty0 = std::max(0, (int)((min_y - bounds_min_y) * inv_pitch) / TILE_SZ); - int ty1 = std::min(num_tiles_y - 1, (int)((max_y - bounds_min_y) * inv_pitch) / TILE_SZ); - - for (int ty = ty0; ty <= ty1; ++ty) { - for (int tx = tx0; tx <= tx1; ++tx) { - tile_segs[ty * num_tiles_x + tx].push_back(s); - } - } - } - - for (int ty = 0; ty < num_tiles_y; ++ty) { - int y_start = ty * TILE_SZ; - int y_end = std::min(rows, y_start + TILE_SZ); - for (int tx = 0; tx < num_tiles_x; ++tx) { - int x_start = tx * TILE_SZ; - int x_end = std::min(cols, x_start + TILE_SZ); - const auto& candidates = tile_segs[ty * num_tiles_x + tx]; - - if (candidates.empty()) { - for (int y = y_start; y < y_end; ++y) { - for (int x = x_start; x < x_end; ++x) { - uint8_t m = layer_mask[(size_t)y * cols + x]; - layer_sdf[(size_t)y * cols + x] = (m > 0) ? SEARCH_R : -SEARCH_R; - } - } - continue; - } - - for (int y = y_start; y < y_end; ++y) { - float py = bounds_min_y + (float)y * pitch; - for (int x = x_start; x < x_end; ++x) { - float px = bounds_min_x + (float)x * pitch; - float min_d2 = 1e20f; - for (int s_idx : candidates) { - float sx0 = layer_segs[s_idx].x0; - float sy0 = layer_segs[s_idx].y0; - float sx1 = layer_segs[s_idx].x1; - float sy1 = layer_segs[s_idx].y1; - float dx = sx1 - sx0; - float dy = sy1 - sy0; - float l2 = dx * dx + dy * dy; - float t = (l2 > 1e-12f) ? ((px - sx0) * dx + (py - sy0) * dy) / l2 : 0.0f; - t = std::max(0.0f, std::min(1.0f, t)); - float qx = sx0 + t * dx; - float qy = sy0 + t * dy; - float d2 = (px - qx) * (px - qx) + (py - qy) * (py - qy); - if (d2 < min_d2) min_d2 = d2; - } - float d = (min_d2 < 1e19f) ? sqrtf(min_d2) : SEARCH_R; - uint8_t m = layer_mask[(size_t)y * cols + x]; - layer_sdf[(size_t)y * cols + x] = (m > 0) ? d : -d; - } + if (inside) has_inside = true; } } } @@ -724,62 +765,244 @@ int cpu_offset_slice_mesh_and_offset( 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 contours; - vectorize_sdf_layer(layer_sdf.data(), cols, rows, bounds_min_x, bounds_min_y, pitch, - iso_thresh, config->simplify_tolerance_mm, config->subpixel_lerp, - next_edge_buf, active_edges_buf, visited_buf, contours); - - res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o]; - res_array[l].perimeters[o].num_contours = contours.size(); - 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 { + 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; } - } - - // Infill pass - float infill_iso_thresh = resolved_offsets[P]; - std::vector infill_contours; - vectorize_sdf_layer(layer_sdf.data(), cols, rows, bounds_min_x, bounds_min_y, pitch, - infill_iso_thresh, config->simplify_tolerance_mm, config->subpixel_lerp, - next_edge_buf, active_edges_buf, visited_buf, infill_contours); - - res_array[l].infill.offset_distance_mm = 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.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 (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; + 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(); + + int num_segs = (int)layer_segs.size(); + for (int s = 0; s < num_segs; ++s) { + const auto& seg = layer_segs[s]; + int tx0 = std::max(0, (int)((seg.min_x - SEARCH_R - bounds_min_x) * inv_pitch) / TILE_SZ); + int tx1 = std::min(num_tiles_x - 1, (int)((seg.max_x + SEARCH_R - bounds_min_x) * inv_pitch) / TILE_SZ); + 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); + + for (int ty = ty0; ty <= ty1; ++ty) { + int row_tile_idx = ty * num_tiles_x; + for (int tx = tx0; tx <= tx1; ++tx) { + 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; + } + } } - res_array[l].infill.mask = infill_mask; + + 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); + + 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]; + layer_sdf[(size_t)y * cols + x] = (m > 0) ? SEARCH_R : -SEARCH_R; + } + } + 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 = 1e10f; + + for (int k = 0; k < num_cand; ++k) { + int s_idx = candidate_buf[k]; + const auto& seg = layer_segs[s_idx]; + float wx = px - seg.x0; + float wy = py - seg.y0; + float t = (wx * seg.vx + wy * seg.vy) * seg.inv_len_sq; + t = (t < 0.0f) ? 0.0f : ((t > 1.0f) ? 1.0f : t); + float qx = seg.x0 + t * seg.vx; + float qy = seg.y0 + t * seg.vy; + float dx = px - qx; + float dy = py - qy; + float d2 = dx * dx + dy * dy; + if (d2 < min_d2) min_d2 = d2; + } + + 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; + } + } + } + } + + // Perimeters with continuous SDF + for (int o = 0; o < P; ++o) { + float iso_thresh = resolved_offsets[o]; + std::vector contours; + 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); + + res_array[l].perimeters[o].offset_distance_mm = resolved_offsets[o]; + res_array[l].perimeters[o].num_contours = contours.size(); + if (!contours.empty()) { + res_array[l].perimeters[o].contours = (GpuContour*)malloc(contours.size() * sizeof(GpuContour)); + memcpy(res_array[l].perimeters[o].contours, contours.data(), contours.size() * sizeof(GpuContour)); + } else { + res_array[l].perimeters[o].contours = nullptr; + } + } + + // Infill pass + float infill_thresh = resolved_offsets[P]; + std::vector infill_contours; + 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); + + 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_sdf[i] >= infill_thresh) ? 1 : 0; + } + res_array[l].infill.mask = infill_mask; + } else { + res_array[l].infill.mask = nullptr; + } + } else { - res_array[l].infill.mask = nullptr; + // 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 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 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(t_total_end - t_total_start).count(); diff --git a/src/vectorizer.h b/src/vectorizer.h index 312851d..3344de2 100644 --- a/src/vectorizer.h +++ b/src/vectorizer.h @@ -7,6 +7,7 @@ #include #include #include +#include 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& pts, size_t i, size_t j, + size_t n, float eps_sq, - std::vector& keep + std::vector& keep ) { if (j <= i + 1) return; - float p1x = pts[i].x, p1y = pts[i].y; - float p2x = pts[j].x, p2y = pts[j].y; + const GpuPoint2D& p1 = pts[i % n]; + 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; - rdp_simplify_recursive(pts, i, max_idx, eps_sq, keep); - rdp_simplify_recursive(pts, max_idx, j, eps_sq, keep); + keep[max_idx] = 1; + rdp_simplify_recursive(pts, i, max_idx, n, eps_sq, keep); + rdp_simplify_recursive(pts, max_idx, j, n, eps_sq, keep); } } @@ -144,16 +149,13 @@ static inline std::vector simplify_contour_rdp( } } - std::vector closed_pts = loop_pts; - closed_pts.push_back(loop_pts[0]); + std::vector keep(n + 1, 0); + keep[0] = 1; + keep[split_idx] = 1; + keep[n] = 1; - std::vector keep(n + 1, false); - keep[0] = true; - 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); + rdp_simplify_recursive(loop_pts, 0, split_idx, n, eps_sq, keep); + rdp_simplify_recursive(loop_pts, split_idx, n, n, eps_sq, keep); std::vector simplified; simplified.reserve(n); @@ -164,9 +166,9 @@ static inline std::vector simplify_contour_rdp( } // C++ Seam Alignment: positions starting vertex at sharpest concave corner (or sharpest angle) -static inline std::vector align_contour_seam_sharpest_cpp(const std::vector& loop) { +static inline void align_contour_seam_sharpest_cpp_inplace(std::vector& 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 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 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 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]); +static inline std::vector align_contour_seam_sharpest_cpp(const std::vector& loop) { + std::vector aligned = loop; + align_contour_seam_sharpest_cpp_inplace(aligned); return aligned; }