up2

CMakeLists.txt

@@ -7,12 +7,13 @@ set(CMAKE_CUDA_STANDARD 14)
 set(CMAKE_CUDA_STANDARD_REQUIRED ON)
 
 # Cross-platform compiler settings
+find_package(OpenMP REQUIRED)
 if(MSVC)
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /O2 /W3 /MP")
-    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler=/MD,/W3")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /O2 /W3 /MP /openmp")
+    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler=/MD,/W3,/openmp")
 else()
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -pthread -fPIC")
-    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler -fPIC,-Wall")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -pthread -fopenmp -fPIC")
+    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -O3 --use_fast_math -Xcompiler -fPIC,-Wall,-fopenmp")
 endif()
 
 # Set broad CUDA GPU architecture support (Turing, Ampere, Ada Lovelace, Hopper/Blackwell)
@@ -83,7 +84,16 @@ execute_process(
     COMMAND "${Python3_EXECUTABLE}" -c "import pybind11; print(pybind11.get_include())"
     OUTPUT_VARIABLE _pybind11_include_dir
     OUTPUT_STRIP_TRAILING_WHITESPACE
+    ERROR_QUIET
 )
+if(NOT _pybind11_include_dir OR NOT EXISTS "${_pybind11_include_dir}")
+    execute_process(
+        COMMAND python3 -c "import pybind11; print(pybind11.get_include())"
+        OUTPUT_VARIABLE _pybind11_include_dir
+        OUTPUT_STRIP_TRAILING_WHITESPACE
+        ERROR_QUIET
+    )
+endif()
 execute_process(
     COMMAND "${Python3_EXECUTABLE}" -c "import sysconfig; print(sysconfig.get_config_var('EXT_SUFFIX') or '.so')"
     OUTPUT_VARIABLE _ext_suffix

Makefile

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

benchmark.py

@@ -0,0 +1,439 @@
+#!/usr/bin/env python3
+"""
+CUDASlicer Benchmark Suite
+==========================
+Comprehensive performance benchmark comparing:
+1. Baseline: PrusaSlicer CLI
+2. CUDASlicer Matrix Configurations:
+   - Backends: CUDA (GPU) vs CPU (Multi-threaded)
+   - Interpolation: Sub-Pixel LERP vs No LERP (Discrete Marching Squares)
+   - Grid Resolutions: 0.10 mm vs 0.05 mm Pitch
+"""
+
+from __future__ import annotations
+
+import argparse
+import os
+import shutil
+import struct
+import subprocess
+import sys
+import time
+from pathlib import Path
+from typing import Any, Dict, List, Optional, Sequence, Tuple
+
+import numpy as np
+
+# Ensure local cudaslicer python package is discoverable
+REPO_ROOT = Path(__file__).resolve().parent
+PYTHON_PKG_DIR = REPO_ROOT / "python"
+if PYTHON_PKG_DIR.exists() and str(PYTHON_PKG_DIR) not in sys.path:
+    sys.path.insert(0, str(PYTHON_PKG_DIR))
+
+try:
+    import cudaslicer
+except ImportError:
+    cudaslicer = None
+
+
+# ==============================================================================
+# Fast Binary / ASCII STL Loader
+# ==============================================================================
+def load_stl(file_path: str | Path) -> np.ndarray:
+    """
+    Load an STL file into a NumPy array of shape (N, 3, 3) representing 3D triangles.
+    Supports both Binary and ASCII STL files without external dependencies.
+    """
+    p = Path(file_path)
+    if not p.exists():
+        raise FileNotFoundError(f"STL file not found: {p}")
+
+    file_size = p.stat().st_size
+    if file_size < 84:
+        raise ValueError(f"File {p} is too small to be a valid STL.")
+
+    with open(p, "rb") as f:
+        header = f.read(80)
+        num_triangles_bytes = f.read(4)
+        num_triangles = struct.unpack("<I", num_triangles_bytes)[0]
+        expected_size = 84 + num_triangles * 50
+
+        # Check if binary STL matches expected file size
+        if expected_size == file_size:
+            f.seek(84)
+            record_dtype = np.dtype([
+                ("normal", np.float32, 3),
+                ("v0", np.float32, 3),
+                ("v1", np.float32, 3),
+                ("v2", np.float32, 3),
+                ("attr", np.uint16),
+            ])
+            data = np.fromfile(f, dtype=record_dtype, count=num_triangles)
+            return np.stack([data["v0"], data["v1"], data["v2"]], axis=1)
+
+    # Fallback to ASCII parser
+    vertices: list[list[float]] = []
+    with open(p, "r", encoding="utf-8", errors="ignore") as f:
+        for line in f:
+            line_str = line.strip()
+            if line_str.startswith("vertex"):
+                parts = line_str.split()
+                if len(parts) >= 4:
+                    vertices.append([float(parts[1]), float(parts[2]), float(parts[3])])
+
+    if len(vertices) % 3 != 0 or len(vertices) == 0:
+        raise ValueError(f"Invalid or corrupted STL file: {p}")
+
+    arr = np.array(vertices, dtype=np.float32)
+    return arr.reshape(-1, 3, 3)
+
+
+# ==============================================================================
+# PrusaSlicer Baseline Runner
+# ==============================================================================
+def find_prusaslicer_bin(custom_path: Optional[str] = None) -> Optional[str]:
+    """Find a functional PrusaSlicer executable."""
+    candidates = [
+        custom_path,
+        str(REPO_ROOT / "tests_and_tools" / "prusa_binaries" / "squashfs-root" / "usr" / "bin" / "prusa-slicer"),
+        str(REPO_ROOT / "tests_and_tools" / "prusa_binaries" / "prusa-slicer.AppImage"),
+        str(REPO_ROOT / "bin" / "prusa-slicer"),
+        "prusa-slicer",
+        "prusa-slicer-console",
+    ]
+    for cand in candidates:
+        if cand and (os.path.isfile(cand) or shutil.which(cand)):
+            return cand
+    return None
+
+
+def run_prusaslicer_benchmark(
+    stl_path: str,
+    prusa_bin: str,
+    layer_height: float = 0.20,
+    perimeters: int = 2,
+    layer_width: float = 0.45,
+    runs: int = 3,
+    expected_layers: Optional[int] = None,
+) -> Optional[Dict[str, Any]]:
+    """Execute PrusaSlicer CLI benchmark and measure execution time."""
+    tmp_gcode = REPO_ROOT / "tests_and_tools" / "gcode_outputs" / "temp_bench_prusa.gcode"
+    tmp_gcode.parent.mkdir(parents=True, exist_ok=True)
+
+    cmd = [
+        prusa_bin,
+        "--export-gcode",
+        f"--layer-height={layer_height}",
+        f"--perimeters={perimeters}",
+        "--fill-density=0%",
+        f"--extrusion-width={layer_width}",
+        f"--output={tmp_gcode}",
+        stl_path,
+    ]
+
+    timings = []
+    for r in range(runs):
+        t0 = time.perf_counter()
+        res = subprocess.run(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True)
+        t1 = time.perf_counter()
+        if res.returncode != 0:
+            print(f"\n[Warning] PrusaSlicer failed with exit code {res.returncode}:")
+            print(res.stderr)
+            return None
+        timings.append((t1 - t0) * 1000.0)
+
+    # Clean up temp file
+    if tmp_gcode.exists():
+        tmp_gcode.unlink()
+
+    return {
+        "name": "PrusaSlicer CLI (Baseline)",
+        "backend": "CPU",
+        "pitch_mm": None,
+        "subpixel_lerp": None,
+        "time_ms": float(np.mean(timings)),
+        "std_ms": float(np.std(timings)),
+        "timings": timings,
+        "layers_count": expected_layers,
+        "total_vertices": None,
+        "total_contours": None,
+    }
+
+
+# ==============================================================================
+# CUDASlicer Benchmark Runner
+# ==============================================================================
+def run_cudaslicer_combination(
+    tris: np.ndarray,
+    backend: str,
+    pitch_mm: float,
+    subpixel_lerp: bool,
+    layer_height: float = 0.20,
+    perimeters_count: int = 2,
+    layer_width: float = 0.45,
+    pad_mm: float = 1.5,
+    runs: int = 3,
+    warmup: int = 1,
+) -> Dict[str, Any]:
+    """Benchmark a specific CUDASlicer parameter combination."""
+    # Warmup
+    for _ in range(warmup):
+        _ = cudaslicer.slice_mesh_and_offset(
+            mesh_or_triangles=tris,
+            layer_height=layer_height,
+            perimeters_count=perimeters_count,
+            layer_width=layer_width,
+            pitch_mm=pitch_mm,
+            pad_mm=pad_mm,
+            backend=backend,
+            subpixel_lerp=subpixel_lerp,
+        )
+
+    timings = []
+    last_layers = None
+    for _ in range(runs):
+        t0 = time.perf_counter()
+        last_layers = cudaslicer.slice_mesh_and_offset(
+            mesh_or_triangles=tris,
+            layer_height=layer_height,
+            perimeters_count=perimeters_count,
+            layer_width=layer_width,
+            pitch_mm=pitch_mm,
+            pad_mm=pad_mm,
+            backend=backend,
+            subpixel_lerp=subpixel_lerp,
+        )
+        t1 = time.perf_counter()
+        timings.append((t1 - t0) * 1000.0)
+
+    total_vertices = 0
+    total_contours = 0
+    if last_layers:
+        for layer in last_layers:
+            for perim in getattr(layer, "perimeters", []):
+                for cnt in perim.get("contours", []):
+                    total_contours += 1
+                    total_vertices += len(cnt)
+
+    return {
+        "name": f"CUDASlicer [{backend.upper()}] {pitch_mm:.2f}mm {'LERP' if subpixel_lerp else 'No LERP'}",
+        "backend": backend.upper(),
+        "pitch_mm": pitch_mm,
+        "subpixel_lerp": subpixel_lerp,
+        "time_ms": float(np.mean(timings)),
+        "std_ms": float(np.std(timings)),
+        "timings": timings,
+        "layers_count": len(last_layers) if last_layers else 0,
+        "total_vertices": total_vertices,
+        "total_contours": total_contours,
+    }
+
+
+# ==============================================================================
+# Console Reporting & Formatting
+# ==============================================================================
+def print_header(title: str) -> None:
+    width = 98
+    print("\n" + "=" * width)
+    print(f" {title.center(width - 2)} ")
+    print("=" * width)
+
+
+def print_table(
+    results: List[Dict[str, Any]],
+    baseline_time_ms: Optional[float] = None,
+    cpu_ref_map: Optional[Dict[Tuple[float, bool], float]] = None,
+) -> None:
+    """Print a clean formatted benchmark summary table to console."""
+    # Header format
+    hdr = f"{'Configuration':<38} | {'Pitch':<6} | {'LERP':<7} | {'Time (ms)':<11} | {'ms/Layer':<9} | {'Speedup (Base)':<15} | {'GPU vs CPU':<11}"
+    sep = "-" * len(hdr)
+
+    print("\n" + sep)
+    print(hdr)
+    print(sep)
+
+    for r in results:
+        name = r["name"]
+        pitch_str = f"{r['pitch_mm']:.2f} mm" if r["pitch_mm"] is not None else "-"
+        lerp_str = "Yes" if r["subpixel_lerp"] is True else ("No" if r["subpixel_lerp"] is False else "-")
+        t_ms = r["time_ms"]
+        num_layers = r["layers_count"] or 1
+        ms_per_layer = t_ms / num_layers if num_layers > 0 else 0.0
+
+        # Baseline speedup
+        if baseline_time_ms is not None and baseline_time_ms > 0:
+            speedup_base = baseline_time_ms / t_ms
+            speedup_base_str = f"{speedup_base:6.2f}x"
+        else:
+            speedup_base_str = "-"
+
+        # GPU vs CPU speedup for same pitch & lerp
+        gpu_vs_cpu_str = "-"
+        if r["backend"] == "CUDA" and cpu_ref_map:
+            key = (r["pitch_mm"], r["subpixel_lerp"])
+            if key in cpu_ref_map and cpu_ref_map[key] > 0:
+                gpu_speedup = cpu_ref_map[key] / t_ms
+                gpu_vs_cpu_str = f"{gpu_speedup:6.2f}x"
+
+        print(
+            f"{name:<38} | {pitch_str:<6} | {lerp_str:<7} | {t_ms:9.2f} ± {r['std_ms']:3.1f} | {ms_per_layer:8.2f} | {speedup_base_str:^15} | {gpu_vs_cpu_str:^11}"
+        )
+
+    print(sep + "\n")
+
+
+# ==============================================================================
+# Main Benchmark CLI
+# ==============================================================================
+def main() -> None:
+    parser = argparse.ArgumentParser(
+        description="CUDASlicer Performance Benchmark Suite (PrusaSlicer Baseline vs CPU/GPU, LERP/No LERP, 0.05mm/0.10mm pitch)"
+    )
+    default_stl = REPO_ROOT / "benchy.stl"
+    parser.add_argument("--stl", type=str, default=str(default_stl), help="Path to STL mesh model (default: benchy.stl)")
+    parser.add_argument("--layer-height", type=float, default=0.20, help="Slicing layer height in mm (default: 0.20)")
+    parser.add_argument("--perimeters", type=int, default=2, help="Number of perimeter passes (default: 2)")
+    parser.add_argument("--layer-width", type=float, default=0.45, help="Extrusion line width in mm (default: 0.45)")
+    parser.add_argument("--pitches", nargs="+", type=float, default=[0.10, 0.05], help="Grid pixel pitches in mm (default: 0.10 0.05)")
+    parser.add_argument("--runs", type=int, default=3, help="Number of benchmark iterations to average (default: 3)")
+    parser.add_argument("--warmup", type=int, default=1, help="Number of warmup runs (default: 1)")
+    parser.add_argument("--no-prusa", action="store_true", help="Skip PrusaSlicer baseline benchmark")
+    parser.add_argument("--prusa-bin", type=str, default=None, help="Path to PrusaSlicer executable")
+    parser.add_argument("--gpu-only", action="store_true", help="Run only GPU (CUDA) benchmarks")
+    parser.add_argument("--cpu-only", action="store_true", help="Run only CPU benchmarks")
+    parser.add_argument("--lerp-only", action="store_true", help="Run only Sub-pixel LERP passes")
+    parser.add_argument("--no-lerp-only", action="store_true", help="Run only No-LERP (discrete) passes")
+
+    args = parser.parse_args()
+
+    print_header("CUDASLICER PERFORMANCE BENCHMARK SUITE")
+
+    # 1. System & Hardware Details
+    stl_path = Path(args.stl).resolve()
+    if not stl_path.exists():
+        print(f"[Error] Specified STL file does not exist: {stl_path}")
+        sys.exit(1)
+
+    print(f"Loading Model: {stl_path.name}")
+    tris = load_stl(stl_path)
+    num_tris = len(tris)
+    bounds_min = tris.min(axis=(0, 1))
+    bounds_max = tris.max(axis=(0, 1))
+    dimensions = bounds_max - bounds_min
+
+    print(f"  • Triangles:      {num_tris:,}")
+    print(f"  • Dimensions:     {dimensions[0]:.2f} x {dimensions[1]:.2f} x {dimensions[2]:.2f} mm")
+    print(f"  • Layer Height:   {args.layer_height:.2f} mm")
+    print(f"  • Perimeters:     {args.perimeters} passes @ {args.layer_width:.2f} mm line width")
+    print(f"  • Iterations:     {args.runs} runs (warmup: {args.warmup})")
+
+    # Check cudaslicer availability
+    if not cudaslicer or not cudaslicer.is_available():
+        print("[Error] cudaslicer native backend is not loaded or available.")
+        sys.exit(1)
+
+    has_gpu = cudaslicer.has_cuda()
+    print(f"  • CUDA GPU:       {'Available' if has_gpu else 'Not Available'}")
+    if not has_gpu and not args.cpu_only:
+        print("    (Notice: No CUDA GPU detected, forcing CPU execution)")
+
+    results: List[Dict[str, Any]] = []
+
+    expected_layers = int(np.ceil(dimensions[2] / args.layer_height))
+
+    # 2. Benchmark Baseline: PrusaSlicer
+    prusa_res = None
+    if not args.no_prusa:
+        prusa_bin = find_prusaslicer_bin(args.prusa_bin)
+        if prusa_bin:
+            print(f"\n[1/3] Benchmarking Baseline: PrusaSlicer CLI ({Path(prusa_bin).name})...")
+            prusa_res = run_prusaslicer_benchmark(
+                stl_path=str(stl_path),
+                prusa_bin=prusa_bin,
+                layer_height=args.layer_height,
+                perimeters=args.perimeters,
+                layer_width=args.layer_width,
+                runs=args.runs,
+                expected_layers=expected_layers,
+            )
+            if prusa_res:
+                print(f"  -> PrusaSlicer Baseline: {prusa_res['time_ms']:.2f} ± {prusa_res['std_ms']:.2f} ms")
+                results.append(prusa_res)
+        else:
+            print("\n[Warning] PrusaSlicer binary not found. Skipping baseline comparison.")
+
+    # 3. Build combination matrix
+    backends: list[str] = []
+    if args.cpu_only:
+        backends = ["cpu"]
+    elif args.gpu_only:
+        backends = ["cuda"] if has_gpu else []
+    else:
+        backends = (["cuda"] if has_gpu else []) + ["cpu"]
+
+    lerp_modes: list[bool] = []
+    if args.lerp_only:
+        lerp_modes = [True]
+    elif args.no_lerp_only:
+        lerp_modes = [False]
+    else:
+        lerp_modes = [True, False]
+
+    pitches = args.pitches
+
+    combo_idx = 1
+    cpu_ref_map: Dict[Tuple[float, bool], float] = {}
+
+    for backend in backends:
+        for pitch in pitches:
+            if backend.lower() == "cpu" and pitch < 0.09 and not args.cpu_only:
+                # Skip 0.05mm on CPU by default as it operates on 200M+ dense grid on CPU
+                continue
+
+            for lerp in lerp_modes:
+                lerp_tag = "LERP" if lerp else "No LERP"
+                print(
+                    f"  [{combo_idx}] Running {backend.upper()} | Pitch: {pitch:.2f}mm | {lerp_tag} ...",
+                    end="",
+                    flush=True,
+                )
+                res = run_cudaslicer_combination(
+                    tris=tris,
+                    backend=backend,
+                    pitch_mm=pitch,
+                    subpixel_lerp=lerp,
+                    layer_height=args.layer_height,
+                    perimeters_count=args.perimeters,
+                    layer_width=args.layer_width,
+                    runs=args.runs,
+                    warmup=args.warmup,
+                )
+                print(f" {res['time_ms']:.2f} ms ({res['time_ms']/res['layers_count']:.2f} ms/layer)")
+                results.append(res)
+                if backend.lower() == "cpu":
+                    cpu_ref_map[(pitch, lerp)] = res["time_ms"]
+                combo_idx += 1
+
+    # 4. Print Summary Report
+    print_header("BENCHMARK RESULTS & SPEEDUP COMPARISON")
+    baseline_time = prusa_res["time_ms"] if prusa_res else None
+    print_table(results, baseline_time_ms=baseline_time, cpu_ref_map=cpu_ref_map)
+
+    # Key Highlights
+    gpu_results = [r for r in results if r["backend"] == "CUDA"]
+    if gpu_results and baseline_time:
+        fastest_gpu = min(gpu_results, key=lambda x: x["time_ms"])
+        speedup = baseline_time / fastest_gpu["time_ms"]
+        print("Performance Highlights:")
+        print(f"  ★ Fastest GPU Configuration: {fastest_gpu['name']} in {fastest_gpu['time_ms']:.2f} ms")
+        print(f"  ★ GPU Speedup vs PrusaSlicer Baseline: {speedup:.2f}x faster")
+        if cpu_ref_map:
+            key = (fastest_gpu["pitch_mm"], fastest_gpu["subpixel_lerp"])
+            if key in cpu_ref_map:
+                gpu_vs_cpu = cpu_ref_map[key] / fastest_gpu["time_ms"]
+                print(f"  ★ GPU Acceleration vs Multi-threaded CPU: {gpu_vs_cpu:.2f}x faster")
+
+
+if __name__ == "__main__":
+    main()

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 y = blockIdx.y * blockDim.y + threadIdx.y;
-    int l = blockIdx.z;
+    int tx = threadIdx.x;
+    int ty = threadIdx.y;
+    int tid = ty * blockDim.x + tx;
 
-    if (x >= cols || y >= rows || l >= num_layers) return;
+    int x = blockIdx.x * blockDim.x + tx;
+    int y = blockIdx.y * blockDim.y + ty;
+    int l = blockIdx.z;
 
     size_t layer_pixel_offset = (size_t)l * rows * cols;
     size_t idx = layer_pixel_offset + (size_t)y * cols + x;
@@ -715,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) {
-        GpuSegment2D seg = segs[s];
+    __shared__ GpuSegment2D s_active_segs[256];
+    __shared__ int s_active_count;
+
+    for (int chunk_start = 0; chunk_start < seg_count; chunk_start += 256) {
+        if (tid == 0) {
+            s_active_count = 0;
+        }
+        __syncthreads();
+
+        int s_idx = chunk_start + tid;
+        if (s_idx < seg_count) {
+            GpuSegment2D seg = segs[s_idx];
+            float sx_min = fminf(seg.x0, seg.x1);
+            float sx_max = fmaxf(seg.x0, seg.x1);
+            float sy_min = fminf(seg.y0, seg.y1);
+            float sy_max = fmaxf(seg.y0, seg.y1);
+
+            if (sx_max >= bx0 && sx_min <= bx1 && sy_max >= by0 && sy_min <= by1) {
+                int pos = atomicAdd(&s_active_count, 1);
+                if (pos < 256) {
+                    s_active_segs[pos] = seg;
+                }
+            }
+        }
+        __syncthreads();
+
+        int num_active = min(s_active_count, 256);
+        if (num_active > 0 && x < cols && y < rows) {
+            for (int k = 0; k < num_active; ++k) {
+                GpuSegment2D seg = s_active_segs[k];
                 float vx = seg.x1 - seg.x0;
                 float vy = seg.y1 - seg.y0;
                 float wx = px - seg.x0;
@@ -742,10 +809,20 @@ __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(
@@ -805,26 +882,26 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_mesh(
     if (max_segs_per_layer < 256) max_segs_per_layer = 256;
     size_t total_segs = (size_t)num_layers * max_segs_per_layer;
 
-    GpuTriangle3D* d_triangles = nullptr;
-    unsigned int* d_parity = nullptr;
-    uint8_t* d_masks = nullptr;
-    GpuSegment2D* d_layer_segments = nullptr;
-    int* d_layer_seg_counts = nullptr;
+    size_t tri_bytes_aligned = (tri_bytes + 255) & ~255ULL;
+    size_t parity_bytes_aligned = (total_words * sizeof(unsigned int) + 255) & ~255ULL;
+    size_t masks_bytes_aligned = (total_volume * sizeof(uint8_t) + 255) & ~255ULL;
+    size_t segs_bytes_aligned = (total_segs * sizeof(GpuSegment2D) + 255) & ~255ULL;
+    size_t counts_bytes_aligned = (num_layers * sizeof(int) + 255) & ~255ULL;
 
-    if (cudaMalloc(&d_triangles, tri_bytes) != cudaSuccess ||
-        cudaMalloc(&d_parity, total_words * sizeof(unsigned int)) != cudaSuccess ||
-        cudaMalloc(&d_masks, total_volume * sizeof(uint8_t)) != cudaSuccess ||
-        cudaMalloc(&d_layer_segments, total_segs * sizeof(GpuSegment2D)) != cudaSuccess ||
-        cudaMalloc(&d_layer_seg_counts, num_layers * sizeof(int)) != cudaSuccess) {
-        if (d_triangles) cudaFree(d_triangles);
-        if (d_parity) cudaFree(d_parity);
-        if (d_masks) cudaFree(d_masks);
-        if (d_layer_segments) cudaFree(d_layer_segments);
-        if (d_layer_seg_counts) cudaFree(d_layer_seg_counts);
+    size_t workspace_total = tri_bytes_aligned + parity_bytes_aligned + masks_bytes_aligned + segs_bytes_aligned + counts_bytes_aligned;
+
+    uint8_t* d_workspace = nullptr;
+    if (cudaMalloc(&d_workspace, workspace_total) != cudaSuccess) {
         gpu_sdf_destroy(stack);
         return nullptr;
     }
 
+    GpuTriangle3D* d_triangles = (GpuTriangle3D*)(d_workspace);
+    unsigned int* d_parity = (unsigned int*)(d_workspace + tri_bytes_aligned);
+    uint8_t* d_masks = (uint8_t*)(d_workspace + tri_bytes_aligned + parity_bytes_aligned);
+    GpuSegment2D* d_layer_segments = (GpuSegment2D*)(d_workspace + tri_bytes_aligned + parity_bytes_aligned + masks_bytes_aligned);
+    int* d_layer_seg_counts = (int*)(d_workspace + tri_bytes_aligned + parity_bytes_aligned + masks_bytes_aligned + segs_bytes_aligned);
+
     cudaMemsetAsync(d_parity, 0, total_words * sizeof(unsigned int), stack->stream);
     cudaMemsetAsync(d_layer_seg_counts, 0, num_layers * sizeof(int), stack->stream);
     cudaMemcpyAsync(d_triangles, triangles, tri_bytes, cudaMemcpyHostToDevice, stack->stream);
@@ -851,12 +928,7 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_from_mesh(
     );
 
     cudaStreamSynchronize(stack->stream);
-
-    cudaFree(d_triangles);
-    cudaFree(d_parity);
-    cudaFree(d_masks);
-    cudaFree(d_layer_segments);
-    cudaFree(d_layer_seg_counts);
+    cudaFree(d_workspace);
 
     return stack;
 }
@@ -1102,6 +1174,284 @@ GPU_SDF_API GpuSdfStackHandle gpu_sdf_shell(
     return out_stack;
 }
 
+struct GpuDirectedSegment {
+    int src_edge;
+    int dst_edge;
+    GpuPoint2D pt;
+};
+
+__device__ inline GpuPoint2D interpolate_gpu_edge(
+    int edge_idx, int x, int y,
+    float c0, float c1, float c2, float c3,
+    float iso_val, float min_x, float min_y, float pitch_mm, int subpixel_lerp
+) {
+    if (!subpixel_lerp) {
+        if (edge_idx == 0) return GpuPoint2D{ min_x + (x + 0.5f) * pitch_mm, min_y + (float)y * pitch_mm };
+        if (edge_idx == 1) return GpuPoint2D{ min_x + (float)(x + 1) * pitch_mm, min_y + (y + 0.5f) * pitch_mm };
+        if (edge_idx == 2) return GpuPoint2D{ min_x + (x + 0.5f) * pitch_mm, min_y + (float)(y + 1) * pitch_mm };
+        /* edge_idx == 3 */return GpuPoint2D{ min_x + (float)x * pitch_mm, min_y + (y + 0.5f) * pitch_mm };
+    }
+
+    float t = 0.5f;
+    if (edge_idx == 0) { // Bottom (c0 -> c1)
+        float denom = c1 - c0;
+        if (fabsf(denom) > 1e-7f) {
+            t = (iso_val - c0) / denom;
+            t = fminf(fmaxf(t, 0.0f), 1.0f);
+        }
+        return GpuPoint2D{ min_x + ((float)x + t) * pitch_mm, min_y + (float)y * pitch_mm };
+    } else if (edge_idx == 1) { // Right (c1 -> c2)
+        float denom = c2 - c1;
+        if (fabsf(denom) > 1e-7f) {
+            t = (iso_val - c1) / denom;
+            t = fminf(fmaxf(t, 0.0f), 1.0f);
+        }
+        return GpuPoint2D{ min_x + (float)(x + 1) * pitch_mm, min_y + ((float)y + t) * pitch_mm };
+    } else if (edge_idx == 2) { // Top (c3 -> c2)
+        float denom = c2 - c3;
+        if (fabsf(denom) > 1e-7f) {
+            t = (iso_val - c3) / denom;
+            t = fminf(fmaxf(t, 0.0f), 1.0f);
+        }
+        return GpuPoint2D{ min_x + ((float)x + t) * pitch_mm, min_y + (float)(y + 1) * pitch_mm };
+    } else { // Left (c0 -> c3)
+        float denom = c3 - c0;
+        if (fabsf(denom) > 1e-7f) {
+            t = (iso_val - c0) / denom;
+            t = fminf(fmaxf(t, 0.0f), 1.0f);
+        }
+        return GpuPoint2D{ min_x + (float)x * pitch_mm, min_y + ((float)y + t) * pitch_mm };
+    }
+}
+
+__global__ void gpu_ms_extract_segments_kernel(
+    const float* __restrict__ d_sdf,
+    int cols, int rows, int num_layers,
+    float bounds_min_x, float bounds_min_y, float pitch_mm,
+    const float* __restrict__ d_offsets,
+    int num_offsets,
+    int subpixel_lerp,
+    int max_segs_per_pass,
+    GpuDirectedSegment* __restrict__ d_out_segments,
+    int* __restrict__ d_out_segment_counts
+) {
+    int x = blockIdx.x * blockDim.x + threadIdx.x;
+    int y = blockIdx.y * blockDim.y + threadIdx.y;
+    int pass_id = blockIdx.z;
+
+    int total_passes = num_layers * num_offsets;
+    if (pass_id >= total_passes) return;
+
+    int l = pass_id / num_offsets;
+    int o = pass_id % num_offsets;
+    float iso_val = d_offsets[o];
+
+    int stride_w = cols + 1;
+    size_t layer_offset = (size_t)l * rows * cols;
+    const float* layer_sdf = d_sdf + layer_offset;
+
+    if (x >= cols - 1 || y >= rows - 1) return;
+
+    const float* row0 = layer_sdf + (size_t)y * cols;
+    const float* row1 = layer_sdf + (size_t)(y + 1) * cols;
+
+    float c0 = row0[x];
+    float c1 = row0[x + 1];
+    float c2 = row1[x + 1];
+    float c3 = row1[x];
+
+    int b0 = (c0 >= iso_val) ? 1 : 0;
+    int b1 = (c1 >= iso_val) ? 1 : 0;
+    int b2 = (c2 >= iso_val) ? 1 : 0;
+    int b3 = (c3 >= iso_val) ? 1 : 0;
+
+    int case_idx = b0 | (b1 << 1) | (b2 << 2) | (b3 << 3);
+    if (case_idx == 0 || case_idx == 15) return;
+
+    int e_bottom = 2 * (y * stride_w + x) + 0;
+    int e_top    = 2 * ((y + 1) * stride_w + x) + 0;
+    int e_left   = 2 * (y * stride_w + x) + 1;
+    int e_right  = 2 * (y * stride_w + (x + 1)) + 1;
+
+    auto emit_segment = [&](int local_src, int global_src, int global_dst) {
+        GpuPoint2D pt = interpolate_gpu_edge(local_src, x, y, c0, c1, c2, c3, iso_val, bounds_min_x, bounds_min_y, pitch_mm, subpixel_lerp);
+        int pos = atomicAdd(&d_out_segment_counts[pass_id], 1);
+        if (pos < max_segs_per_pass) {
+            d_out_segments[(size_t)pass_id * max_segs_per_pass + pos] = GpuDirectedSegment{ global_src, global_dst, pt };
+        }
+    };
+
+    switch (case_idx) {
+        case 1:  emit_segment(3, e_left, e_bottom); break;
+        case 2:  emit_segment(0, e_bottom, e_right); break;
+        case 3:  emit_segment(3, e_left, e_right); break;
+        case 4:  emit_segment(1, e_right, e_top); break;
+        case 5: {
+            float center = (c0 + c1 + c2 + c3) * 0.25f;
+            if (center >= iso_val) {
+                emit_segment(3, e_left, e_top);
+                emit_segment(0, e_bottom, e_right);
+            } else {
+                emit_segment(3, e_left, e_bottom);
+                emit_segment(1, e_right, e_top);
+            }
+            break;
+        }
+        case 6:  emit_segment(0, e_bottom, e_top); break;
+        case 7:  emit_segment(3, e_left, e_top); break;
+        case 8:  emit_segment(2, e_top, e_left); break;
+        case 9:  emit_segment(2, e_top, e_bottom); break;
+        case 10: {
+            float center = (c0 + c1 + c2 + c3) * 0.25f;
+            if (center >= iso_val) {
+                emit_segment(0, e_bottom, e_left);
+                emit_segment(2, e_top, e_right);
+            } else {
+                emit_segment(2, e_top, e_left);
+                emit_segment(0, e_bottom, e_right);
+            }
+            break;
+        }
+        case 11: emit_segment(2, e_top, e_right); break;
+        case 12: emit_segment(1, e_right, e_left); break;
+        case 13: emit_segment(1, e_right, e_bottom); break;
+        case 14: emit_segment(0, e_bottom, e_left); break;
+    }
+}
+
+__global__ void gpu_prefix_scan_counts_kernel(
+    const int* __restrict__ d_counts,
+    int* __restrict__ d_offsets,
+    int total_passes
+) {
+    if (threadIdx.x == 0 && blockIdx.x == 0) {
+        int sum = 0;
+        for (int i = 0; i < total_passes; ++i) {
+            d_offsets[i] = sum;
+            sum += d_counts[i];
+        }
+        d_offsets[total_passes] = sum;
+    }
+}
+
+__global__ void gpu_compact_segments_kernel(
+    const GpuDirectedSegment* __restrict__ d_src_segments,
+    const int* __restrict__ d_counts,
+    const int* __restrict__ d_pass_offsets,
+    int max_segs_per_pass,
+    int total_passes,
+    GpuDirectedSegment* __restrict__ d_dst_segments
+) {
+    int p = blockIdx.x;
+    int tid = threadIdx.x;
+    if (p >= total_passes) return;
+
+    int count = d_counts[p];
+    int dst_offset = d_pass_offsets[p];
+    const GpuDirectedSegment* src = d_src_segments + (size_t)p * max_segs_per_pass;
+    GpuDirectedSegment* dst = d_dst_segments + dst_offset;
+
+    for (int i = tid; i < count; i += blockDim.x) {
+        dst[i] = src[i];
+    }
+}
+
+struct CompactEdgeMap {
+    int cap;
+    int mask;
+    int* keys;
+    int* dsts;
+    GpuPoint2D* pts;
+    uint8_t* visited;
+
+    CompactEdgeMap() : cap(0), mask(0), keys(nullptr), dsts(nullptr), pts(nullptr), visited(nullptr) {}
+
+    void ensure_capacity(int required) {
+        int target = 128;
+        while (target < required * 2) target <<= 1;
+        if (target <= cap) return;
+
+        if (keys) { free(keys); free(dsts); free(pts); free(visited); }
+        cap = target;
+        mask = cap - 1;
+        keys = (int*)malloc(cap * sizeof(int));
+        dsts = (int*)malloc(cap * sizeof(int));
+        pts = (GpuPoint2D*)malloc(cap * sizeof(GpuPoint2D));
+        visited = (uint8_t*)malloc(cap * sizeof(uint8_t));
+        memset(keys, -1, cap * sizeof(int));
+        memset(visited, 0, cap * sizeof(uint8_t));
+    }
+    void free_mem() {
+        if (keys) { free(keys); free(dsts); free(pts); free(visited); }
+        keys = nullptr; cap = 0;
+    }
+    inline void insert(int src, int dst, GpuPoint2D pt) {
+        int idx = ((unsigned int)src * 2654435761u) & mask;
+        while (keys[idx] != -1 && keys[idx] != src) {
+            idx = (idx + 1) & mask;
+        }
+        keys[idx] = src;
+        dsts[idx] = dst;
+        pts[idx] = pt;
+        visited[idx] = 0;
+    }
+    inline int find_slot(int src) const {
+        int idx = ((unsigned int)src * 2654435761u) & mask;
+        while (keys[idx] != -1) {
+            if (keys[idx] == src) return idx;
+            idx = (idx + 1) & mask;
+        }
+        return -1;
+    }
+    inline void clear_slot(int idx) {
+        keys[idx] = -1;
+        visited[idx] = 0;
+    }
+};
+
+struct GpuVecWorkspace {
+    GpuDirectedSegment* d_segments = nullptr;
+    GpuDirectedSegment* d_compact_segments = nullptr;
+    int* d_counts = nullptr;
+    int* d_offsets = nullptr;
+    float* d_iso_offsets = nullptr;
+    int* h_offsets = nullptr;
+    GpuDirectedSegment* h_compact_segments = nullptr;
+    size_t max_total_segs = 0;
+    int total_passes_cap = 0;
+    int num_offsets_cap = 0;
+    std::mutex mtx;
+
+    void ensure(int total_passes, int max_segs_per_pass, int num_offsets) {
+        std::lock_guard<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;
-    size_t sdf_bytes = total_voxels * sizeof(float);
+    g_vec_workspace.ensure(total_passes, max_segs_per_pass, (int)num_offsets);
 
-    float* h_pinned_sdf = nullptr;
-    CUDA_CHECK_INT(cudaMallocHost(&h_pinned_sdf, sdf_bytes));
-    CUDA_CHECK_INT(cudaMemcpyAsync(h_pinned_sdf, stack->d_sdf, sdf_bytes,
-                                   cudaMemcpyDeviceToHost, stack->stream));
-    CUDA_CHECK_INT(cudaStreamSynchronize(stack->stream));
+    cudaMemcpyAsync(g_vec_workspace.d_iso_offsets, resolved_offsets.data(), num_offsets * sizeof(float), cudaMemcpyHostToDevice, stack->stream);
+    cudaMemsetAsync(g_vec_workspace.d_counts, 0, total_passes * sizeof(int), stack->stream);
+
+    auto t_k_0 = std::chrono::high_resolution_clock::now();
+
+    dim3 block(16, 16);
+    dim3 grid((cols + 15) / 16, (rows + 15) / 16, total_passes);
+    gpu_ms_extract_segments_kernel<<<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());
-    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) {
+    for (size_t l = 0; l < num_layers; ++l) {
         res_array[l].z_height = (float)z_heights_um[l] * 0.001f;
         res_array[l].num_perimeters = P;
         res_array[l].perimeters = (GpuOffsetPerimeter*)malloc(P * sizeof(GpuOffsetPerimeter));
-
-            const float* layer_sdf_ptr = h_pinned_sdf + l * ((size_t)rows * cols);
-
+        memset(res_array[l].perimeters, 0, P * sizeof(GpuOffsetPerimeter));
         for (int o = 0; o < P; ++o) {
-                float iso_thresh = 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].perimeters[o].offset_distance_mm = resolved_offsets[o];
         }
-            }
-
-            // 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.offset_distance_mm = resolved_offsets[P];
         res_array[l].infill.width = cols;
         res_array[l].infill.height = rows;
-            res_array[l].infill.num_contours = infill_contours.size();
-            if (!infill_contours.empty()) {
-                res_array[l].infill.contours = (GpuContour*)malloc(infill_contours.size() * sizeof(GpuContour));
-                memcpy(res_array[l].infill.contours, infill_contours.data(), infill_contours.size() * sizeof(GpuContour));
-            } else {
-                res_array[l].infill.contours = nullptr;
     }
 
-            if (config->export_infill_mask) {
-                uint8_t* infill_mask = (uint8_t*)malloc((size_t)cols * rows);
-                for (int i = 0; i < cols * rows; ++i) {
-                    infill_mask[i] = (layer_sdf_ptr[i] >= infill_iso_thresh) ? 1 : 0;
-                }
-                res_array[l].infill.mask = infill_mask;
-            } else {
-                res_array[l].infill.mask = nullptr;
-            }
-        }
-    };
+    #pragma omp parallel
+    {
+        CompactEdgeMap map;
+        std::vector<GpuPoint2D> loop_pts;
+        loop_pts.reserve(2000);
+        std::vector<int> active_slots;
+        active_slots.reserve(4000);
+        std::vector<GpuContour> pass_contours;
+        pass_contours.reserve(32);
 
-    size_t chunk = (num_layers + num_threads - 1) / num_threads;
-    for (unsigned int t = 0; t < num_threads; ++t) {
-        size_t s = t * chunk;
-        size_t e = std::min(s + chunk, num_layers);
-        if (s < e) workers.emplace_back(vec_worker, s, e);
+        #pragma omp for schedule(dynamic, 8)
+        for (int p = 0; p < total_passes; ++p) {
+            size_t l = (size_t)(p / (int)num_offsets);
+            int o = p % (int)num_offsets;
+
+            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 (x >= 0 && x < width) d0 = sdf[y * width + x];
-            if (x + 1 >= 0 && x + 1 < width) d1 = sdf[y * width + (x + 1)];
+        if ((unsigned int)y < (unsigned int)height) {
+            size_t row_idx = (size_t)y * width;
+            if ((unsigned int)x < (unsigned int)width) d0 = sdf[row_idx + x];
+            if ((unsigned int)(x + 1) < (unsigned int)width) d1 = sdf[row_idx + (x + 1)];
         }
         float t = 0.5f;
         float denom = d1 - d0;
@@ -62,9 +63,9 @@ inline GpuPoint2D get_sdf_edge_point(
         // Vertical edge: from (x, y) to (x, y + 1)
         float d0 = -1e5f;
         float d3 = -1e5f;
-        if (x >= 0 && x < width) {
-            if (y >= 0 && y < height) d0 = sdf[y * width + x];
-            if (y + 1 >= 0 && y + 1 < height) d3 = sdf[(y + 1) * width + x];
+        if ((unsigned int)x < (unsigned int)width) {
+            if ((unsigned int)y < (unsigned int)height) d0 = sdf[(size_t)y * width + x];
+            if ((unsigned int)(y + 1) < (unsigned int)height) d3 = sdf[(size_t)(y + 1) * width + x];
         }
         float t = 0.5f;
         float denom = d3 - d0;
@@ -162,6 +163,70 @@ static inline std::vector<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);
             }