How To Create a 3d Printed Jigsaw Puzzle

I wanted to make a 3d printed jigsaw puzzle as a gift for my Mom who’s both a puzzle and map lover but it ended up being a little bit harder than I was expecting so here’s how I ended up doing it which hopefully will save you some trouble if you want to make one too!

There’s lots of different sites and ways to download map data into a model but the best one I found that lets you customize colors, details, and a whole lot more is Mini Skyline. It does even have a “project mode” that theoretically lets you subdivide the model into the puzzle pieces but as of the time of writing it’s still bugged out as far as I can tell (and with all the combinations of parameters I’ve tried) and won’t let you export. It also has a limit of 20 x 20 pieces which my method doesn’t have.

So to get started, the first step is to select the area you’re interested in, customize the colors to something close to what you can print, make sure to set the intended model longest side length (which informs the scale and will be used later). I exported with the terrain tile mode and high quality, but neither of those are required for this method.

The next step is to create your puzzle piece shapes using JigsawMake. Choose the layout option and then customize your piece count with the rows and columns and match the dimensions of the puzzle with the model you generated. Then use the export “Pieces SVG (Zip)” option which has each piece as its own svg file in a zip file. When you unzip it you should see a piece_r1_c1.svg file along with all the other pieces.

The next step is to get Blender set up for cutting the model into pieces. First go into the preferences > Add-ons and add the 3MF Import/Export plugin from Clonephaze. Delete the default cube and then from File > Import > Import 3MF file and select your model. Then I renamed my model in the hierarchy to something shorter (mt) since the script will need to reference it with the ORIGINAL_NAME.

After that go to the script tab and paste in the script below. The important things to configure are at the top. I would start by putting a TEST_LIMIT of 5 pieces, then run the script by clicking the play icon. Verify that you’re getting pieces generated and that they import into your slicer and look right after slicing. After you get it looking good, then set the TEST_LIMIT back to none and generate the whole puzzle. This script is vibecoded from Claude and hasn’t been optimized so it took quite awhile (probably 15-20min) to generate a 400 piece puzzle on a fairly beefy machine.

cut-pieces-blender-script.py

SVG_FOLDER = "/Users/you/Downloads/3D Prints/puzzle_pieces_svg_20x20"      # folder containing piece_000.svg etc.
ORIGINAL_NAME = "mt"          # name of your original colored object in the blender scene
OUTPUT_DIR = "/Users/you/Downloads/3D Prints/puzzle-out"      # folder to write piece_000.3mf, etc.
EXTRUDE_DEPTH = 100.0                     # mm; bigger than your model's thickness
TEST_LIMIT = None                           # set to None to process every piece
GRID_ROWS = 20                           # rows in your puzzle grid (e.g. 20x20=400)
GRID_COLS = 20                           # columns in your puzzle grid

# The config below shouldn't need to be touched
ONLY_FILES = None             # e.g. {"r1_c3.svg"} to re-cut specific pieces
# Fraction of a piece's footprint its upward-facing surface must cover.
# Border pieces of a non-rectangular model can legitimately fall below it.
MIN_TOP_COVERAGE = 0.95
# Material for every piece's side walls and bottom. None: the material on
# the model's bottom face (normally the base).
BASE_MATERIAL_NAME = None
# Gap between neighboring pieces, like the width of the cutting line. Each
# piece gives up half of it along every edge it shares with a neighbor;
# the puzzle's outer edges aren't moved.
PIECE_TOLERANCE_MM = 0.07
# Taller pieces need a bigger gap. On top of PIECE_TOLERANCE_MM, a piece's
# gap grows by EXTRA_TOLERANCE_MM for every EXTRA_TOLERANCE_PER_HEIGHT_MM of
# that piece's own height, measured from the model's bottom to the highest
# point of the model inside the piece's outline (scaled linearly, so a 50 mm
# tall piece adds 0.025 mm with the defaults). Set EXTRA_TOLERANCE_MM to 0
# to turn this off.
EXTRA_TOLERANCE_MM = 0.01
EXTRA_TOLERANCE_PER_HEIGHT_MM = 20.0
"""
Batch jigsaw cutter for Blender -- v7 (handles multi-object hierarchies)

Cuts a colored model into N interlocking jigsaw pieces and exports each
as its own .3mf. Handles models whose 3MF import created a parent Empty
with the colored geometry split across multiple child mesh objects
(common for multi-material 3MF imports, one mesh per color group).

HOW IT WORKS
------------
1. Every SVG is read directly (not through Blender's SVG importer, which
   shifts each piece by its file's cropped viewBox). The piece paths are
   drawn in full-puzzle coordinates, so all outlines together give the
   puzzle's extent. That extent is centered over the model's XY bounds
   with one uniform scale (SVG units -> model units), flipping SVG's
   downward Y to Blender's upward Y. Nothing is stretched.
2. Each color mesh is repaired on a copy (the originals are left untouched):
   duplicate vertices are merged, small holes filled, and normals made
   consistent. Common 3MF imports leave duplicate vertices along color
   boundaries, which shows up as open edges. Each outline becomes a closed
   prism spanning the model's full height; border edges are pushed just
   past the model so they never coincide with the model's own side walls,
   and shared edges are pulled in by half of the piece gap so neighbors end
   up that far apart. Each piece's gap is PIECE_TOLERANCE_MM plus an extra
   EXTRA_TOLERANCE_MM per EXTRA_TOLERANCE_PER_HEIGHT_MM of that piece's own
   height, found by casting rays down onto the model inside its outline.
3. Each color mesh is intersected with that prism. Afterwards:
   - Cut walls are found geometrically (vertical faces lying on the
     outline). Color meshes that aren't closed shells can make the boolean
     keep walls outside the mesh; a wall is kept only if its source mesh
     is both above and below it.
   - Cut walls, downward-facing faces, and the model's outer side walls
     all get one base material, so every piece has the same side/bottom
     color no matter which color regions it contains. Upward-facing
     surfaces keep their original colors.
   - If the cut leaves open edges where the source mesh was closed, the
     same repair is run on the cut piece and the edges are counted again.
   - The result is checked (within the model's height, within the outline,
     top surface covers the footprint, no open edges where the source mesh
     was closed) and retried without hole-tolerant mode if needed.
4. Every piece is exported with the whole model's filament palette, in the
   same order, even colors the piece doesn't use. The 3MF add-on normally
   numbers filaments by only the colors in each file, so importing one
   piece into another in Bambu Studio / Orca would shuffle colors.
5. Output files are named from the SVG's rN_cN (1-based), else the SVG's
   own name. The piece's object name inside the 3MF (what Bambu Studio /
   Orca show as the model name) is set to the same name, without .3mf.

HOW TO USE
----------
1. Export puzzle pieces as individual SVGs into one folder.
2. Import the colored .3mf (File > Import > 3D Manufacturing Format) and
   note the TOP-LEVEL object's name in the Outliner.
3. Edit CONFIG below and run in the Scripting workspace (Alt+P).
4. Test with TEST_LIMIT or ONLY_FILES first; watch the System Console.
"""

