arachne

include/gpu_offset.h

@@ -83,12 +83,24 @@ typedef struct {
 } GpuInfillRegion;
 
 /**
- * @brief All offset perimeters and infill mask for a single sliced layer
+ * @brief Variable-width toolpath bead (e.g. Medial Axis Transform / Arachne centerline or perimeter pass)
+ */
+typedef struct {
+    GpuPoint2D* points;       // Polyline vertices (X, Y in mm)
+    float* widths;            // Local extrusion line width at each vertex in mm
+    size_t num_points;        // Number of vertices
+    int is_closed;            // 0 = open bead (centerline / gap fill), 1 = closed loop
+} GpuBead;
+
+/**
+ * @brief All offset perimeters, Arachne variable-width beads, and infill mask for a single sliced layer
  */
 typedef struct {
     GpuOffsetPerimeter* perimeters; // Array of perimeters (length = num_perimeters)
     size_t num_perimeters;         // perimeters_count
     GpuInfillRegion infill;        // Infill boundary contours & 2D mask
+    GpuBead* arachne_beads;        // Array of variable-width Arachne centerline / gap-fill beads
+    size_t num_arachne_beads;      // Count of variable-width beads
     float z_height;
 } GpuLayerResult;
 
@@ -108,6 +120,10 @@ typedef struct {
     float simplify_tolerance_mm; // RDP simplification tolerance in mm (default: 0.02 mm, 0 to disable)
     int subpixel_lerp;         // 1 = Continuous SDF sub-pixel LERP vectorization (default: 1), 0 = discrete bitpack
     int backend;               // GpuOffsetBackend: 0 = AUTO, 1 = CUDA, 2 = CPU (default: 0)
+    int enable_arachne;        // 1 = Enable Medial Axis Transform (Arachne) variable line width for collapsing walls (default: 0)
+    float min_feature_size_mm; // Minimum printable feature thickness (e.g. 0.20 mm, default: 0.20 mm)
+    float min_bead_width_mm;   // Minimum extrusion width for variable bead (e.g. 0.25 mm, default: 0.25 mm)
+    float max_bead_width_mm;   // Maximum extrusion width for variable bead (e.g. 0.65 mm, default: 0.65 mm)
 } GpuOffsetConfig;
 
 /**

python/cudaslicer/__init__.py

@@ -177,19 +177,20 @@ def offset_polygons(
     return_infill_mask: bool = True,
     offset_distances: Optional[Sequence[float]] = None,
     backend: Optional[str] = "auto",
-    output_3d: bool = False
+    output_3d: bool = False,
+    enable_arachne: bool = False,
+    min_feature_size_mm: float = 0.20,
+    min_bead_width_mm: float = 0.25,
+    max_bead_width_mm: float = 0.65,
+    simplify_tolerance_mm: float = 0.02
 ) -> Union[LayerResult, List[LayerResult]]:
     """
-    Offset 2D polygon loops across one or more layers (CUDA or CPU fallback).
+    Offset 2D polygon loops across one or more layers (CUDA or CPU fallback) with optional Arachne MAT.
 
     Parameters
     ----------
     layers : list of np.ndarray OR list of dict
-        Input layer data. Can be:
-        1. Single layer with a list of (K, 2) NumPy polygon contours:
-           `[poly1, poly2, ...]`
-        2. Multi-layer list of layer dicts:
-           `[{"z": 0.2, "contours": [poly1, ...]}, {"z": 0.4, "contours": [...]}]`
+        Input layer data.
     perimeters_count : int, default=2
         Number of perimeter passes (Perimeter 1: -0.5*width, Perimeter 2: -1.5*width, etc.).
     layer_width : float, default=0.45
@@ -204,12 +205,16 @@ def offset_polygons(
         Custom offset distances.
     backend : str, optional, default='auto'
         Compute backend: 'auto', 'cuda', or 'cpu'.
-
-    Returns
-    -------
-    LayerResult or list of LayerResult
-        If single layer: returns single LayerResult object.
-        If multi-layer: returns list of LayerResult objects.
+    enable_arachne : bool, default=False
+        Enable Medial Axis Transform (MAT) variable line width for collapsing walls / thin features.
+    min_feature_size_mm : float, default=0.20
+        Minimum printable wall thickness in millimeters.
+    min_bead_width_mm : float, default=0.25
+        Minimum variable extrusion bead width in millimeters.
+    max_bead_width_mm : float, default=0.65
+        Maximum variable extrusion bead width in millimeters.
+    simplify_tolerance_mm : float, default=0.02
+        RDP polyline simplification tolerance in millimeters.
     """
     if not _HAS_BACKEND or _backend is None:
         raise RuntimeError("CUDASlicer backend extension (_cudaslicer_backend) is not available.")
@@ -235,7 +240,7 @@ def offset_polygons(
                 all_pts.append(c)
 
     if not all_pts:
-        empty_res = LayerResult({"z": 0.0, "perimeters": [[] for _ in range(perimeters_count)], "infill_contours": [], "infill_mask": None})
+        empty_res = LayerResult({"z": 0.0, "perimeters": [[] for _ in range(perimeters_count)], "infill_contours": [], "infill_mask": None, "arachne_beads": []})
         return empty_res if is_single_layer_format else []
 
     stacked = np.concatenate(all_pts, axis=0)
@@ -251,7 +256,12 @@ def offset_polygons(
         float(pitch_mm), float(pad_mm),
         bool(return_infill_mask), offsets_float,
         backend_str,
-        bool(output_3d)
+        bool(output_3d),
+        bool(enable_arachne),
+        float(min_feature_size_mm),
+        float(min_bead_width_mm),
+        float(max_bead_width_mm),
+        float(simplify_tolerance_mm)
     )
 
     layer_results = [LayerResult(r) for r in raw_results]
@@ -273,10 +283,14 @@ def slice_mesh_and_offset(
     offset_distances: Optional[Sequence[float]] = None,
     backend: Optional[str] = "auto",
     simplify_tolerance_mm: float = 0.0,
-    subpixel_lerp: bool = True
+    subpixel_lerp: bool = True,
+    enable_arachne: bool = False,
+    min_feature_size_mm: float = 0.20,
+    min_bead_width_mm: float = 0.25,
+    max_bead_width_mm: float = 0.65
 ) -> List[LayerResult]:
     """
-    3D mesh slicing, rasterization, distance field, and perimeter offsetting (CUDA or CPU).
+    3D mesh slicing, rasterization, distance field, perimeter offsetting, and optional Arachne MAT.
 
     Parameters
     ----------
@@ -304,15 +318,14 @@ def slice_mesh_and_offset(
         Ramer-Douglas-Peucker (RDP) polygon simplification tolerance in mm (0.0 to disable).
     subpixel_lerp : bool, default=True
         Continuous distance-field sub-pixel linear interpolation for smooth arbitrary angles.
-
-    Returns
-    -------
-    list of LayerResult
-        List of LayerResult dictionaries containing:
-        - 'z': Layer Z height
-        - 'perimeters': List of perimeters (each containing list of (K, 2) contour arrays)
-        - 'infill_contours': Infill boundary contour loops
-        - 'infill_mask': 2D binary NumPy array (H, W) of the infill region
+    enable_arachne : bool, default=False
+        Enable Medial Axis Transform (MAT) variable line width for collapsing walls / thin features.
+    min_feature_size_mm : float, default=0.20
+        Minimum printable wall thickness in millimeters.
+    min_bead_width_mm : float, default=0.25
+        Minimum variable extrusion bead width in millimeters.
+    max_bead_width_mm : float, default=0.65
+        Maximum variable extrusion bead width in millimeters.
     """
     if not _HAS_BACKEND or _backend is None:
         raise RuntimeError("CUDASlicer backend extension (_cudaslicer_backend) is not available.")
@@ -348,11 +361,18 @@ def slice_mesh_and_offset(
         bool(return_infill_mask), offsets_float,
         backend_str,
         float(simplify_tolerance_mm),
-        bool(subpixel_lerp)
+        bool(subpixel_lerp),
+        False,
+        bool(enable_arachne),
+        float(min_feature_size_mm),
+        float(min_bead_width_mm),
+        float(max_bead_width_mm)
     )
 
     return [LayerResult(r) for r in raw_results]
 
+    return [LayerResult(r) for r in raw_results]
+
 
 # Alias for convenience
 slice_mesh = slice_mesh_and_offset
@@ -431,10 +451,14 @@ class SdfStack:
         export_infill_mask: bool = True,
         simplify_tolerance_mm: float = 0.02,
         subpixel_lerp: bool = True,
-        offset_distances: Optional[Sequence[float]] = None
+        offset_distances: Optional[Sequence[float]] = None,
+        enable_arachne: bool = False,
+        min_feature_size_mm: float = 0.20,
+        min_bead_width_mm: float = 0.25,
+        max_bead_width_mm: float = 0.65
     ) -> List[LayerResult]:
         """
-        Vectorize GPU SDF stack directly to smooth polygon perimeter contours and infill masks.
+        Vectorize GPU SDF stack directly to smooth polygon perimeter contours, variable-width Arachne beads, and infill masks.
 
         Parameters
         ----------
@@ -450,6 +474,14 @@ class SdfStack:
             Continuous SDF linear interpolation for sub-pixel boundary placement at arbitrary angles.
         offset_distances : sequence of float, optional
             Custom offset distances.
+        enable_arachne : bool, default=False
+            Enable Medial Axis Transform (MAT) variable line width beads for collapsing walls.
+        min_feature_size_mm : float, default=0.20
+            Minimum printable feature thickness in millimeters.
+        min_bead_width_mm : float, default=0.25
+            Minimum variable extrusion bead width in millimeters.
+        max_bead_width_mm : float, default=0.65
+            Maximum variable extrusion bead width in millimeters.
         """
         offsets_float = [float(d) for d in offset_distances] if offset_distances is not None else []
         raw_res = self._stack.vectorize(
@@ -458,7 +490,12 @@ class SdfStack:
             bool(export_infill_mask),
             float(simplify_tolerance_mm),
             bool(subpixel_lerp),
-            offsets_float
+            offsets_float,
+            False,
+            bool(enable_arachne),
+            float(min_feature_size_mm),
+            float(min_bead_width_mm),
+            float(max_bead_width_mm)
         )
         return [LayerResult(r) for r in raw_res]
 
@@ -560,6 +597,73 @@ def masks_to_sdf(
     return SdfStack(backend_sdf)
 
 
+def extract_arachne_mat(
+    sdf_or_stack: Any,
+    pitch_mm: float = 0.10,
+    min_x_mm: float = 0.0,
+    min_y_mm: float = 0.0,
+    nominal_line_width_mm: float = 0.45,
+    min_feature_size_mm: float = 0.20,
+    min_bead_width_mm: float = 0.25,
+    max_bead_width_mm: float = 0.65,
+    simplify_tolerance_mm: float = 0.02
+) -> Union[List[Dict[str, Any]], List[List[Dict[str, Any]]]]:
+    """
+    Extract Medial Axis Transform (MAT) variable line width beads directly from
+    a Signed Distance Field (SDF) 2D/3D numpy array or SdfStack.
+
+    Parameters
+    ----------
+    sdf_or_stack : SdfStack or np.ndarray
+        2D (H, W) or 3D (N, H, W) continuous float32 Signed Distance Field array or SdfStack.
+    pitch_mm : float, default=0.10
+        Pixel grid pitch in millimeters.
+    min_x_mm : float, default=0.0
+        Origin X in millimeters.
+    min_y_mm : float, default=0.0
+        Origin Y in millimeters.
+    nominal_line_width_mm : float, default=0.45
+        Standard nozzle extrusion line width in millimeters.
+    min_feature_size_mm : float, default=0.20
+        Minimum printable wall thickness in millimeters.
+    min_bead_width_mm : float, default=0.25
+        Minimum variable extrusion bead width in millimeters.
+    max_bead_width_mm : float, default=0.65
+        Maximum variable extrusion bead width in millimeters.
+    simplify_tolerance_mm : float, default=0.02
+        RDP polyline simplification tolerance in millimeters.
+
+    Returns
+    -------
+    list of dict or list of list of dict
+        Variable-width beads with 'points' (K, 2), 'widths' (K,), 'is_closed', and 'num_points'.
+    """
+    if not _HAS_BACKEND or _backend is None:
+        raise RuntimeError("CUDASlicer backend extension (_cudaslicer_backend) is not available.")
+
+    if isinstance(sdf_or_stack, SdfStack):
+        arr = sdf_or_stack.numpy()
+        meta = sdf_or_stack.metadata
+        p_mm = float(meta["pitch_um"]) * 0.001
+        mx_mm = float(meta["min_x_um"]) * 0.001
+        my_mm = float(meta["min_y_um"]) * 0.001
+    else:
+        arr = np.ascontiguousarray(sdf_or_stack, dtype=np.float32)
+        p_mm = float(pitch_mm)
+        mx_mm = float(min_x_mm)
+        my_mm = float(min_y_mm)
+
+    return _backend.extract_arachne_mat(
+        arr, p_mm, mx_mm, my_mm,
+        float(nominal_line_width_mm), float(min_feature_size_mm),
+        float(min_bead_width_mm), float(max_bead_width_mm),
+        float(simplify_tolerance_mm)
+    )
+
+
+arachne_mat = extract_arachne_mat
+
+
 __all__ = [
     "offset_masks",
     "offset_polygons",
@@ -567,6 +671,8 @@ __all__ = [
     "slice_mesh",
     "slice_mesh_to_sdf",
     "masks_to_sdf",
+    "extract_arachne_mat",
+    "arachne_mat",
     "SdfStack",
     "LayerResult",
     "is_available",

python/cudaslicer/bindings.cpp

@@ -2,6 +2,7 @@
 #include <pybind11/numpy.h>
 #include <pybind11/stl.h>
 #include "gpu_offset.h"
+#include "arachne_mat.h"
 #include <vector>
 #include <iostream>
 #include <string>
@@ -112,7 +113,12 @@ py::list py_offset_polygon_layers(
     bool export_infill_mask,
     const std::vector<float>& offset_distances,
     const std::string& backend,
-    bool output_3d = false
+    bool output_3d = false,
+    bool enable_arachne = false,
+    float min_feature_size_mm = 0.20f,
+    float min_bead_width_mm = 0.25f,
+    float max_bead_width_mm = 0.65f,
+    float simplify_tolerance_mm = 0.02f
 ) {
     size_t num_layers = py_layers.size();
     if (num_layers == 0) return py::list();
@@ -161,6 +167,11 @@ py::list py_offset_polygon_layers(
     config.export_infill_mask = export_infill_mask ? 1 : 0;
     config.use_gpu_rasterizer = 1;
     config.backend = parse_backend_str(backend);
+    config.enable_arachne = enable_arachne ? 1 : 0;
+    config.min_feature_size_mm = min_feature_size_mm;
+    config.min_bead_width_mm = min_bead_width_mm;
+    config.max_bead_width_mm = max_bead_width_mm;
+    config.simplify_tolerance_mm = simplify_tolerance_mm;
 
     if (!offset_distances.empty()) {
         config.offset_distances_mm = offset_distances.data();
@@ -271,6 +282,29 @@ py::list py_offset_polygon_layers(
             py_layer_dict["infill_mask"] = mask_arr;
         }
 
+        // Arachne Variable-Width Beads
+        py::list py_arachne;
+        for (size_t b = 0; b < results[l].num_arachne_beads; ++b) {
+            const auto& bead = results[l].arachne_beads[b];
+            py::dict py_bead;
+            py_bead["is_closed"] = (bead.is_closed != 0);
+            py_bead["num_points"] = bead.num_points;
+
+            py::array_t<float> pts_arr({ (ssize_t)bead.num_points, (ssize_t)2 });
+            py::array_t<float> widths_arr((ssize_t)bead.num_points);
+            auto r_pts = pts_arr.mutable_unchecked<2>();
+            auto r_w = widths_arr.mutable_unchecked<1>();
+            for (size_t p = 0; p < bead.num_points; ++p) {
+                r_pts(p, 0) = bead.points[p].x;
+                r_pts(p, 1) = bead.points[p].y;
+                r_w(p) = bead.widths[p];
+            }
+            py_bead["points"] = pts_arr;
+            py_bead["widths"] = widths_arr;
+            py_arachne.append(py_bead);
+        }
+        py_layer_dict["arachne_beads"] = py_arachne;
+
         py_layer_results.append(py_layer_dict);
     }
 
@@ -294,7 +328,11 @@ py::list py_slice_mesh_and_offset(
     const std::string& backend,
     float simplify_tolerance_mm,
     bool subpixel_lerp,
-    bool output_3d = false
+    bool output_3d = false,
+    bool enable_arachne = false,
+    float min_feature_size_mm = 0.20f,
+    float min_bead_width_mm = 0.25f,
+    float max_bead_width_mm = 0.65f
 ) {
     py::buffer_info buf = triangles_arr.request();
     size_t num_tris = 0;
@@ -321,6 +359,10 @@ py::list py_slice_mesh_and_offset(
     config.subpixel_lerp = subpixel_lerp ? 1 : 0;
     config.use_gpu_rasterizer = 1;
     config.backend = parse_backend_str(backend);
+    config.enable_arachne = enable_arachne ? 1 : 0;
+    config.min_feature_size_mm = min_feature_size_mm;
+    config.min_bead_width_mm = min_bead_width_mm;
+    config.max_bead_width_mm = max_bead_width_mm;
 
     if (!offset_distances.empty()) {
         config.offset_distances_mm = offset_distances.data();
@@ -433,6 +475,29 @@ py::list py_slice_mesh_and_offset(
             py_layer_dict["infill_mask"] = mask_arr;
         }
 
+        // Arachne Variable-Width Beads
+        py::list py_arachne;
+        for (size_t b = 0; b < results[l].num_arachne_beads; ++b) {
+            const auto& bead = results[l].arachne_beads[b];
+            py::dict py_bead;
+            py_bead["is_closed"] = (bead.is_closed != 0);
+            py_bead["num_points"] = bead.num_points;
+
+            py::array_t<float> pts_arr({ (ssize_t)bead.num_points, (ssize_t)2 });
+            py::array_t<float> widths_arr((ssize_t)bead.num_points);
+            auto r_pts = pts_arr.mutable_unchecked<2>();
+            auto r_w = widths_arr.mutable_unchecked<1>();
+            for (size_t p = 0; p < bead.num_points; ++p) {
+                r_pts(p, 0) = bead.points[p].x;
+                r_pts(p, 1) = bead.points[p].y;
+                r_w(p) = bead.widths[p];
+            }
+            py_bead["points"] = pts_arr;
+            py_bead["widths"] = widths_arr;
+            py_arachne.append(py_bead);
+        }
+        py_layer_dict["arachne_beads"] = py_arachne;
+
         py_layer_results.append(py_layer_dict);
     }
 
@@ -590,7 +655,11 @@ public:
         float simplify_tolerance_mm = 0.02f,
         bool subpixel_lerp = true,
         const std::vector<float>& offset_distances = std::vector<float>(),
-        bool output_3d = false
+        bool output_3d = false,
+        bool enable_arachne = false,
+        float min_feature_size_mm = 0.20f,
+        float min_bead_width_mm = 0.25f,
+        float max_bead_width_mm = 0.65f
     ) const {
         GpuOffsetConfig config = {0};
         config.perimeters_count = perimeters_count;
@@ -600,6 +669,10 @@ public:
         config.subpixel_lerp = subpixel_lerp ? 1 : 0;
         config.use_gpu_rasterizer = 1;
         config.backend = GPU_OFFSET_BACKEND_CUDA;
+        config.enable_arachne = enable_arachne ? 1 : 0;
+        config.min_feature_size_mm = min_feature_size_mm;
+        config.min_bead_width_mm = min_bead_width_mm;
+        config.max_bead_width_mm = max_bead_width_mm;
 
         if (!offset_distances.empty()) {
             config.offset_distances_mm = offset_distances.data();
@@ -698,6 +771,29 @@ public:
                 py_layer_dict["infill_mask"] = mask_arr;
             }
 
+            // Arachne Variable-Width Beads
+            py::list py_arachne;
+            for (size_t b = 0; b < results[l].num_arachne_beads; ++b) {
+                const auto& bead = results[l].arachne_beads[b];
+                py::dict py_bead;
+                py_bead["is_closed"] = (bead.is_closed != 0);
+                py_bead["num_points"] = bead.num_points;
+
+                py::array_t<float> pts_arr({ (ssize_t)bead.num_points, (ssize_t)2 });
+                py::array_t<float> widths_arr((ssize_t)bead.num_points);
+                auto r_pts = pts_arr.mutable_unchecked<2>();
+                auto r_w = widths_arr.mutable_unchecked<1>();
+                for (size_t p = 0; p < bead.num_points; ++p) {
+                    r_pts(p, 0) = bead.points[p].x;
+                    r_pts(p, 1) = bead.points[p].y;
+                    r_w(p) = bead.widths[p];
+                }
+                py_bead["points"] = pts_arr;
+                py_bead["widths"] = widths_arr;
+                py_arachne.append(py_bead);
+            }
+            py_layer_dict["arachne_beads"] = py_arachne;
+
             py_layer_results.append(py_layer_dict);
         }
 
@@ -797,8 +893,80 @@ std::shared_ptr<PyGpuSdfStack> py_masks_to_sdf(
     return std::make_shared<PyGpuSdfStack>(handle);
 }
 
+// Arachne Medial Axis Extraction on NumPy 2D/3D float SDF arrays
+py::list py_extract_arachne_mat(
+    py::array_t<float, py::array::c_style | py::array::forcecast> sdf_arr,
+    float pitch_mm = 0.10f,
+    float min_x_mm = 0.0f,
+    float min_y_mm = 0.0f,
+    float nominal_line_width_mm = 0.45f,
+    float min_feature_size_mm = 0.20f,
+    float min_bead_width_mm = 0.25f,
+    float max_bead_width_mm = 0.65f,
+    float simplify_tolerance_mm = 0.02f
+) {
+    py::buffer_info buf = sdf_arr.request();
+    int H = 0, W = 0, N = 1;
+    if (buf.ndim == 2) {
+        H = (int)buf.shape[0];
+        W = (int)buf.shape[1];
+        N = 1;
+    } else if (buf.ndim == 3) {
+        N = (int)buf.shape[0];
+        H = (int)buf.shape[1];
+        W = (int)buf.shape[2];
+    } else {
+        throw std::runtime_error("SDF array must be 2D (H, W) or 3D (N, H, W) float32.");
+    }
+
+    const float* raw_sdf = static_cast<const float*>(buf.ptr);
+    size_t layer_pixels = (size_t)H * W;
+
+    py::list all_layer_beads;
+    for (int l = 0; l < N; ++l) {
+        const float* l_sdf = raw_sdf + (size_t)l * layer_pixels;
+        std::vector<ArachneBeadInternal> internal_beads;
+        extract_arachne_mat_layer(
+            l_sdf, W, H, min_x_mm, min_y_mm, pitch_mm,
+            nominal_line_width_mm,
+            min_feature_size_mm,
+            min_bead_width_mm,
+            max_bead_width_mm,
+            simplify_tolerance_mm,
+            internal_beads
+        );
+
+        py::list py_bead_list;
+        for (const auto& bead : internal_beads) {
+            py::dict py_bead;
+            py_bead["is_closed"] = (bead.is_closed != 0);
+            py_bead["num_points"] = bead.points.size();
+
+            py::array_t<float> pts_arr({ (ssize_t)bead.points.size(), (ssize_t)2 });
+            py::array_t<float> widths_arr((ssize_t)bead.widths.size());
+            auto r_pts = pts_arr.mutable_unchecked<2>();
+            auto r_w = widths_arr.mutable_unchecked<1>();
+            for (size_t p = 0; p < bead.points.size(); ++p) {
+                r_pts(p, 0) = bead.points[p].x;
+                r_pts(p, 1) = bead.points[p].y;
+                r_w(p) = bead.widths[p];
+            }
+            py_bead["points"] = pts_arr;
+            py_bead["widths"] = widths_arr;
+            py_bead_list.append(py_bead);
+        }
+
+        if (buf.ndim == 2) {
+            return py_bead_list;
+        } else {
+            all_layer_beads.append(py_bead_list);
+        }
+    }
+    return all_layer_beads;
+}
+
 PYBIND11_MODULE(_cudaslicer_backend, m) {
-    m.doc() = "GPU/CPU Slicing, Signed Distance Fields (SDF), CSG Booleans & Vectorizer Extension";
+    m.doc() = "GPU/CPU Slicing, Signed Distance Fields (SDF), CSG Booleans, Vectorizer & Arachne MAT Extension";
 
     m.def("has_cuda", &gpu_offset_has_cuda, "Returns 1 if a functional CUDA GPU is present, 0 otherwise.");
     m.def("is_cuda_available", &gpu_offset_has_cuda, "Returns 1 if a functional CUDA GPU is present, 0 otherwise.");
@@ -838,6 +1006,11 @@ PYBIND11_MODULE(_cudaslicer_backend, m) {
           py::arg("offset_distances") = std::vector<float>(),
           py::arg("backend") = "auto",
           py::arg("output_3d") = false,
+          py::arg("enable_arachne") = false,
+          py::arg("min_feature_size_mm") = 0.20f,
+          py::arg("min_bead_width_mm") = 0.25f,
+          py::arg("max_bead_width_mm") = 0.65f,
+          py::arg("simplify_tolerance_mm") = 0.02f,
           "Compute polygon perimeter offsets and infill mask.");
 
     m.def("slice_mesh_and_offset", &py_slice_mesh_and_offset,
@@ -855,8 +1028,24 @@ PYBIND11_MODULE(_cudaslicer_backend, m) {
           py::arg("simplify_tolerance_mm") = 0.0f,
           py::arg("subpixel_lerp") = true,
           py::arg("output_3d") = false,
+          py::arg("enable_arachne") = false,
+          py::arg("min_feature_size_mm") = 0.20f,
+          py::arg("min_bead_width_mm") = 0.25f,
+          py::arg("max_bead_width_mm") = 0.65f,
           "Mesh slicing, perimeter offsetting, and infill mask generation.");
 
+    m.def("extract_arachne_mat", &py_extract_arachne_mat,
+          py::arg("sdf_arr"),
+          py::arg("pitch_mm") = 0.10f,
+          py::arg("min_x_mm") = 0.0f,
+          py::arg("min_y_mm") = 0.0f,
+          py::arg("nominal_line_width_mm") = 0.45f,
+          py::arg("min_feature_size_mm") = 0.20f,
+          py::arg("min_bead_width_mm") = 0.25f,
+          py::arg("max_bead_width_mm") = 0.65f,
+          py::arg("simplify_tolerance_mm") = 0.02f,
+          "Extract Medial Axis Transform (MAT) variable line width beads directly from 2D/3D float SDF arrays.");
+
     // GPU SDF Stack Factory Functions
     m.def("slice_mesh_to_sdf", &py_slice_mesh_to_sdf,
           py::arg("triangles_arr"), py::arg("min_x"), py::arg("min_y"), py::arg("min_z"),
@@ -901,7 +1090,11 @@ PYBIND11_MODULE(_cudaslicer_backend, m) {
              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.")
+             py::arg("enable_arachne") = false,
+             py::arg("min_feature_size_mm") = 0.20f,
+             py::arg("min_bead_width_mm") = 0.25f,
+             py::arg("max_bead_width_mm") = 0.65f,
+             "Vectorize SDF stack directly with sub-pixel LERP, RDP polygon simplification, and optional Arachne MAT.")
         .def_property_readonly("__cuda_array_interface__", &PyGpuSdfStack::cuda_array_interface,
                                "Zero-copy CUDA array interface for PyTorch / CuPy integration.");
 }

src/arachne_mat.h

@@ -0,0 +1,416 @@
+#ifndef ARACHNE_MAT_H
+#define ARACHNE_MAT_H
+
+#include "gpu_offset.h"
+#include <vector>
+#include <cmath>
+#include <algorithm>
+#include <cstring>
+#include <cstdlib>
+
+struct ArachneRidgeNode {
+    int id;
+    int gx, gy;
+    float px, py;
+    float width;
+    int dir;
+};
+
+struct ArachneBeadInternal {
+    std::vector<GpuPoint2D> points;
+    std::vector<float> widths;
+    int is_closed;
+};
+
+// RDP simplification with per-vertex width retention
+static inline void rdp_simplify_bead_recursive(
+    const std::vector<GpuPoint2D>& pts,
+    const std::vector<float>& widths,
+    size_t i, size_t j,
+    float eps_sq,
+    std::vector<uint8_t>& keep
+) {
+    if (j <= i + 1) return;
+
+    const GpuPoint2D& p1 = pts[i];
+    const GpuPoint2D& p2 = pts[j];
+    float dx = p2.x - p1.x;
+    float dy = p2.y - p1.y;
+    float len_sq = dx * dx + dy * dy;
+
+    float max_dist_sq = 0.0f;
+    size_t max_idx = i;
+
+    for (size_t k = i + 1; k < j; ++k) {
+        const GpuPoint2D& pk = pts[k];
+        float dist_sq = 0.0f;
+        if (len_sq < 1e-10f) {
+            float dpx = pk.x - p1.x;
+            float dpy = pk.y - p1.y;
+            dist_sq = dpx * dpx + dpy * dpy;
+        } else {
+            float cross = dy * (pk.x - p1.x) - dx * (pk.y - p1.y);
+            dist_sq = (cross * cross) / len_sq;
+        }
+
+        if (dist_sq > max_dist_sq) {
+            max_dist_sq = dist_sq;
+            max_idx = k;
+        }
+    }
+
+    if (max_dist_sq > eps_sq) {
+        keep[max_idx] = 1;
+        rdp_simplify_bead_recursive(pts, widths, i, max_idx, eps_sq, keep);
+        rdp_simplify_bead_recursive(pts, widths, max_idx, j, eps_sq, keep);
+    }
+}
+
+static inline void simplify_bead_rdp(
+    ArachneBeadInternal& bead,
+    float tolerance_mm
+) {
+    size_t n = bead.points.size();
+    if (n <= 2 || tolerance_mm <= 1e-6f) return;
+
+    float eps_sq = tolerance_mm * tolerance_mm;
+    std::vector<uint8_t> keep(n, 0);
+    keep[0] = 1;
+    keep[n - 1] = 1;
+
+    rdp_simplify_bead_recursive(bead.points, bead.widths, 0, n - 1, eps_sq, keep);
+
+    std::vector<GpuPoint2D> simp_pts;
+    std::vector<float> simp_w;
+    simp_pts.reserve(n);
+    simp_w.reserve(n);
+
+    for (size_t i = 0; i < n; ++i) {
+        if (keep[i]) {
+            simp_pts.push_back(bead.points[i]);
+            simp_w.push_back(bead.widths[i]);
+        }
+    }
+
+    bead.points = std::move(simp_pts);
+    bead.widths = std::move(simp_w);
+}
+
+/**
+ * @brief Extract Medial Axis Transform (MAT) variable line width beads from continuous 2D Signed Distance Field (SDF).
+ */
+static inline void extract_arachne_mat_layer(
+    const float* sdf,
+    int width, int height,
+    float min_x, float min_y, float pitch_mm,
+    float nominal_line_width_mm,
+    float min_feature_size_mm,
+    float min_bead_width_mm,
+    float max_bead_width_mm,
+    float simplify_tolerance_mm,
+    std::vector<ArachneBeadInternal>& out_beads
+) {
+    if (!sdf || width < 3 || height < 3) return;
+
+    float W_nom = (nominal_line_width_mm > 0.0f) ? nominal_line_width_mm : 0.45f;
+    float T_min = (min_feature_size_mm > 0.0f) ? min_feature_size_mm : 0.20f;
+    float W_min = (min_bead_width_mm > 0.0f) ? min_bead_width_mm : 0.25f;
+    float W_max = (max_bead_width_mm > 0.0f) ? max_bead_width_mm : 0.65f;
+
+    // Collapsing region upper and lower bounds in distance (half-thickness)
+    float d_min = 0.5f * T_min;
+    float d_max = W_nom + 0.5f * pitch_mm;
+
+    std::vector<int> node_grid(width * height, -1);
+    std::vector<ArachneRidgeNode> nodes;
+    nodes.reserve(width * height / 32);
+
+    // 1. Detect Medial Axis Ridges on 2D continuous SDF
+    for (int y = 1; y < height - 1; ++y) {
+        const float* row_prev = sdf + (y - 1) * width;
+        const float* row_cur  = sdf + y * width;
+        const float* row_next = sdf + (y + 1) * width;
+
+        for (int x = 1; x < width - 1; ++x) {
+            float c = row_cur[x];
+            if (c < d_min || c > d_max) continue;
+
+            float c_l  = row_cur[x - 1];
+            float c_r  = row_cur[x + 1];
+            float c_d  = row_prev[x];
+            float c_u  = row_next[x];
+            float c_dl = row_prev[x - 1];
+            float c_ur = row_next[x + 1];
+            float c_ul = row_next[x - 1];
+            float c_dr = row_prev[x + 1];
+
+            float c_ll = (x >= 2) ? row_cur[x - 2] : c_l;
+            float c_rr = (x < width - 2) ? row_cur[x + 2] : c_r;
+            float c_dd = (y >= 2) ? sdf[(y - 2) * width + x] : c_d;
+            float c_uu = (y < height - 2) ? sdf[(y + 2) * width + x] : c_u;
+
+            float D_xx = c_l - 2.0f * c + c_r;
+            float D_yy = c_d - 2.0f * c + c_u;
+            float D_d1 = c_dl - 2.0f * c + c_ur;
+            float D_d2 = c_ul - 2.0f * c + c_dr;
+
+            float curv_thresh = -0.05f * pitch_mm;
+
+            bool is_peak_h = (c > c_l && c > c_r && D_xx <= curv_thresh) ||
+                             (c == c_r && c > c_l && c_r > c_rr && D_xx <= curv_thresh) ||
+                             (c == c_l && c > c_r && c_l > c_ll && D_xx <= curv_thresh);
+            bool is_peak_v = (c > c_d && c > c_u && D_yy <= curv_thresh) ||
+                             (c == c_u && c > c_d && c_u > c_uu && D_yy <= curv_thresh) ||
+                             (c == c_d && c > c_u && c_d > c_dd && D_yy <= curv_thresh);
+            bool is_peak_d1 = (c > c_dl && c > c_ur && D_d1 <= curv_thresh);
+            bool is_peak_d2 = (c > c_ul && c > c_dr && D_d2 <= curv_thresh);
+
+            if (!is_peak_h && !is_peak_v && !is_peak_d1 && !is_peak_d2) continue;
+
+            // Pick direction of strongest ridge crest
+            float min_curv = curv_thresh;
+            int best_dir = -1;
+            if (is_peak_h && D_xx < min_curv) { min_curv = D_xx; best_dir = 0; }
+            if (is_peak_v && D_yy < min_curv) { min_curv = D_yy; best_dir = 1; }
+            if (is_peak_d1 && D_d1 < min_curv) { min_curv = D_d1; best_dir = 2; }
+            if (is_peak_d2 && D_d2 < min_curv) { min_curv = D_d2; best_dir = 3; }
+
+            if (best_dir == -1) continue;
+
+            float sub_px = (float)x;
+            float sub_py = (float)y;
+            float d_peak = c;
+
+            if (best_dir == 0) {
+                float denom = 2.0f * (c_l - 2.0f * c + c_r);
+                float delta = (std::abs(denom) > 1e-7f) ? (c_l - c_r) / denom : 0.0f;
+                delta = std::max(-0.5f, std::min(0.5f, delta));
+                sub_px += delta;
+                d_peak = c - ((c_l - c_r) * (c_l - c_r)) / (8.0f * (c_l - 2.0f * c + c_r));
+            } else if (best_dir == 1) {
+                float denom = 2.0f * (c_d - 2.0f * c + c_u);
+                float delta = (std::abs(denom) > 1e-7f) ? (c_d - c_u) / denom : 0.0f;
+                delta = std::max(-0.5f, std::min(0.5f, delta));
+                sub_py += delta;
+                d_peak = c - ((c_d - c_u) * (c_d - c_u)) / (8.0f * (c_d - 2.0f * c + c_u));
+            } else if (best_dir == 2) {
+                float denom = 2.0f * (c_dl - 2.0f * c + c_ur);
+                float delta = (std::abs(denom) > 1e-7f) ? (c_dl - c_ur) / denom : 0.0f;
+                delta = std::max(-0.5f, std::min(0.5f, delta));
+                sub_px += delta;
+                sub_py += delta;
+                d_peak = c - ((c_dl - c_ur) * (c_dl - c_ur)) / (8.0f * (c_dl - 2.0f * c + c_ur));
+            } else if (best_dir == 3) {
+                float denom = 2.0f * (c_ul - 2.0f * c + c_dr);
+                float delta = (std::abs(denom) > 1e-7f) ? (c_ul - c_dr) / denom : 0.0f;
+                delta = std::max(-0.5f, std::min(0.5f, delta));
+                sub_px += delta;
+                sub_py -= delta;
+                d_peak = c - ((c_ul - c_dr) * (c_ul - c_dr)) / (8.0f * (c_ul - 2.0f * c + c_dr));
+            }
+
+            d_peak = std::max(d_peak, c);
+            float thickness = 2.0f * d_peak;
+            float bead_w = std::min(std::max(thickness, W_min), W_max);
+
+            int node_id = (int)nodes.size();
+            nodes.push_back({
+                node_id,
+                x, y,
+                min_x + sub_px * pitch_mm,
+                min_y + sub_py * pitch_mm,
+                bead_w,
+                best_dir
+            });
+            node_grid[y * width + x] = node_id;
+        }
+    }
+
+    if (nodes.empty()) return;
+
+    // 2. Build Adjacency Graph of Ridge Skeleton with Tangent Compatibility
+    size_t num_nodes = nodes.size();
+    std::vector<std::vector<int>> adj(num_nodes);
+
+    const int dx8[8] = { 1, 1, 0, -1, -1, -1,  0,  1 };
+    const int dy8[8] = { 0, 1, 1,  1,  0, -1, -1, -1 };
+
+    for (size_t i = 0; i < num_nodes; ++i) {
+        int gx = nodes[i].gx;
+        int gy = nodes[i].gy;
+        int bdir = nodes[i].dir;
+        for (int k = 0; k < 8; ++k) {
+            int dx = dx8[k];
+            int dy = dy8[k];
+            // Directional tangent filter:
+            // best_dir == 0: cross-curvature in X -> ridge runs along Y -> dy != 0
+            if (bdir == 0 && dy == 0) continue;
+            // best_dir == 1: cross-curvature in Y -> ridge runs along X -> dx != 0
+            if (bdir == 1 && dx == 0) continue;
+            // best_dir == 2: cross-curvature in DL-UR -> ridge runs along UL-DR -> reject dx == dy
+            if (bdir == 2 && dx == dy) continue;
+            // best_dir == 3: cross-curvature in UL-DR -> ridge runs along DL-UR -> reject dx == -dy
+            if (bdir == 3 && dx == -dy) continue;
+
+            int nx = gx + dx;
+            int ny = gy + dy;
+            if (nx >= 0 && nx < width && ny >= 0 && ny < height) {
+                int neighbor_id = node_grid[ny * width + nx];
+                if (neighbor_id > (int)i) {
+                    adj[i].push_back(neighbor_id);
+                    adj[neighbor_id].push_back((int)i);
+                }
+            }
+        }
+    }
+
+    // 3. Trace Graph Chains into Continuous Polylines
+    std::vector<uint8_t> visited_nodes(num_nodes, 0);
+    std::vector<std::pair<int, int>> visited_edges; // store traversed (u, v)
+
+    auto is_edge_visited = [&](int u, int v) {
+        int a = std::min(u, v);
+        int b = std::max(u, v);
+        for (const auto& e : visited_edges) {
+            if (e.first == a && e.second == b) return true;
+        }
+        return false;
+    };
+
+    auto mark_edge_visited = [&](int u, int v) {
+        visited_edges.push_back({ std::min(u, v), std::max(u, v) });
+    };
+
+    // Priority 1: Trace from Endpoints (degree 1) or Junctions (degree >= 3)
+    for (size_t start_pass = 0; start_pass < 2; ++start_pass) {
+        for (size_t i = 0; i < num_nodes; ++i) {
+            size_t deg = adj[i].size();
+            if (deg == 0) continue;
+
+            if (start_pass == 0 && deg == 2) continue; // First pass: only endpoints/junctions
+
+            for (int neighbor : adj[i]) {
+                if (is_edge_visited((int)i, neighbor)) continue;
+
+                std::vector<GpuPoint2D> poly_pts;
+                std::vector<float> poly_w;
+
+                poly_pts.push_back({ nodes[i].px, nodes[i].py });
+                poly_w.push_back(nodes[i].width);
+
+                int prev = (int)i;
+                int cur = neighbor;
+                mark_edge_visited(prev, cur);
+
+                poly_pts.push_back({ nodes[cur].px, nodes[cur].py });
+                poly_w.push_back(nodes[cur].width);
+
+                // Follow chain of degree-2 nodes
+                while (adj[cur].size() == 2 && cur != (int)i) {
+                    int next = (adj[cur][0] == prev) ? adj[cur][1] : adj[cur][0];
+                    if (is_edge_visited(cur, next)) break;
+
+                    mark_edge_visited(cur, next);
+                    prev = cur;
+                    cur = next;
+
+                    poly_pts.push_back({ nodes[cur].px, nodes[cur].py });
+                    poly_w.push_back(nodes[cur].width);
+                }
+
+                // Filter out tiny micro-segments & noise
+                size_t n_pts = poly_pts.size();
+                if (n_pts < 2) continue;
+
+                float arc_len = 0.0f;
+                for (size_t k = 0; k < n_pts - 1; ++k) {
+                    float dx = poly_pts[k + 1].x - poly_pts[k].x;
+                    float dy = poly_pts[k + 1].y - poly_pts[k].y;
+                    arc_len += std::sqrt(dx * dx + dy * dy);
+                }
+
+                float min_bead_len = std::max(1.2f * W_nom, 0.50f);
+                if (arc_len < min_bead_len) continue; // Drop micro-slivers and short corner spurs
+
+                // Apply moving average smoothing (preserve endpoints)
+                if (n_pts >= 3) {
+                    std::vector<GpuPoint2D> smooth_pts = poly_pts;
+                    std::vector<float> smooth_w = poly_w;
+                    for (size_t k = 1; k < n_pts - 1; ++k) {
+                        smooth_pts[k].x = 0.25f * poly_pts[k - 1].x + 0.5f * poly_pts[k].x + 0.25f * poly_pts[k + 1].x;
+                        smooth_pts[k].y = 0.25f * poly_pts[k - 1].y + 0.5f * poly_pts[k].y + 0.25f * poly_pts[k + 1].y;
+                        smooth_w[k] = 0.25f * poly_w[k - 1] + 0.5f * poly_w[k] + 0.25f * poly_w[k + 1];
+                    }
+                    poly_pts = std::move(smooth_pts);
+                    poly_w = std::move(smooth_w);
+                }
+
+                ArachneBeadInternal bead;
+                bead.points = std::move(poly_pts);
+                bead.widths = std::move(poly_w);
+                bead.is_closed = (cur == (int)i) ? 1 : 0;
+
+                if (simplify_tolerance_mm > 1e-6f) {
+                    simplify_bead_rdp(bead, simplify_tolerance_mm);
+                }
+
+                if (bead.points.size() >= 2) {
+                    out_beads.push_back(std::move(bead));
+                }
+            }
+        }
+    }
+
+    // 4. Stitch adjacent open polylines that connect at endpoints
+    if (out_beads.size() >= 2) {
+        std::vector<ArachneBeadInternal> merged_beads;
+        std::vector<uint8_t> used(out_beads.size(), 0);
+
+        for (size_t i = 0; i < out_beads.size(); ++i) {
+            if (used[i]) continue;
+            used[i] = 1;
+
+            ArachneBeadInternal cur = std::move(out_beads[i]);
+            if (cur.is_closed) {
+                merged_beads.push_back(std::move(cur));
+                continue;
+            }
+
+            bool extended = true;
+            while (extended) {
+                extended = false;
+                const GpuPoint2D& tail = cur.points.back();
+
+                for (size_t j = 0; j < out_beads.size(); ++j) {
+                    if (used[j] || out_beads[j].is_closed) continue;
+
+                    const GpuPoint2D& j_head = out_beads[j].points.front();
+                    const GpuPoint2D& j_tail = out_beads[j].points.back();
+
+                    float d_tail_jhead = std::hypot(tail.x - j_head.x, tail.y - j_head.y);
+                    float d_tail_jtail = std::hypot(tail.x - j_tail.x, tail.y - j_tail.y);
+
+                    if (d_tail_jhead < 2.0f * pitch_mm) {
+                        used[j] = 1;
+                        extended = true;
+                        cur.points.insert(cur.points.end(), out_beads[j].points.begin() + 1, out_beads[j].points.end());
+                        cur.widths.insert(cur.widths.end(), out_beads[j].widths.begin() + 1, out_beads[j].widths.end());
+                        break;
+                    } else if (d_tail_jtail < 2.0f * pitch_mm) {
+                        used[j] = 1;
+                        extended = true;
+                        for (size_t k = out_beads[j].points.size() - 1; k-- > 0; ) {
+                            cur.points.push_back(out_beads[j].points[k]);
+                            cur.widths.push_back(out_beads[j].widths[k]);
+                        }
+                        break;
+                    }
+                }
+            }
+            merged_beads.push_back(std::move(cur));
+        }
+        out_beads = std::move(merged_beads);
+    }
+}
+
+#endif // ARACHNE_MAT_H