import bpy
import bmesh
import os
import re
import importlib
import addon_utils
import xml.etree.ElementTree as ElementTree
import numpy as np
from mathutils import Vector
from mathutils.bvhtree import BVHTree



FORCED_MODULE_NAME = "bl_ext.blender_org.ThreeMF_io"
# ------------------------------------------------------


def _get_export_3mf_func(module_name):
    """
    Tries module_name.api.export_3mf, then module_name.export_3mf.
    Requires a callable: the package may also have an export_3mf
    submodule that would otherwise be mistaken for the function.
    """
    for suffix in (".api", ""):
        try:
            mod = importlib.import_module(module_name + suffix)
        except ImportError:
            continue
        candidate = getattr(mod, "export_3mf", None)
        if callable(candidate):
            return candidate
    return None


def _find_3mf_api():
    candidates = []
    for mod in addon_utils.modules():
        bl_info = getattr(mod, "bl_info", {}) or {}
        display_name = str(bl_info.get("name", ""))
        module_name = mod.__name__
        if "3mf" in display_name.lower() or "3mf" in module_name.lower():
            candidates.append(module_name)

    for name in candidates:
        func = _get_export_3mf_func(name)
        if func is not None:
            return func

    all_modules = [
        (mod.__name__, str(getattr(mod, "bl_info", {}).get("name", "")))
        for mod in addon_utils.modules()
    ]
    raise RuntimeError(
        "Could not find a callable export_3mf function.\n"
        f"3mf-like candidates checked: {candidates or '(none)'}\n"
        "All installed add-ons (module name | display name):\n" +
        "\n".join(f"  {n} | {d}" for n, d in all_modules) +
        "\nSet FORCED_MODULE_NAME to the correct one and rerun."
    )


if FORCED_MODULE_NAME:
    export_3mf = _get_export_3mf_func(FORCED_MODULE_NAME)
    if export_3mf is None:
        raise RuntimeError(
            f"FORCED_MODULE_NAME '{FORCED_MODULE_NAME}' didn't yield a "
            "callable export_3mf (checked both '.api' and the bare module)."
        )
else:
    export_3mf = _find_3mf_api()


# ---------------- SVG outlines (plain Python, no bpy) ----------------
# Points are (x, y) tuples in SVG units.

def load_orca_color_hooks():
    """
    Returns (orca exporter module, its collect_face_colors, material_to_hex_color)
    from the 3MF add-on. Raises RuntimeError if this add-on version is laid
    out differently.
    """
    module = export_3mf.__module__
    for root in (module, module.rsplit(".", 1)[0]):
        try:
            orca = importlib.import_module(root + ".export_3mf.orca")
        except ImportError:
            continue
        collect = getattr(orca, "collect_face_colors", None)
        to_hex = getattr(orca, "material_to_hex_color", None)
        if callable(collect) and callable(to_hex):
            return orca, collect, to_hex
    raise RuntimeError(
        "Couldn't find the 3MF add-on's export_3mf.orca.collect_face_colors, which "
        "this script replaces during export to give every piece the same filament "
        "palette. This add-on version may be structured differently.")


def export_with_shared_palette(path, objects):
    """
    Exports with every file using filament_palette. Swaps the add-on's
    per-file color collection only for the duration of this call, so other
    exports in the Blender session are unaffected.
    """
    def shared_palette(blender_objects, *args, **kwargs):
        used = original_collect_face_colors(blender_objects, *args, **kwargs)
        missing = sorted(set(used) - set(filament_palette))
        if missing:
            print(f"  [WARN] {os.path.basename(path)}: colors {missing} aren't in the "
                  f"shared palette and won't be painted.")
        return dict(filament_palette)

    exporter_orca.collect_face_colors = shared_palette
    try:
        return export_3mf(path, objects=objects)
    finally:
        exporter_orca.collect_face_colors = original_collect_face_colors


def natural_sort_key(filename):
    """Sorts piece_2.svg before piece_10.svg instead of alphabetically."""
    return [int(chunk) if chunk.isdigit() else chunk
            for chunk in re.split(r'(\d+)', filename)]


def output_filename(svg_filename):
    """piece_rRR_cCC.3mf from a name containing rN_cN, else the SVG's own stem."""
    match = re.search(r'r(\d+)[_-]?c(\d+)', svg_filename, re.IGNORECASE)
    if match is None:
        return os.path.splitext(svg_filename)[0] + ".3mf"
    return f"piece_r{int(match.group(1)):02d}_c{int(match.group(2)):02d}.3mf"


def polygon_area(points):
    """Signed shoelace area; positive for counter-clockwise in a Y-up frame."""
    return 0.5 * sum(ax * by - bx * ay
                     for (ax, ay), (bx, by) in zip(points, points[1:] + points[:1]))


def polygon_bounds(points):
    """(min_x, min_y, max_x, max_y)"""
    xs = [p[0] for p in points]
    ys = [p[1] for p in points]
    return min(xs), min(ys), max(xs), max(ys)


def _segments_intersect(p1, p2, p3, p4):
    def orient(a, b, c):
        return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])

    def on_segment(a, b, c):
        return (min(a[0], b[0]) <= c[0] <= max(a[0], b[0])
                and min(a[1], b[1]) <= c[1] <= max(a[1], b[1]))

    d1, d2 = orient(p3, p4, p1), orient(p3, p4, p2)
    d3, d4 = orient(p1, p2, p3), orient(p1, p2, p4)
    if ((d1 > 0) != (d2 > 0)) and d1 != 0 and d2 != 0 \
            and ((d3 > 0) != (d4 > 0)) and d3 != 0 and d4 != 0:
        return True
    return ((d1 == 0 and on_segment(p3, p4, p1)) or (d2 == 0 and on_segment(p3, p4, p2))
            or (d3 == 0 and on_segment(p1, p2, p3)) or (d4 == 0 and on_segment(p1, p2, p4)))


def first_self_intersection(points):
    """Returns a point where non-adjacent edges touch or cross, or None."""
    count = len(points)
    edges = [(points[k], points[(k + 1) % count]) for k in range(count)]
    boxes = [(min(a[0], b[0]), min(a[1], b[1]), max(a[0], b[0]), max(a[1], b[1]))
             for a, b in edges]
    for k in range(count):
        for j in range(k + 2, count):
            if k == 0 and j == count - 1:
                continue  # these two edges share the closing vertex
            bk, bj = boxes[k], boxes[j]
            if bk[2] < bj[0] or bj[2] < bk[0] or bk[3] < bj[1] or bj[3] < bk[1]:
                continue
            if _segments_intersect(*edges[k], *edges[j]):
                return edges[k][0]
    return None


def extend_border_points(outline, puzzle_bounds, margin, tolerance):
    """
    Moves points lying on the puzzle's outer boundary outward by margin.
    Border pieces' outer cuts then land outside the model instead of exactly
    on its side walls, where coplanar faces confuse the boolean. Neighbors
    share their border corners, so they move identically and still tile.
    Returns (points, sides), where sides[k] is the set of boundary sides
    ("min_x", "max_x", "min_y", "max_y") point k lies on.
    """
    min_x, min_y, max_x, max_y = puzzle_bounds
    extended, sides = [], []
    for x, y in outline:
        on = set()
        if abs(x - min_x) <= tolerance:
            x -= margin
            on.add("min_x")
        elif abs(x - max_x) <= tolerance:
            x += margin
            on.add("max_x")
        if abs(y - min_y) <= tolerance:
            y -= margin
            on.add("min_y")
        elif abs(y - max_y) <= tolerance:
            y += margin
            on.add("max_y")
        extended.append((x, y))
        sides.append(on)
    return extended, sides


def inset_polygon(points, edge_insets):
    """
    Moves edge k (points[k] -> points[k+1]) inward by edge_insets[k], joining
    neighboring edges where their moved lines meet. Insets should be small
    relative to the outline's curvature; the caller checks the result for
    self-intersection.
    """
    count = len(points)
    orientation = 1.0 if polygon_area(points) > 0 else -1.0
    lines = []  # (point on moved line, direction, inward unit normal, inset)
    for k in range(count):
        (ax, ay), (bx, by) = points[k], points[(k + 1) % count]
        dx, dy = bx - ax, by - ay
        length = (dx * dx + dy * dy) ** 0.5
        nx, ny = -dy / length * orientation, dx / length * orientation
        inset = edge_insets[k]
        lines.append(((ax + nx * inset, ay + ny * inset), (dx, dy), (nx, ny), inset))

    result = []
    for k in range(count):
        (p1, d1, n1, inset1), (p2, d2, n2, inset2) = lines[k - 1], lines[k]
        px, py = points[k]
        cross = d1[0] * d2[1] - d1[1] * d2[0]
        largest = max(inset1, inset2)
        joined = None
        if abs(cross) > 1e-12 * (d1[0] ** 2 + d1[1] ** 2 + d2[0] ** 2 + d2[1] ** 2):
            t = ((p2[0] - p1[0]) * d2[1] - (p2[1] - p1[1]) * d2[0]) / cross
            joined = (p1[0] + d1[0] * t, p1[1] + d1[1] * t)
            # Nearly-reversing edges put the meeting point far away; cap it.
            if (joined[0] - px) ** 2 + (joined[1] - py) ** 2 > (4 * largest) ** 2:
                joined = None
        if joined is None:
            nx, ny = n1[0] + n2[0], n1[1] + n2[1]
            length = (nx * nx + ny * ny) ** 0.5 or 1.0
            joined = (px + nx / length * largest, py + ny / length * largest)
        result.append(joined)
    return result