src/cpu_offset.cpp

@@ -360,10 +360,11 @@ int cpu_offset_polygon_layers(
         active_edges_buf.reserve(cols * rows / 16);
         std::vector<uint8_t> visited_buf(max_edges, 0);
 
+        bool need_sdf = subpixel_lerp || config->enable_arachne;
         std::vector<float> dist_sq_ext(layer_pixels);
         std::vector<float> dist_sq_int;
         std::vector<float> layer_sdf;
-        if (subpixel_lerp) {
+        if (need_sdf) {
             dist_sq_int.resize(layer_pixels);
             layer_sdf.resize(layer_pixels);
         }
@@ -383,7 +384,7 @@ int cpu_offset_polygon_layers(
             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) {
+            if (need_sdf) {
                 edt_2d_layer(mask_ptr, cols, rows, dist_sq_int.data(),
                              g_col_buf, row_f, row_d, v_buf, z_buf, true);
 
@@ -510,6 +511,39 @@ int cpu_offset_polygon_layers(
             } else {
                 res_array[l].infill.mask = nullptr;
             }
+
+            // Arachne Medial Axis Transform (MAT) variable line width beads
+            if (config->enable_arachne) {
+                std::vector<ArachneBeadInternal> arachne_beads_internal;
+                extract_arachne_mat_layer(
+                    layer_sdf.data(), cols, rows, bounds_min_x, bounds_min_y, pitch,
+                    config->layer_width,
+                    config->min_feature_size_mm,
+                    config->min_bead_width_mm,
+                    config->max_bead_width_mm,
+                    config->simplify_tolerance_mm,
+                    arachne_beads_internal
+                );
+
+                res_array[l].num_arachne_beads = arachne_beads_internal.size();
+                if (!arachne_beads_internal.empty()) {
+                    res_array[l].arachne_beads = (GpuBead*)malloc(arachne_beads_internal.size() * sizeof(GpuBead));
+                    for (size_t b = 0; b < arachne_beads_internal.size(); ++b) {
+                        const auto& src_bead = arachne_beads_internal[b];
+                        res_array[l].arachne_beads[b].num_points = src_bead.points.size();
+                        res_array[l].arachne_beads[b].is_closed = src_bead.is_closed;
+                        res_array[l].arachne_beads[b].points = (GpuPoint2D*)malloc(src_bead.points.size() * sizeof(GpuPoint2D));
+                        res_array[l].arachne_beads[b].widths = (float*)malloc(src_bead.widths.size() * sizeof(float));
+                        memcpy(res_array[l].arachne_beads[b].points, src_bead.points.data(), src_bead.points.size() * sizeof(GpuPoint2D));
+                        memcpy(res_array[l].arachne_beads[b].widths, src_bead.widths.data(), src_bead.widths.size() * sizeof(float));
+                    }
+                } else {
+                    res_array[l].arachne_beads = nullptr;
+                }
+            } else {
+                res_array[l].arachne_beads = nullptr;
+                res_array[l].num_arachne_beads = 0;
+            }
         }
     }
 
@@ -1001,6 +1035,52 @@ int cpu_offset_slice_mesh_and_offset(
                     res_array[l].infill.mask = nullptr;
                 }
             }
+
+            // Arachne Medial Axis Transform (MAT) variable line width beads
+            if (config->enable_arachne) {
+                const float* sdf_data = nullptr;
+                std::vector<float> temp_sdf;
+                if (subpixel_lerp) {
+                    sdf_data = layer_sdf.data();
+                } else {
+                    temp_sdf.resize(layer_pixels);
+                    for (size_t i = 0; i < layer_pixels; ++i) {
+                        float d = std::sqrt(dist_sq_ext[i]) * pitch;
+                        temp_sdf[i] = (layer_mask[i] > 0) ? d : -d;
+                    }
+                    sdf_data = temp_sdf.data();
+                }
+
+                std::vector<ArachneBeadInternal> arachne_beads_internal;
+                extract_arachne_mat_layer(
+                    sdf_data, cols, rows, bounds_min_x, bounds_min_y, pitch,
+                    config->layer_width,
+                    config->min_feature_size_mm,
+                    config->min_bead_width_mm,
+                    config->max_bead_width_mm,
+                    config->simplify_tolerance_mm,
+                    arachne_beads_internal
+                );
+
+                res_array[l].num_arachne_beads = arachne_beads_internal.size();
+                if (!arachne_beads_internal.empty()) {
+                    res_array[l].arachne_beads = (GpuBead*)malloc(arachne_beads_internal.size() * sizeof(GpuBead));
+                    for (size_t b = 0; b < arachne_beads_internal.size(); ++b) {
+                        const auto& src_bead = arachne_beads_internal[b];
+                        res_array[l].arachne_beads[b].num_points = src_bead.points.size();
+                        res_array[l].arachne_beads[b].is_closed = src_bead.is_closed;
+                        res_array[l].arachne_beads[b].points = (GpuPoint2D*)malloc(src_bead.points.size() * sizeof(GpuPoint2D));
+                        res_array[l].arachne_beads[b].widths = (float*)malloc(src_bead.widths.size() * sizeof(float));
+                        memcpy(res_array[l].arachne_beads[b].points, src_bead.points.data(), src_bead.points.size() * sizeof(GpuPoint2D));
+                        memcpy(res_array[l].arachne_beads[b].widths, src_bead.widths.data(), src_bead.widths.size() * sizeof(float));
+                    }
+                } else {
+                    res_array[l].arachne_beads = nullptr;
+                }
+            } else {
+                res_array[l].arachne_beads = nullptr;
+                res_array[l].num_arachne_beads = 0;
+            }
         }
     }
 