def make_outline_locator(outline, tolerance):
    """Returns on_outline(x, y): whether a point is within tolerance of the outline's edges."""
    count = len(outline)
    segments = [(outline[k], outline[(k + 1) % count]) for k in range(count)]
    perimeter = sum(((b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2) ** 0.5 for a, b in segments)
    cell = max(perimeter / count, tolerance * 4)

    grid = {}
    for segment in segments:
        (ax, ay), (bx, by) = segment
        for gx in range(int((min(ax, bx) - tolerance) // cell), int((max(ax, bx) + tolerance) // cell) + 1):
            for gy in range(int((min(ay, by) - tolerance) // cell), int((max(ay, by) + tolerance) // cell) + 1):
                grid.setdefault((gx, gy), []).append(segment)

    def on_outline(x, y):
        for (ax, ay), (bx, by) in grid.get((int(x // cell), int(y // cell)), ()):
            dx, dy = bx - ax, by - ay
            length_sq = dx * dx + dy * dy
            t = 0.0 if length_sq == 0 else max(0.0, min(1.0, ((x - ax) * dx + (y - ay) * dy) / length_sq))
            if (x - ax - t * dx) ** 2 + (y - ay - t * dy) ** 2 <= tolerance * tolerance:
                return True
        return False

    return on_outline


_PATH_TOKEN = re.compile(r'[A-Za-z]|[-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?')


def parse_path_polygons(d):
    """
    Closed polygons from an SVG path made of straight segments (M, L, H, V,
    Z, absolute or relative). Raises ValueError for curves or arcs.
    """
    tokens = _PATH_TOKEN.findall(d)
    unsupported = sorted({t for t in tokens if t.isalpha() and t not in "MmLlHhVvZz"})
    if unsupported:
        raise ValueError(f"path uses unsupported commands {unsupported} "
                         f"(only straight-line paths are supported)")

    polygons, current = [], []
    x = y = start_x = start_y = 0.0
    command = None
    i = 0

    def number():
        nonlocal i
        if i >= len(tokens) or tokens[i].isalpha():
            raise ValueError(f"path is missing a number after '{command}'")
        value = float(tokens[i])
        i += 1
        return value

    while i < len(tokens):
        if tokens[i].isalpha():
            command = tokens[i]
            i += 1
            if command in "Zz":
                if current:
                    polygons.append(current)
                current = []
                x, y = start_x, start_y
                continue
        elif command is None:
            raise ValueError("path doesn't start with a command")

        relative = command.islower()
        if command in "Mm":
            if current:
                polygons.append(current)
            dx, dy = number(), number()
            x, y = (x + dx, y + dy) if relative else (dx, dy)
            start_x, start_y = x, y
            current = [(x, y)]
            command = "l" if relative else "L"  # extra coordinate pairs are line-tos
            continue
        if command in "Ll":
            dx, dy = number(), number()
            x, y = (x + dx, y + dy) if relative else (dx, dy)
        elif command in "Hh":
            dx = number()
            x = x + dx if relative else dx
        elif command in "Vv":
            dy = number()
            y = y + dy if relative else dy
        if not current:
            current = [(start_x, start_y)]
        current.append((x, y))
    if current:
        polygons.append(current)

    cleaned_polygons = []
    for polygon in polygons:
        cleaned = []
        for point in polygon:
            if not cleaned or point != cleaned[-1]:
                cleaned.append(point)
        while len(cleaned) > 1 and cleaned[0] == cleaned[-1]:
            cleaned.pop()
        if len(cleaned) >= 3:
            cleaned_polygons.append(cleaned)
    return cleaned_polygons


def load_svg_outline(svg_path):
    """
    Largest closed polygon among the SVG's <path> elements, in the path's
    own coordinates (the SVG's viewBox cropping is ignored). Raises
    ValueError if there is none or the file uses transforms.
    """
    root = ElementTree.parse(svg_path).getroot()
    polygons = []
    for element in root.iter():
        if "transform" in element.attrib:
            raise ValueError("SVG uses transform attributes, which aren't supported")
        if element.tag.rsplit("}", 1)[-1] == "path" and element.get("d"):
            polygons.extend(parse_path_polygons(element.get("d")))
    if not polygons:
        raise ValueError("no closed path found")
    return max(polygons, key=lambda p: abs(polygon_area(p)))


# ---------------- model analysis ----------------

UP = Vector((0.0, 0.0, 1.0))
DOWN = Vector((0.0, 0.0, -1.0))
# |normal.z| below this counts as a vertical face. Loose enough for float32
# rounding on small wall slivers; the on-outline test is the strict part.
VERTICAL_NZ = 1e-2


def collect_hierarchy(obj):
    """Returns [obj] + every descendant, recursively, depth-first."""
    result = [obj]
    for child in obj.children:
        result.extend(collect_hierarchy(child))
    return result


def world_vertices(mesh, matrix_world):
    """(N, 3) float64 array of the mesh's vertices in world space."""
    count = len(mesh.vertices)
    local = np.empty(count * 3, dtype=np.float32)
    mesh.vertices.foreach_get("co", local)
    matrix = np.array(matrix_world, dtype=np.float64)
    return local.reshape(count, 3) @ matrix[:3, :3].T + matrix[:3, 3]


def world_vertex_bounds(objects):
    """
    Axis-aligned world bounding box (min, max) of the objects' mesh
    vertices. Uses the vertices rather than bound_box, whose transformed
    corners overstate the extent of rotated objects.
    """
    all_verts = [world_vertices(obj.data, obj.matrix_world)
                 for obj in objects if len(obj.data.vertices)]
    return (Vector(np.min([v.min(axis=0) for v in all_verts], axis=0)),
            Vector(np.max([v.max(axis=0) for v in all_verts], axis=0)))


def world_bvh(mesh, matrix_world):
    """BVH tree of the mesh in world space, for inside/outside ray tests."""
    verts = [Vector(v) for v in world_vertices(mesh, matrix_world)]
    polygons = [tuple(p.vertices) for p in mesh.polygons]
    return BVHTree.FromPolygons(verts, polygons)


def _open_edges(bm):
    return sum(1 for e in bm.edges if not e.is_manifold)


def repair_bmesh(bm, merge_distance):
    """
    Merges duplicate vertices, fills holes of up to 64 edges, and, if that
    leaves no open edges, recalculates normals (on an open surface that can
    flip it the wrong way). Returns (new_faces, open_edges_after), where
    new_faces are the faces created by filling holes.
    """
    bmesh.ops.remove_doubles(bm, verts=list(bm.verts), dist=merge_distance)
    existing = set(bm.faces)
    boundary = [e for e in bm.edges if e.is_boundary]
    if boundary:
        filled = bmesh.ops.holes_fill(bm, edges=boundary, sides=64)["faces"]
        if filled:
            bmesh.ops.triangulate(bm, faces=filled)
    new_faces = [f for f in bm.faces if f not in existing]
    after = _open_edges(bm)
    if after == 0:
        bmesh.ops.recalc_face_normals(bm, faces=list(bm.faces))
    return new_faces, after


def prepare_source_mesh(obj, merge_distance_world):
    """
    Returns (repaired copy of the object's mesh, open_edges_before,
    open_edges_after). See repair_bmesh for what the repair does.
    """
    bm = bmesh.new()
    try:
        bm.from_mesh(obj.data)
        before = _open_edges(bm)

        scale = max(abs(s) for s in obj.matrix_world.to_scale()) or 1.0
        _, after = repair_bmesh(bm, merge_distance_world / scale)

        mesh = obj.data.copy()
        mesh.name = f"_repaired_{obj.name}"
        bm.to_mesh(mesh)
        return mesh, before, after
    finally:
        bm.free()


def bottom_material(objects, z_min, z_tolerance):
    """The material covering the most area on downward faces at the model's lowest Z."""
    totals = {}
    for obj in objects:
        mesh = obj.data
        count = len(mesh.polygons)
        if count == 0:
            continue
        verts = world_vertices(mesh, obj.matrix_world)
        for poly in mesh.polygons:
            if poly.material_index >= len(obj.material_slots):
                continue
            material = obj.material_slots[poly.material_index].material
            if material is None:
                continue
            zs = verts[list(poly.vertices), 2]
            if zs.max() <= z_min + z_tolerance:
                totals[material.name] = totals.get(material.name, 0.0) + poly.area
    if not totals:
        return None
    return bpy.data.materials[max(totals, key=totals.get)]


# ---------------- cutting ----------------

def build_cutter_mesh(outline, z_min, z_max, name):
    """
    Closed prism over a world-space XY outline spanning z_min..z_max, with
    outward normals. Raises ValueError if it can't be a watertight solid.
    """
    crossing = first_self_intersection(outline)
    if crossing is not None:
        raise ValueError(f"outline touches or crosses itself near "
                         f"({crossing[0]:.4g}, {crossing[1]:.4g})")

    bm = bmesh.new()
    try:
        bottom = [bm.verts.new((x, y, z_min)) for x, y in outline]
        top = [bm.verts.new((x, y, z_max)) for x, y in outline]
        try:
            caps = [bm.faces.new(top), bm.faces.new(list(reversed(bottom)))]
            count = len(outline)
            for k in range(count):
                j = (k + 1) % count
                bm.faces.new((bottom[k], bottom[j], top[j], top[k]))
        except ValueError as error:
            raise ValueError(f"outline has repeated points ({error})")

        bmesh.ops.triangulate(bm, faces=caps, quad_method='BEAUTY', ngon_method='BEAUTY')
        # Outline winding varies (and the SVG->Blender Y flip mirrors it);
        # an inside-out cutter makes the intersect discard the model.
        bmesh.ops.recalc_face_normals(bm, faces=list(bm.faces))

        bad_edges = sum(1 for e in bm.edges if not e.is_manifold)
        if bad_edges:
            raise ValueError(f"{bad_edges} non-manifold edges")
        if bm.calc_volume(signed=True) <= 0:
            raise ValueError("cutter has no volume or is inside-out")

        mesh = bpy.data.meshes.new(name)
        bm.to_mesh(mesh)
        return mesh
    finally:
        bm.free()


def material_slot_index(obj, material):
    """Index of the object's slot holding material, adding a slot if needed."""
    for index, slot in enumerate(obj.material_slots):
        if slot.material == material:
            return index
    obj.data.materials.append(material)
    return len(obj.material_slots) - 1


def _newell_normal(points):
    """Area-weighted (unnormalized) polygon normal."""
    normal = Vector((0.0, 0.0, 0.0))
    for a, b in zip(points, points[1:] + points[:1]):
        normal.x += (a.y - b.y) * (a.z + b.z)
        normal.y += (a.z - b.z) * (a.x + b.x)
        normal.z += (a.x - b.x) * (a.y + b.y)
    return normal


def _enclosed_by(bvh, points):
    """
    Whether the source mesh lies both above and below the face these points
    belong to. Two sample points guard against a ray slipping between
    triangles.
    """
    centroid = sum(points, Vector((0.0, 0.0, 0.0))) / len(points)
    for sample in (centroid, (centroid + points[0]) / 2):
        if bvh.ray_cast(sample, UP)[0] is not None and bvh.ray_cast(sample, DOWN)[0] is not None:
            return True
    return False


def finish_cut_object(obj, source, context):
    """
    Post-processes one boolean result in place: deletes cut walls the
    boolean kept outside the source mesh, and gives cut walls, downward
    faces, and the model's outer side walls the base material. Returns a
    list of problem strings.
    """
    base_index = material_slot_index(obj, context["base_material"])
    on_outline = context["on_outline"]
    z_min, z_max, z_tolerance = context["z_min"], context["z_max"], context["z_tolerance"]
    model_min, model_max = context["model_min"], context["model_max"]
    xy_tolerance = context["xy_tolerance"]

    def on_model_side(points):
        return any(all(abs(p[axis] - bound[axis]) <= xy_tolerance for p in points)
                   for axis in (0, 1) for bound in (model_min, model_max))

    bm = bmesh.new()
    try:
        bm.from_mesh(obj.data)
        matrix = obj.matrix_world
        world = [matrix @ v.co for v in bm.verts]
        vertex_on_outline = {}

        leaked = []
        for face in bm.faces:
            indices = [v.index for v in face.verts]
            points = [world[i] for i in indices]
            if any(not (z_min - z_tolerance <= p.z <= z_max + z_tolerance) for p in points):
                leaked.append(face)  # only the cutter reaches past the model's height
                continue
            normal = _newell_normal(points)
            if normal.length == 0:
                continue
            nz = normal.z / normal.length

            if abs(nz) < VERTICAL_NZ:
                for i in indices:
                    if i not in vertex_on_outline:
                        vertex_on_outline[i] = on_outline(world[i].x, world[i].y)
                if all(vertex_on_outline[i] for i in indices):
                    if not _enclosed_by(source["bvh"], points):
                        leaked.append(face)
                        continue
                    face.material_index = base_index
                elif on_model_side(points):
                    face.material_index = base_index
            elif nz < -0.5:
                face.material_index = base_index

        if leaked:
            bmesh.ops.delete(bm, geom=leaked, context='FACES')
        open_edges = _open_edges(bm)

        # The source was closed, so the cut should be too. Try the same
        # repair the source meshes got before reporting anything.
        open_before_repair = open_edges
        if source["closed"] and open_edges:
            new_faces, open_edges = repair_bmesh(bm, context["merge_distance"])
            for face in new_faces:
                points = [matrix @ v.co for v in face.verts]
                normal = _newell_normal(points)
                nz = normal.z / normal.length if normal.length else 0.0
                if nz > 0.5:
                    # Upward faces keep the color of the surface around them.
                    neighbors = [f for e in face.edges for f in e.link_faces if f is not face]
                    if neighbors:
                        face.material_index = neighbors[0].material_index
                else:
                    face.material_index = base_index
            if open_edges == 0:
                print(f"    {source['object'].name}: repair closed {open_before_repair} "
                      f"open edge(s) in the cut")
        bm.to_mesh(obj.data)
    finally:
        bm.free()

    problems = []
    if source["closed"] and open_edges:
        problems.append(f"{source['object'].name}: {open_edges} open edge(s) in the cut, "
                        f"so walls or surfaces are missing "
                        f"(repair left {open_edges} of {open_before_repair})")
    return problems


def discard(objects, meshes):
    """Removes temporary objects, then any of the meshes left without users."""
    for obj in objects:
        try:
            bpy.data.objects.remove(obj, do_unlink=True)
        except ReferenceError:
            pass  # already consumed by a join
    for mesh in meshes:
        try:
            if mesh.users == 0:
                bpy.data.meshes.remove(mesh)
        except ReferenceError:
            pass


def cut_model(sources, cutter_obj, collection, hole_tolerant, context,
              created_objects, created_meshes):
    """
    Intersects each source mesh with the cutter, post-processes each result,
    and joins them into one object. Returns (merged object or None, problems).
    Temporaries are appended to created_objects/created_meshes as they're
    made so the caller can discard them even if this raises.
    """
    results = []
    problems = []
    for source in sources:
        # The copy shares the source's mesh data, which is never modified:
        # the modifier result is baked into a new mesh and swapped in.
        temp = source["object"].copy()
        temp.data = source["mesh"]
        collection.objects.link(temp)
        created_objects.append(temp)
        temp.parent = None
        temp.matrix_world = source["object"].matrix_world.copy()
        temp.hide_viewport = False
        temp.hide_select = False
        temp.hide_set(False)

        mod = temp.modifiers.new(name="PuzzleCut", type='BOOLEAN')
        mod.operation = 'INTERSECT'
        mod.object = cutter_obj
        mod.solver = 'EXACT'
        mod.use_hole_tolerant = hole_tolerant

        depsgraph = bpy.context.evaluated_depsgraph_get()
        baked = bpy.data.meshes.new_from_object(temp.evaluated_get(depsgraph))
        created_meshes.append(baked)
        temp.modifiers.clear()
        temp.data = baked
        if baked.polygons:
            problems.extend(finish_cut_object(temp, source, context))
        if baked.polygons:
            results.append(temp)

    if not results:
        return None, problems
    if len(results) > 1:
        # Join rather than Boolean Union: union gets ambiguous where color
        # regions share coincident faces and can drop whole layers.
        bpy.ops.object.select_all(action='DESELECT')
        for obj in results:
            obj.select_set(True)
        bpy.context.view_layer.objects.active = results[0]
        bpy.ops.object.join()
    return results[0], problems


def inspect_piece(obj, outline, cutter_outline, z_min, z_max, z_tolerance):
    """
    Checks a cut result against the piece outline (for top coverage) and the
    cutter outline (for bounds). Returns {"problems": [...], "bounds_ok": bool,
    "coverage": float}.
    """
    if obj is None or not obj.data.polygons:
        return {"problems": ["cut result is empty"], "bounds_ok": False, "coverage": 0.0}

    world = [obj.matrix_world @ v.co for v in obj.data.vertices]
    problems = []
    low = min(c.z for c in world)
    high = max(c.z for c in world)
    height_ok = low >= z_min - z_tolerance and high <= z_max + z_tolerance
    if not height_ok:
        problems.append(f"spans z {low:.4g}..{high:.4g} but the model spans "
                        f"{z_min:.4g}..{z_max:.4g} (cutter geometry leaked in)")

    min_x, min_y, max_x, max_y = polygon_bounds(cutter_outline)
    xy_tolerance = ((max_x - min_x) ** 2 + (max_y - min_y) ** 2) ** 0.5 * 1e-3
    footprint_ok = all(min_x - xy_tolerance <= c.x <= max_x + xy_tolerance
                       and min_y - xy_tolerance <= c.y <= max_y + xy_tolerance
                       for c in world)
    if not footprint_ok:
        problems.append("extends outside the piece outline (model geometry leaked in)")

    # Projected XY area of upward-facing polygons.
    up_area = 0.0
    for poly in obj.data.polygons:
        up_area += max(_newell_normal([world[i] for i in poly.vertices]).z / 2, 0.0)
    footprint_area = abs(polygon_area(outline))
    coverage = up_area / footprint_area if footprint_area > 0 else 0.0
    if coverage < MIN_TOP_COVERAGE:
        problems.append(f"top surface covers only {coverage:.0%} of the footprint")

    return {"problems": problems, "bounds_ok": height_ok and footprint_ok,
            "coverage": coverage}


def attempt_score(attempt):
    report = attempt["report"]
    return (report["bounds_ok"], len(report["problems"]) == 0, report["coverage"])


def point_in_polygon(x, y, points):
    """Even-odd test of whether (x, y) is inside the polygon."""
    inside = False
    for (ax, ay), (bx, by) in zip(points, points[1:] + points[:1]):
        if (ay > y) != (by > y) and x < ax + (y - ay) / (by - ay) * (bx - ax):
            inside = not inside
    return inside


def piece_top_z(outline, sources, ray_start_z, samples=20):
    """
    Highest point of the source meshes inside the outline, found by casting
    rays straight down from above at a grid of points within the outline
    plus the outline's own corners. Returns None if nothing was hit.
    """
    min_x, min_y, max_x, max_y = polygon_bounds(outline)
    points = list(outline)
    for ix in range(samples):
        x = min_x + (ix + 0.5) / samples * (max_x - min_x)
        for iy in range(samples):
            y = min_y + (iy + 0.5) / samples * (max_y - min_y)
            if point_in_polygon(x, y, outline):
                points.append((x, y))
    top = None
    for x, y in points:
        for source in sources:
            hit = source["bvh"].ray_cast(Vector((x, y, ray_start_z)), DOWN)[0]
            if hit is not None and (top is None or hit.z > top):
                top = hit.z
    return top


# ---------------- setup ----------------

os.makedirs(OUTPUT_DIR, exist_ok=True)

original = bpy.data.objects.get(ORIGINAL_NAME)
if original is None:
    raise RuntimeError(f"Couldn't find object named '{ORIGINAL_NAME}'. "
                       f"Check the Outliner for the exact name.")
model_meshes = [obj for obj in collect_hierarchy(original)
                if obj.type == 'MESH' and len(obj.data.polygons)]
if not model_meshes:
    raise RuntimeError(f"'{ORIGINAL_NAME}' has no mesh objects in its hierarchy.")

model_min, model_max = world_vertex_bounds(model_meshes)
model_height = model_max.z - model_min.z
z_tolerance = max(model_height * 1e-3, 1e-7)
print(f"Model bounds: x {model_min.x:.4g}..{model_max.x:.4g}, "
      f"y {model_min.y:.4g}..{model_max.y:.4g}, z {model_min.z:.4g}..{model_max.z:.4g}")

if BASE_MATERIAL_NAME is not None:
    base_material = bpy.data.materials.get(BASE_MATERIAL_NAME)
    if base_material is None:
        raise RuntimeError(f"No material named '{BASE_MATERIAL_NAME}'. Materials: "
                           f"{', '.join(sorted(bpy.data.materials.keys()))}")
else:
    base_material = bottom_material(model_meshes, model_min.z, z_tolerance)
    if base_material is None:
        raise RuntimeError("Couldn't find a material on the model's bottom face; "
                           "set BASE_MATERIAL_NAME.")
print(f"Base material for piece sides and bottoms: '{base_material.name}'")

# Shared filament palette, computed with the add-on's own color logic so
# the hex values match what it writes: sorted by hex, numbered from 1.
exporter_orca, original_collect_face_colors, material_to_hex_color = load_orca_color_hooks()
palette_colors = set(original_collect_face_colors(model_meshes, False, lambda *a, **k: None,
                                                  export_hidden=True, include_disabled=True))
base_color = material_to_hex_color(base_material)
if base_color:
    palette_colors.add(base_color)
filament_palette = {color: index + 1 for index, color in enumerate(sorted(palette_colors))}

material_names_by_color = {}
for obj in model_meshes:
    for slot in obj.material_slots:
        if slot.material is not None:
            material_names_by_color.setdefault(
                material_to_hex_color(slot.material), set()).add(slot.material.name)
print("Filament palette shared by every piece:")
for color, index in filament_palette.items():
    names = ", ".join(sorted(material_names_by_color.get(color, ())))
    print(f"  filament {index}: {color} ({names})")

print("Preparing model meshes (repairing copies)...")
merge_distance = max(model_max.x - model_min.x, model_max.y - model_min.y) * 1e-6
sources = []
for obj in model_meshes:
    mesh, open_before, open_after = prepare_source_mesh(obj, merge_distance)
    sources.append({"object": obj, "mesh": mesh, "bvh": world_bvh(mesh, obj.matrix_world),
                    "closed": open_after == 0})
    state = "closed" if open_after == 0 else f"{open_after} open edges"
    if open_before != open_after:
        state += f" (was {open_before})"
    print(f"  {obj.name}: {len(mesh.polygons)} faces, {state}")

all_svg_files = sorted(
    (f for f in os.listdir(SVG_FOLDER) if f.lower().endswith('.svg')),
    key=natural_sort_key
)

# Pass 1: read every outline to find the whole puzzle's extent.
print(f"Reading {len(all_svg_files)} SVG outline(s)...")
svg_outlines = {}
for filename in all_svg_files:
    try:
        svg_outlines[filename] = load_svg_outline(os.path.join(SVG_FOLDER, filename))
    except (ValueError, ElementTree.ParseError) as error:
        print(f"  [WARN] {filename}: {error}, skipping.")
if not svg_outlines:
    raise RuntimeError(f"No usable SVG outlines in {SVG_FOLDER}")

all_bounds = [polygon_bounds(o) for o in svg_outlines.values()]
puzzle_min_x = min(b[0] for b in all_bounds)
puzzle_min_y = min(b[1] for b in all_bounds)
puzzle_width = max(b[2] for b in all_bounds) - puzzle_min_x
puzzle_height = max(b[3] for b in all_bounds) - puzzle_min_y
puzzle_bounds = (puzzle_min_x, puzzle_min_y,
                 puzzle_min_x + puzzle_width, puzzle_min_y + puzzle_height)
border_margin = max(puzzle_width, puzzle_height) * 0.01
border_tolerance = max(puzzle_width, puzzle_height) * 1e-6

# Pieces that tile the puzzle exactly have areas summing to its rectangle.
area_ratio = (sum(abs(polygon_area(o)) for o in svg_outlines.values())
              / (puzzle_width * puzzle_height))
print(f"SVG puzzle: {puzzle_width:.4g} x {puzzle_height:.4g} SVG units; pieces cover "
      f"{area_ratio:.2%} of it (100% = no gaps or overlaps, with every piece present)")
if area_ratio > 1.005:
    print("  [WARN] the SVG outlines overlap each other; neighboring pieces "
          "will share a sliver of the model.")

model_width = model_max.x - model_min.x
model_depth = model_max.y - model_min.y
scale_x = model_width / puzzle_width
scale_y = model_depth / puzzle_height
# One scale for both axes. Averaging the two only matters if the sizes
# disagree, which is warned about below.
puzzle_scale = (scale_x + scale_y) / 2
model_center_x = (model_min.x + model_max.x) / 2
model_center_y = (model_min.y + model_max.y) / 2
puzzle_center_x = puzzle_min_x + puzzle_width / 2
puzzle_center_y = puzzle_min_y + puzzle_height / 2
print(f"Model XY size {model_width:.4g} x {model_depth:.4g}; scale {puzzle_scale:.6g} "
      f"model units per SVG unit")
if abs(scale_x / scale_y - 1) > 0.01:
    print(f"  [WARN] model and puzzle aspect ratios differ by {abs(scale_x / scale_y - 1):.1%} "
          f"(x scale {scale_x:.6g}, y scale {scale_y:.6g}). Pieces aren't stretched, so "
          f"the puzzle won't line up exactly with the model's edges.")


def to_model_xy(point):
    # SVG Y points down; Blender's Y points up (top of the puzzle = +Y).
    return (model_center_x + (point[0] - puzzle_center_x) * puzzle_scale,
            model_center_y - (point[1] - puzzle_center_y) * puzzle_scale)


output_names = {f: output_filename(f) for f in svg_outlines}
seen_outputs = {}
for filename, out_name in output_names.items():
    if out_name in seen_outputs:
        print(f"  [WARN] {filename} and {seen_outputs[out_name]} both export to "
              f"{out_name}; one will overwrite the other.")
    seen_outputs[out_name] = filename

files_to_cut = [f for f in all_svg_files if f in svg_outlines]
if ONLY_FILES is not None:
    files_to_cut = [f for f in files_to_cut if f in ONLY_FILES]
if TEST_LIMIT is not None:
    files_to_cut = files_to_cut[:TEST_LIMIT]

# The cutter extends past the model's top and bottom so the cut is clean.
z_padding = max(model_height * 0.1, 1e-4)
cutter_z_min = model_min.z - z_padding
cutter_z_max = model_max.z + z_padding
xy_tolerance = max(model_width, model_depth) * 1e-6

# Blender units -> meters is the scene's unit scale.
unit_scale = bpy.context.scene.unit_settings.scale_length
print(f"Gap between pieces: {PIECE_TOLERANCE_MM} mm, plus {EXTRA_TOLERANCE_MM} mm per "
      f"{EXTRA_TOLERANCE_PER_HEIGHT_MM} mm of each piece's own height "
      f"(model is {model_width * unit_scale * 1000:.1f} x {model_depth * unit_scale * 1000:.1f} x "
      f"{model_height * unit_scale * 1000:.1f} mm)")

# ---------------- pass 2: cut and export ----------------

print(f"Cutting {len(files_to_cut)} piece(s)...")
work_collection = bpy.context.collection
manifest_rows = []   # (sequential_index, svg_filename, out_filename, status)
suspect_pieces = []  # (svg_filename, problem description)

for i, filename in enumerate(files_to_cut):
    out_name = output_names[filename]
    out_path = os.path.join(OUTPUT_DIR, out_name)
    outline = [to_model_xy(p) for p in svg_outlines[filename]]
    extended, border_sides = extend_border_points(
        svg_outlines[filename], puzzle_bounds, border_margin, border_tolerance)
    cutter_outline = [to_model_xy(p) for p in extended]
    top_z = piece_top_z(outline, sources, model_max.z + z_padding)
    piece_height_mm = (max(top_z - model_min.z, 0.0) if top_z is not None
                       else model_height) * unit_scale * 1000
    extra_tolerance_mm = 0.0
    if EXTRA_TOLERANCE_PER_HEIGHT_MM > 0:
        extra_tolerance_mm = EXTRA_TOLERANCE_MM * piece_height_mm / EXTRA_TOLERANCE_PER_HEIGHT_MM
    gap_mm = PIECE_TOLERANCE_MM + extra_tolerance_mm
    piece_gap = gap_mm * 0.001 / unit_scale
    if piece_gap > 0:
        count = len(cutter_outline)
        # Edges along the puzzle's outer boundary keep their position.
        edge_insets = [0.0 if border_sides[k] & border_sides[(k + 1) % count] else piece_gap / 2
                       for k in range(count)]
        cutter_outline = inset_polygon(cutter_outline, edge_insets)

    try:
        cutter_mesh = build_cutter_mesh(cutter_outline, cutter_z_min, cutter_z_max,
                                        f"_cutter_data_{i}")
    except ValueError as error:
        print(f"  [WARN] {filename}: can't build cutter: {error}. Skipping; check it by hand.")
        suspect_pieces.append((filename, f"skipped, cutter: {error}"))
        continue

    cutter_obj = bpy.data.objects.new(f"_cutter_tmp_{i}", cutter_mesh)
    work_collection.objects.link(cutter_obj)
    cutter_obj.hide_render = True
    context = {
        "base_material": base_material,
        "on_outline": make_outline_locator(cutter_outline, xy_tolerance * 10),
        "z_min": model_min.z,
        "z_max": model_max.z,
        "z_tolerance": z_tolerance,
        "model_min": model_min,
        "model_max": model_max,
        "xy_tolerance": xy_tolerance * 10,
        "merge_distance": merge_distance,
    }

    best = None
    try:
        # Cut, check, and retry without hole-tolerant mode if the result
        # fails; keep only the best attempt.
        for hole_tolerant in (True, False):
            label = "EXACT hole-tolerant" if hole_tolerant else "EXACT"
            attempt = {"label": label, "objects": [], "meshes": []}
            try:
                attempt["merged"], cut_problems = cut_model(
                    sources, cutter_obj, work_collection, hole_tolerant, context,
                    attempt["objects"], attempt["meshes"])
            except Exception as error:
                print(f"  [RETRY] {filename}: {label} failed: {error}")
                discard(attempt["objects"], attempt["meshes"])
                continue
            attempt["report"] = inspect_piece(attempt["merged"], outline, cutter_outline,
                                              model_min.z, model_max.z, z_tolerance)
            attempt["report"]["problems"] = cut_problems + attempt["report"]["problems"]

            if best is None or attempt_score(attempt) > attempt_score(best):
                if best is not None:
                    discard(best["objects"], best["meshes"])
                best = attempt
            else:
                discard(attempt["objects"], attempt["meshes"])

            if not attempt["report"]["problems"]:
                break
            print(f"  [RETRY] {filename}: {label}: "
                  f"{'; '.join(attempt['report']['problems'])}")

        if best is None or best["merged"] is None:
            print(f"  [WARN] {filename}: every boolean attempt came out empty, not exported.")
            suspect_pieces.append((filename, "not exported, empty cut"))
            continue

        # Without this the exported model keeps the source object's name
        # (e.g. "Base.001"), which is what Bambu Studio / Orca list.
        piece_name = os.path.splitext(out_name)[0]
        best["merged"].name = piece_name
        best["merged"].data.name = piece_name
        if best["merged"].name != piece_name:
            print(f"  [WARN] {filename}: an object named '{piece_name}' already exists, "
                  f"so the model is named '{best['merged'].name}' instead.")

        result = export_with_shared_palette(out_path, [best["merged"]])
        problems = best["report"]["problems"]
        status = "ok" if not problems else "suspect"
        print(f"  [{i:03d}] {filename}: {result.status} -> {out_name} "
              f"({best['label']}, top coverage {best['report']['coverage']:.0%}, "
              f"height {piece_height_mm:.1f} mm, gap {gap_mm:.3f} mm)")
        if problems:
            print(f"  [WARN] {filename}: best attempt still has problems: "
                  f"{'; '.join(problems)}")
            suspect_pieces.append((filename, '; '.join(problems)))
        manifest_rows.append((i, filename, out_name, status))

    finally:
        if best is not None:
            discard(best["objects"], best["meshes"])
        discard([cutter_obj], [cutter_mesh])

for source in sources:
    if source["mesh"] is not source["object"].data and source["mesh"].users == 0:
        bpy.data.meshes.remove(source["mesh"])

print("Done cutting and exporting.")

if manifest_rows:
    manifest_path = os.path.join(OUTPUT_DIR, "piece_manifest.csv")
    with open(manifest_path, "w") as manifest:
        manifest.write("sequential_index,svg_filename,filename,status\n")
        for row in manifest_rows:
            manifest.write(",".join(str(value) for value in row) + "\n")
    print(f"Manifest written to {manifest_path}")

if suspect_pieces:
    print(f"{len(suspect_pieces)} piece(s) need a manual check:")
    for name, reason in suspect_pieces:
        print(f"  {name}: {reason}")

print("Done.")