src/gpu_offset.cu

@@ -581,6 +581,10 @@ GPU_OFFSET_API int gpu_offset_polygon_layers(
         eff_backend = gpu_offset_has_cuda() ? GPU_OFFSET_BACKEND_CUDA : GPU_OFFSET_BACKEND_CPU;
     }
 
+    if (eff_backend == GPU_OFFSET_BACKEND_CPU || config->subpixel_lerp || config->enable_arachne) {
+        return cpu_offset_polygon_layers(in_layers, num_layers, min_x, min_y, max_x, max_y, config, out_results, stats);
+    }
+
     if (eff_backend == GPU_OFFSET_BACKEND_CUDA) {
         if (!gpu_offset_has_cuda()) {
             return GPU_OFFSET_ERR_NO_CUDA;
@@ -816,6 +820,13 @@ GPU_OFFSET_API void gpu_offset_free_results(GpuLayerResult* results, size_t num_
         if (results[l].infill.mask) {
             free(results[l].infill.mask);
         }
+        if (results[l].arachne_beads) {
+            for (size_t b = 0; b < results[l].num_arachne_beads; ++b) {
+                if (results[l].arachne_beads[b].points) free(results[l].arachne_beads[b].points);
+                if (results[l].arachne_beads[b].widths) free(results[l].arachne_beads[b].widths);
+            }
+            free(results[l].arachne_beads);
+        }
     }
     free(results);
 }

src/gpu_sdf.cu

@@ -14,27 +14,6 @@
 #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; \
@@ -1651,6 +1630,45 @@ GPU_SDF_API int gpu_sdf_vectorize(
         map.free_mem();
     }
 
+    // Arachne Medial Axis Transform (MAT) variable line width beads
+    if (config->enable_arachne) {
+        size_t layer_pixels = (size_t)cols * rows;
+        size_t total_floats = layer_pixels * num_layers;
+        std::vector<float> host_sdf(total_floats);
+        cudaMemcpy(host_sdf.data(), stack->d_sdf, total_floats * sizeof(float), cudaMemcpyDeviceToHost);
+
+        #pragma omp parallel for schedule(dynamic, 2)
+        for (int l = 0; l < (int)num_layers; ++l) {
+            const float* layer_sdf = host_sdf.data() + (size_t)l * layer_pixels;
+            std::vector<ArachneBeadInternal> arachne_beads_internal;
+            extract_arachne_mat_layer(
+                layer_sdf, cols, rows, bounds_min_x, bounds_min_y, pitch_mm,
+                config->layer_width,
+                config->min_feature_size_mm,
+                config->min_bead_width_mm,
+                config->max_bead_width_mm,
+                config->simplify_tolerance_mm,
+                arachne_beads_internal
+            );
+
+            res_array[l].num_arachne_beads = arachne_beads_internal.size();
+            if (!arachne_beads_internal.empty()) {
+                res_array[l].arachne_beads = (GpuBead*)malloc(arachne_beads_internal.size() * sizeof(GpuBead));
+                for (size_t b = 0; b < arachne_beads_internal.size(); ++b) {
+                    const auto& src_bead = arachne_beads_internal[b];
+                    res_array[l].arachne_beads[b].num_points = src_bead.points.size();
+                    res_array[l].arachne_beads[b].is_closed = src_bead.is_closed;
+                    res_array[l].arachne_beads[b].points = (GpuPoint2D*)malloc(src_bead.points.size() * sizeof(GpuPoint2D));
+                    res_array[l].arachne_beads[b].widths = (float*)malloc(src_bead.widths.size() * sizeof(float));
+                    memcpy(res_array[l].arachne_beads[b].points, src_bead.points.data(), src_bead.points.size() * sizeof(GpuPoint2D));
+                    memcpy(res_array[l].arachne_beads[b].widths, src_bead.widths.data(), src_bead.widths.size() * sizeof(float));
+                }
+            } else {
+                res_array[l].arachne_beads = nullptr;
+            }
+        }
+    }
+
     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();

src/vectorizer.h

@@ -2,6 +2,7 @@
 #define VECTORIZER_H
 
 #include "gpu_offset.h"
+#include "arachne_mat.h"
 #include <vector>
 #include <cmath>
 #include <cstring>
@@ -386,6 +387,17 @@ static inline void vectorize_sdf_layer(
         }
 
         if (loop_pts.size() >= 3) {
+            // Area check to discard micro-sliver artifacts at pinch-off points
+            double loop_area = 0.0;
+            size_t n_pts = loop_pts.size();
+            for (size_t i = 0; i < n_pts; ++i) {
+                size_t next_i = (i + 1 == n_pts) ? 0 : i + 1;
+                loop_area += (double)loop_pts[i].x * (double)loop_pts[next_i].y - (double)loop_pts[next_i].x * (double)loop_pts[i].y;
+            }
+            if (std::abs(loop_area) * 0.5 < 0.25 * pitch_mm * pitch_mm) {
+                continue; // Discard degenerate artifact loop
+            }
+
             std::vector<GpuPoint2D> simplified = (simplify_tolerance_mm > 1e-6f)
                 ? simplify_contour_rdp(loop_pts, simplify_tolerance_mm)
                 : loop_pts;