bl_info = {
    "name": "Celtic Knot Shape",
    "author": "Robert J. Wallace",
    "version": (1, 0),
    "blender": (4, 0, 0),
    "location": "View3D > Add > Mesh > Celtic Knot Shape",
    "description": "Builds a design from the Celtic Knot Shapes page: knotwork filling any outline",
    "category": "Add Mesh",
}

# How to use: install this file as an add-on (Edit > Preferences > Add-ons >
# Install from Disk), or open it in the Text Editor and press Run Script.
# On the Celtic Knot Shapes page press "Copy for Blender", then in Blender
# choose Add > Mesh > Celtic Knot Shape. If the clipboard doesn't hold a
# design it asks for a file saved with "Download for Blender" instead.
#
# The page sends the outline and each strand's path (its control points and
# which crossings it goes over or under); this script draws the curves, the
# dips at crossings and the solids itself, at full detail, by the same rules
# as the page. Everything is in mm, flat on the XY plane, the top facing +Z.
# Each run replaces the Shape_* objects from the run before.
#
# Headless:
#   blender -b --factory-startup --python make_knot_shape.py -- design.json [out.stl]

import json
import math
import os
import sys
import bpy
import bmesh
import numpy as np
from mathutils import Vector
from mathutils.geometry import tessellate_polygon

# ---------------------------------------------------------------- settings (mm)
SAMPLES_PER_SPAN = 12      # curve points from one crossing (or turn) to the next; more is smoother
CARVE_RESOLUTION = 0.035   # Carved style: surface grid spacing (smaller is finer, and slower)
SIMPLIFY_ANGLE = 0.0       # degrees: above 0, faces this close to flat are merged (Decimate, Planar), so the
                           # flat tops, floor and base take far fewer triangles; about 1 suits printing
MAKE_PRINT_SOLID = False   # also make Shape_Print: everything fused into one solid (voxel remesh)
PRINT_VOXEL = 0.03         # voxel size for that solid (mm)
WEIGH_SAMPLES = 250000     # how many places the piece is measured at to weigh it
DESIGN_FILE = ''           # a design saved with "Download for Blender", used when the clipboard has none

ALLOYS = {'STERLING': ('sterling silver', 10.36), '10K': ('10k gold', 11.6), '14K': ('14k gold', 13.1),
          '18K': ('18k gold', 15.5), 'PT950': ('platinum 950', 20.1), 'PLA': ('PLA plastic', 1.24)}
METAL_COLORS = {'GOLD': (1.0, 0.72, 0.28), 'SILVER': (0.97, 0.96, 0.94),
                'BRONZE': (0.80, 0.45, 0.25), 'JADE': (0.25, 0.60, 0.48)}
HEADER = '# Exported from Celtic Knot Shapes'

# cross-sections across a swept strand: (across -1..1, up 0..1, sharp corner)
PROFILES = {
    'RIBBON': [(-1, 0, 1), (-1, 0.7, 1), (-0.75, 1, 1), (0.75, 1, 1), (1, 0.7, 1), (1, 0, 1)],
    'ROUNDED': [(-1, 0, 1)] + [(math.cos(math.pi * (1 - k / 12)), 0.25 + 0.75 * math.sin(math.pi * (1 - k / 12)), 0)
                               for k in range(13)] + [(1, 0, 1)],
}
LAST_REPORT = ['']


# DOUBLE: each strand is two rims (double_rim: each one's share of the
# strand's width) with the middle open down to the plate, or to a low floor
# when there is no plate. Where a strand goes under, its rims stop inside the
# strand on top, so that strand's open middle stays clear, as in fretwork.
RIM_PROFILE = [(-1, 0, 1), (-1, 0.9, 1), (-0.5, 1, 1), (0.5, 1, 1), (1, 0.9, 1), (1, 0, 1)]
FLOOR_PROFILE = [(-1, 0, 1), (-1, 1, 1), (1, 1, 1), (1, 0, 1)]
FLOOR = 0.25               # the floor's height, as a share of the strand's, with no plate


# ---------------------------------------------------------------- the design
def read_design(text):
    text = text.lstrip()
    if text.startswith(HEADER):
        text = text[len(HEADER):]
    d = json.loads(text)
    if d.get('format') != 'celtic-knot-shape':
        raise ValueError("That isn't a design from the Celtic Knot Shapes page")
    return d


def catmull_rom(ctrl, S):
    """Closed centripetal Catmull-Rom through the control points, S samples a
    span (as the page draws it)."""
    P = np.asarray(ctrl, float)
    n = len(P)
    out = np.empty((n * S, 2))
    s = np.arange(S) / S
    for k in range(n):
        p0, p1, p2, p3 = P[(k - 1) % n], P[k], P[(k + 1) % n], P[(k + 2) % n]
        d = lambda a, b: max(math.sqrt(math.hypot(*(b - a))), 1e-4)
        t0 = 0.0
        t1 = t0 + d(p0, p1)
        t2 = t1 + d(p1, p2)
        t3 = t2 + d(p2, p3)
        t = (t1 + (t2 - t1) * s)[:, None]
        A1 = (p0 * (t1 - t) + p1 * (t - t0)) / (t1 - t0)
        A2 = (p1 * (t2 - t) + p2 * (t - t1)) / (t2 - t1)
        A3 = (p2 * (t3 - t) + p3 * (t - t2)) / (t3 - t2)
        B1 = (A1 * (t2 - t) + A2 * (t - t0)) / (t2 - t0)
        B2 = (A2 * (t3 - t) + A3 * (t - t1)) / (t3 - t1)
        out[k * S:(k + 1) * S] = (B1 * (t2 - t) + B2 * (t - t1)) / (t2 - t1)
    return out


def polyline(ctrl, S):
    P = np.asarray(ctrl, float)
    Q = np.roll(P, -1, axis=0)
    s = (np.arange(S) / S)[None, :, None]
    return (P[:, None] * (1 - s) + Q[:, None] * s).reshape(-1, 2)


class Strand:
    def __init__(self, data, S, straight, settings):
        # its own look (panels can differ), or the design's
        self.style = data.get('style', settings['style'])
        self.width = data.get('width', settings['strand_width'])
        self.double_rim = data.get('double_rim', settings.get('double_rim', 0.2))
        self.height = data.get('height', settings['strand_height'])
        self.dip = data.get('dip', settings['dip'])
        # stands on a plate (its panel's), and so isn't woven
        self.on_plate = data.get('plate', settings['plate'])
        self.woven = settings.get('woven', False) and not self.on_plate
        self.ctrl = data['points']
        self.crossing = data['crossing']
        self.over = data['over']
        self.S = S
        # a spiral: its path itself, already sampled, open unless it closes
        self.open = data.get('open', False)
        if data.get('path'):
            self.pts = np.asarray(self.ctrl, float)
        else:
            self.pts = polyline(self.ctrl, S) if straight else catmull_rom(self.ctrl, S)
        # a tapered spiral: each point's share of the width (else all 1)
        self.scale = np.asarray(data['width_scale'], float) if data.get('width_scale') else np.ones(len(self.pts))
        seg = np.linalg.norm(np.roll(self.pts, -1, axis=0) - self.pts, axis=1)
        self.arc = np.concatenate([[0.0], np.cumsum(seg)])
        self.length = self.arc[-1]

    def along(self, j):
        """distance along the strand from sample j to every sample"""
        d = np.abs(self.arc[:-1] - self.arc[j])
        return np.minimum(d, self.length - d)

    def tangent(self, j):
        n = len(self.pts)
        t = self.pts[(j + 1) % n] - self.pts[(j - 1) % n]
        return t / (np.linalg.norm(t) or 1)


def crossings_of(strands):
    """[(over strand, sample, under strand, sample, sin of the crossing angle)]"""
    by_edge = {}
    for si, s in enumerate(strands):
        for k, e in enumerate(s.crossing):
            if e >= 0:
                by_edge.setdefault(e, {})['over' if s.over[k] else 'under'] = (si, k * s.S)
    out = []
    for e, p in by_edge.items():
        if 'over' in p and 'under' in p:
            (oi, oj), (ui, uj) = p['over'], p['under']
            a, b = strands[oi].tangent(oj), strands[ui].tangent(uj)
            out.append((oi, oj, ui, uj, abs(a[0] * b[1] - a[1] * b[0])))
    return out


# ---------------------------------------------------------------- meshes
def mesh_object(name, verts, faces, material, smooth_angle=None):
    """A mesh object from an (N, 3) array of points and a list of faces."""
    mesh = bpy.data.meshes.new(name)
    mesh.from_pydata(np.asarray(verts, float).tolist(), [], faces if isinstance(faces, list) else faces.tolist())
    mesh.validate(clean_customdata=False)
    mesh.materials.append(material)
    if smooth_angle is not None:
        for p in mesh.polygons:
            p.use_smooth = True
        if hasattr(mesh, "set_sharp_from_angle"):
            mesh.set_sharp_from_angle(angle=smooth_angle)
    obj = bpy.data.objects.new(name, mesh)
    bpy.context.scene.collection.objects.link(obj)
    return obj


def signed_volume(verts, tris):
    v = np.asarray(verts, float)
    t = np.asarray(tris)
    a, b, c = v[t[:, 0]], v[t[:, 1]], v[t[:, 2]]
    return float(np.einsum('ij,ij->i', a, np.cross(b, c)).sum() / 6)


def sweep(pts, top, bottom, width, profile, closed=True):
    """A strand swept along its samples; the profile is scaled to the strand's
    height at each sample, so it dips under and rises over. An open one gets
    flat ends."""
    n = len(pts)
    cols, col_of = [], []
    for k, p in enumerate(profile):
        col_of.append((len(cols), len(cols) + (1 if p[2] else 0)))
        cols += [k, k] if p[2] else [k]
    m, C = len(profile), len(cols)
    seg_cols = [(col_of[k][1], col_of[(k + 1) % m][0]) for k in range(m)]
    t = np.roll(pts, -1, axis=0) - np.roll(pts, 1, axis=0)
    if not closed:
        t[0], t[-1] = pts[1] - pts[0], pts[-1] - pts[-2]
    t /= np.maximum(np.linalg.norm(t, axis=1, keepdims=True), 1e-12)
    nrm = np.stack([-t[:, 1], t[:, 0]], axis=1)
    prof = np.array([profile[c][:2] for c in cols], float)
    # (width: one for all, or each sample's)
    wd = np.asarray(width, float).reshape(-1, 1, 1) if np.ndim(width) else width
    xy = pts[:, None, :] + nrm[:, None, :] * (prof[None, :, 0:1] * wd / 2)
    bottom = np.broadcast_to(np.asarray(bottom, float), top.shape)[:, None]
    z = bottom + prof[None, :, 1] * (top[:, None] - bottom)
    verts = np.concatenate([xy, z[..., None]], axis=2).reshape(-1, 3)
    i = np.arange(n if closed else n - 1)[:, None]
    i2 = (i + 1) % n
    tris = []
    for c0, c1 in seg_cols:
        a, b, d, e = i * C + c0, i * C + c1, i2 * C + c0, i2 * C + c1
        tris.append(np.concatenate([a, d, b], 1))
        tris.append(np.concatenate([b, d, e], 1))
    if not closed:
        # the ends: a fan over the profile, wound against the sides
        e = (n - 1) * C
        for k in range(1, m - 1):
            p0, p1, p2 = col_of[0][0], col_of[k][0], col_of[k + 1][0]
            tris.append(np.array([[p0, p1, p2], [e + p0, e + p2, e + p1]]))
    tris = np.concatenate(tris)
    if signed_volume(verts, tris) < 0:
        tris = tris[:, [0, 2, 1]]
    return verts, tris


def extrude(groups, z0, z1):
    """Outlines (outer loops with their holes) as solids from z0 up to z1."""
    verts, faces = [], []
    for g in groups:
        loops = [g['outer']] + g['holes']
        flat = [p for l in loops for p in l]
        base = len(verts)
        n = len(flat)
        verts += [(x, y, z1) for x, y in flat] + [(x, y, z0) for x, y in flat]
        for tri in tessellate_polygon([[Vector((x, y, 0)) for x, y in l] for l in loops]):
            a, b, c = (flat[i] for i in tri)
            up = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]) > 0
            tri = tri if up else (tri[0], tri[2], tri[1])
            faces.append([base + i for i in tri])
            faces.append([base + n + i for i in (tri[0], tri[2], tri[1])])
        start = 0
        for l in loops:              # outer counter-clockwise, holes clockwise: the solid on the left
            m = len(l)
            for k in range(m):
                a, b = base + start + k, base + start + (k + 1) % m
                faces.append([a, a + n, b + n, b])
            start += m
    return verts, faces


# ---------------------------------------------------------------- swept styles
def swept_heights(strand, crossings, si, height, dip, cell, base, woven):
    """(top, bottom) at each sample. Woven: the strand keeps its thickness and
    moves up over and down under, so the back shows the over and under too."""
    R = 0.42 * cell
    top = np.full(len(strand.pts), height)
    shift = np.zeros(len(strand.pts))
    lift = dip * height / 2
    for oi, oj, ui, uj, sn in crossings:
        for idx, j, over in ((oi, oj, True), (ui, uj, False)):
            if idx != si:
                continue
            # only the samples near this crossing: elsewhere b is 0, and the
            # plain height would undo the other crossings' dips and rises
            d = strand.along(j)
            near = d < R
            b = 0.5 * (1 + np.cos(np.pi * d[near] / R))
            if woven:
                shift[near] = np.maximum(shift[near], lift * b) if over else np.minimum(shift[near], -lift * b)
            else:
                top[near] = np.maximum(top[near], height * (1 + 0.12 * b)) if over else np.minimum(top[near], height * (1 - dip * b))
    if woven:
        top = height + shift
    return top, base + shift


def offset_line(pts, d, open_=False):
    """a line (closed unless open_) moved sideways by d (left of its direction when d > 0)"""
    t = np.roll(pts, -1, axis=0) - np.roll(pts, 1, axis=0)
    if open_:
        t[0], t[-1] = pts[1] - pts[0], pts[-1] - pts[-2]
    t /= np.maximum(np.linalg.norm(t, axis=1, keepdims=True), 1e-12)
    d = np.asarray(d, float).reshape(-1, 1) if np.ndim(d) else d    # (or each point by its own)
    return pts + d * np.stack([-t[:, 1], t[:, 0]], axis=1)


def around(strand, j, r):
    """the samples of a closed strand within distance r (along it) of sample j, in order"""
    n = len(strand.pts)
    k = np.arange(j - n // 2, j + n // 2 + 1)
    k = k[strand.along(j)[k % n] <= r]
    return k % n


def to_line(pts, q):
    """distance from each point q to the line through pts"""
    a, b = pts[:-1], pts[1:]
    ab = b - a
    t = np.clip(((q[:, None, :] - a[None]) * ab[None]).sum(2) / np.maximum((ab * ab).sum(1), 1e-12)[None], 0, 1)
    return np.linalg.norm(a[None] + t[..., None] * ab[None] - q[:, None, :], axis=2).min(1)


def pieces(q, f, top, bot):
    """the runs of a closed line where f > 0, each ended where f crosses 0"""
    n = len(f)
    cut_at = np.nonzero(f <= 0)[0]
    if not len(cut_at):
        return [(q, top, bot, True)]
    s = cut_at[0]
    v = np.concatenate([q, top[:, None], bot[:, None]], axis=1)

    def cut(i, j):
        t = f[i] / (f[i] - f[j]) if np.isfinite(f[i]) else 0.0
        return v[i] + t * (v[j] - v[i])

    out, run = [], None
    for k in range(1, n + 1):
        i, prev = (s + k) % n, (s + k - 1) % n
        if f[i] > 0:
            if run is None:
                run = [cut(i, prev)]
            run.append(v[i])
        elif run is not None:
            run.append(cut(prev, i))
            r = np.array(run)
            out.append((r[:, :2], r[:, 2], r[:, 3], False))
            run = None
    return out


def double_strands(strands, crossings, heights, s):
    """the DOUBLE style's pieces: each strand's two rims, cut where they pass
    under another strand (unless woven, when they pass beneath), and the floor
    between them when there's no plate. heights: each strand's (top, bottom)"""
    rw = [min(st.double_rim, 0.5) * st.width for st in strands]
    inner = [st.width / 2 - 0.6 * r for st, r in zip(strands, rw)]
    # (a tapered strand's rims come in with it, each keeping its share of the width)
    rims = [[[offset_line(st.pts, side * (st.width / 2 - r / 2) * st.scale, st.open), np.full(len(st.pts), np.inf)] for side in (-1, 1)]
            for st, r in zip(strands, rw)]
    for oi, oj, ui, uj, sn in crossings:
        if strands[ui].woven:
            continue
        # no further from the crossing than its cut can reach, so that other
        # places where the two strands meet aren't cut too
        reach = strands[oi].width / max(sn, 0.25)
        over = strands[oi].pts[around(strands[oi], oj, reach)]
        near = around(strands[ui], uj, reach)
        for q, f in rims[ui]:
            f[near] = np.minimum(f[near], to_line(over, q[near]) - inner[oi])
    out = []
    for si, st in enumerate(strands):
        top, bot = heights[si]
        for q, f in rims[si]:
            if np.isinf(f).all():
                # never cut (a spiral): whole, tapered with the strand
                out.append(sweep(q, top, bot, rw[si] * st.scale, RIM_PROFILE, not st.open))
                continue
            for pts, ptop, pbot, closed in pieces(q, f, top, bot):
                out.append(sweep(pts, ptop, pbot, rw[si], RIM_PROFILE, closed and not st.open))
        if not st.on_plate:
            out.append(sweep(st.pts, bot + FLOOR * (top - bot), bot + 0.01, (st.width - rw[si]) * st.scale, FLOOR_PROFILE, not st.open))
    return out


# ---------------------------------------------------------------- carved style
class Outline:
    """Signed distance to the outline (positive inside), sampled on a grid and
    read back with its slope anywhere."""
    def __init__(self, groups, res, reach):
        loops = [np.asarray(l, float) for g in groups for l in [g['outer']] + g['holes']]
        allp = np.concatenate(loops)
        self.x0, self.y0 = allp.min(0) - reach - 2 * res
        x1, y1 = allp.max(0) + reach + 2 * res
        self.res = res
        self.nx = int(math.ceil((x1 - self.x0) / res)) + 1
        self.ny = int(math.ceil((y1 - self.y0) / res)) + 1
        xs = self.x0 + np.arange(self.nx) * res
        ys = self.y0 + np.arange(self.ny) * res
        A = np.concatenate([l for l in loops])
        B = np.concatenate([np.roll(l, -1, axis=0) for l in loops])
        # inside: count the outline crossings to the left of each point, row by row
        inside = np.zeros((self.ny, self.nx), bool)
        for j, y in enumerate(ys):
            hit = (A[:, 1] > y) != (B[:, 1] > y)
            if not hit.any():
                continue
            a, b = A[hit], B[hit]
            xc = np.sort(a[:, 0] + (y - a[:, 1]) * (b[:, 0] - a[:, 0]) / (b[:, 1] - a[:, 1]))
            inside[j] = np.searchsorted(xc, xs) % 2 == 1
        # distance, worked out only near the outline (further in it is "far")
        dist = np.full((self.ny, self.nx), reach + res)
        for a, b in zip(A, B):
            lo, hi = np.minimum(a, b) - reach, np.maximum(a, b) + reach
            i0, i1 = max(0, int((lo[0] - self.x0) / res)), min(self.nx - 1, int((hi[0] - self.x0) / res) + 1)
            j0, j1 = max(0, int((lo[1] - self.y0) / res)), min(self.ny - 1, int((hi[1] - self.y0) / res) + 1)
            if i0 > i1 or j0 > j1:
                continue
            X, Y = np.meshgrid(xs[i0:i1 + 1], ys[j0:j1 + 1])
            e = b - a
            t = np.clip(((X - a[0]) * e[0] + (Y - a[1]) * e[1]) / max(e @ e, 1e-12), 0, 1)
            d = np.hypot(X - a[0] - t * e[0], Y - a[1] - t * e[1])
            blk = dist[j0:j1 + 1, i0:i1 + 1]
            np.minimum(blk, d, out=blk)
        self.sdf = np.where(inside, dist, -dist)

    def sample(self, X, Y):
        fx = np.clip((X - self.x0) / self.res, 0, self.nx - 1.001)
        fy = np.clip((Y - self.y0) / self.res, 0, self.ny - 1.001)
        i, j = fx.astype(int), fy.astype(int)
        tx, ty = fx - i, fy - j
        s = self.sdf
        a, b, c, d = s[j, i], s[j, i + 1], s[j + 1, i], s[j + 1, i + 1]
        v = (a * (1 - tx) + b * tx) * (1 - ty) + (c * (1 - tx) + d * tx) * ty
        gx = ((b - a) * (1 - ty) + (d - c) * ty) / self.res
        gy = ((c - a) * (1 - tx) + (d - b) * tx) / self.res
        return v, gx, gy


def carved_tops(strands, crossings, C):
    """Each strand's top height at each sample: under strands dip near each
    crossing so their lowered top meets the over strand's flank in a straight
    line; over passes stay up (rules of make_celtic_band.py's CARVED style)."""
    k = math.tan(math.radians(C['slope']))
    drop = C['drop']
    slope = min(drop / max(C['ramp'], 1e-6), 0.8 * k)
    depth = [np.zeros(len(s.pts)) for s in strands]
    limit = [np.full(len(s.pts), np.inf) for s in strands]
    for oi, oj, ui, uj, sn in crossings:
        meet = strands[ui].width / 2 + drop / k
        full_u = meet / max(0.35, sn) + 0.05
        full_o = meet / max(0.35, sn)
        depth[ui] = np.maximum(depth[ui], drop - slope * np.maximum(0, strands[ui].along(uj) - full_u))
        limit[oi] = np.minimum(limit[oi], slope * np.maximum(0, strands[oi].along(oj) - full_o))
    return [C['depth'] - np.minimum(d, l) for d, l in zip(depth, limit)]


def carved_field(X, Y, x0, y0, dx, strands, tops, outline, C, walls=None, flat=None, raises=()):
    """Height (and slope) of the carved surface at points X, Y, arrays laid
    out like a grid of spacing dx from (x0, y0) (each point within a cell of
    its place). The background is the pocket floor (z = 0) rising to the rim
    near the edge (with panels, to the walls: the rim and the frame between
    them; flat: the panels that aren't carved, flat floor right up to the
    walls; raises: emblems, (outline, top), flat-topped with the carving
    sloping down from them); each strand has a flat top and flanks at
    C['slope']."""
    k = math.tan(math.radians(C['slope']))
    rim_top = C['rim_top']
    if walls is not None:
        d, sgx, sgy = walls.sample(X, Y)          # > 0 in the walls
        d, sgx, sgy = -d, -sgx, -sgy              # distance out of them
        slope = (d > 0) & (rim_top - k * d > 0)
        if flat is not None:
            slope &= flat.sample(X, Y)[0] < 0
        h = np.where(d <= 0, rim_top, np.where(slope, rim_top - k * d, 0.0))
        gx, gy = np.where(slope, -k * sgx, 0.0), np.where(slope, -k * sgy, 0.0)
    elif C['rim']:
        d, sgx, sgy = outline.sample(X, Y)
        wall = (d > C['rim_w']) & (rim_top - k * (d - C['rim_w']) > 0)
        h = np.where(d <= C['rim_w'], rim_top, np.maximum(0.0, rim_top - k * (d - C['rim_w'])))
        gx, gy = np.where(wall, -k * sgx, 0.0), np.where(wall, -k * sgy, 0.0)
    else:
        h, gx, gy = np.zeros_like(X), np.zeros_like(X), np.zeros_like(X)
    for emblem, top in raises:
        e, egx, egy = emblem.sample(X, Y)
        val = np.where(e >= 0, top, top + k * e)
        up = val > h
        h = np.where(up, val, h)
        gx = np.where(up & (e < 0), k * egx, np.where(up, 0.0, gx))
        gy = np.where(up & (e < 0), k * egy, np.where(up, 0.0, gy))
    nrows, ncols = X.shape
    half = reach = 0.0

    def window(lo, hi):
        c0, c1 = max(0, int(math.floor((lo[0] - reach - x0) / dx))), min(ncols - 1, int(math.ceil((hi[0] + reach - x0) / dx)))
        r0, r1 = max(0, int(math.floor((lo[1] - reach - y0) / dx))), min(nrows - 1, int(math.ceil((hi[1] + reach - y0) / dx)))
        return (slice(r0, r1 + 1), slice(c0, c1 + 1)) if c0 <= c1 and r0 <= r1 else None

    def raise_(w, mask, val, vx, vy):
        nonlocal h, gx, gy
        up = mask & (val > h[w])
        if up.any():
            for arr, new in ((h, val), (gx, vx), (gy, vy)):
                arr[w] = np.where(up, new, arr[w])

    for strand, top in zip(strands, tops):
        half = strand.width / 2
        hs = half * strand.scale             # (a tapered spiral: each sample's half width)
        reach = half + C['depth'] / k + dx
        P, n = strand.pts, len(strand.pts)
        # (an open strand, a spiral: no segment from its end back to its
        # start, and rounded ends)
        for a in range(n - 1 if strand.open else n):
            A, B, Z = P[a], P[(a + 1) % n], P[(a - 1) % n]
            ta, tb = top[a], top[(a + 1) % n]
            ha, hb = hs[a], hs[(a + 1) % n]
            e = B - A
            l2 = max(e @ e, 1e-12)
            # the strip beside the segment
            w = window(np.minimum(A, B), np.maximum(A, B))
            if w is not None:
                px, py = X[w], Y[w]
                t = ((px - A[0]) * e[0] + (py - A[1]) * e[1]) / l2
                lo = -np.inf if strand.open and a == 0 else -1e-6
                hi = np.inf if strand.open and a == n - 2 else 1 + 1e-6
                beside = (t >= lo) & (t <= hi)
                tc = np.clip(t, 0, 1)
                rx, ry = px - A[0] - tc * e[0], py - A[1] - tc * e[1]
                dd = np.hypot(rx, ry)
                hw = ha + (hb - ha) * tc
                fl = dd > hw
                f = np.where(fl, k / np.maximum(dd, 1e-9), 0.0)
                val = ta + (tb - ta) * tc - k * np.where(fl, dd - hw, 0.0)
                raise_(w, beside & (dd < reach), val,
                       (tb - ta) * e[0] / l2 - f * rx, (tb - ta) * e[1] / l2 - f * ry)
            # the wedge outside the bend at A
            if strand.open and a == 0:
                continue
            w = window(A, A)
            if w is not None:
                rx, ry = X[w] - A[0], Y[w] - A[1]
                p = A - Z
                wedge = (rx * p[0] + ry * p[1] > 0) & (rx * e[0] + ry * e[1] < 0)
                dd = np.hypot(rx, ry)
                fl = dd > ha
                f = np.where(fl, k / np.maximum(dd, 1e-9), 0.0)
                raise_(w, wedge & (dd < reach), ta - k * np.where(fl, dd - ha, 0.0), -f * rx, -f * ry)
    return h, gx, gy


def build_carved(strands, crossings, design, s):
    """The carved piece as one closed solid: the surface meshed by dual
    contouring (one point a cell, placed on the creases from the slopes),
    clipped to the outline, with walls down to a flat base plate_thickness
    below the pocket floor. strands: the carved ones."""
    C = dict(s['carve'])
    C['rim'], C['rim_w'] = s['rim'], s['rim_width']
    C.setdefault('rim_top', C['depth'] + s['rim_height'])
    k = math.tan(math.radians(C['slope']))
    res = CARVE_RESOLUTION
    outline = Outline(design['outline'], max(res, 0.05), (s['rim_width'] if s['rim'] else 0) + C['depth'] / k + 0.3)
    # with panels: the walls (the rim and the frame), and the flat panels
    walls = flat = None
    if s.get('walls') and design.get('rim_outline'):
        walls = Outline(design['rim_outline'], max(res, 0.05), C['rim_top'] / k + 0.3)
    if design.get('flat_outline'):
        flat = Outline(design['flat_outline'], max(res, 0.05), 0.3)
    raises = [(Outline(e['outline'], max(res, 0.05), e['top'] / k + 0.3), e['top'])
              for e in design.get('emblems', []) if e['carved']]
    tops = carved_tops(strands, crossings, C)

    allp = np.concatenate([np.asarray(g['outer'], float) for g in design['outline']])
    x0, y0 = allp.min(0) - 2 * res
    x1, y1 = allp.max(0) + 2 * res
    nx, ny = int(math.ceil((x1 - x0) / res)), int(math.ceil((y1 - y0) / res))
    print(f"Carving on a {nx} x {ny} grid ({res} mm)...")
    field = lambda X, Y, ox, oy: carved_field(X, Y, ox, oy, res, strands, tops, outline, C, walls, flat, raises)

    gxs, gys = np.meshgrid(x0 + np.arange(nx + 1) * res, y0 + np.arange(ny + 1) * res)
    H, GX, GY = field(gxs, gys, x0, y0)
    cxs, cys = np.meshgrid(x0 + (np.arange(nx) + 0.5) * res, y0 + (np.arange(ny) + 0.5) * res)
    Hc, GXc, GYc = field(cxs, cys, x0 + res / 2, y0 + res / 2)

    # one point a cell, where the planes through the samples best agree
    samples = [(gxs[di:di + ny, dj:dj + nx], gys[di:di + ny, dj:dj + nx], H[di:di + ny, dj:dj + nx],
                GX[di:di + ny, dj:dj + nx], GY[di:di + ny, dj:dj + nx]) for di in (0, 1) for dj in (0, 1)]
    samples.append((cxs, cys, Hc, GXc, GYc))
    ATA = np.zeros((ny, nx, 3, 3))
    ATb = np.zeros((ny, nx, 3))
    for px, py, ph, sx, sy in samples:
        w = 1.0 / (1.0 + sx * sx + sy * sy)
        row = np.stack([-sx, -sy, np.ones_like(sx)], -1)
        rhs = ph - sx * px - sy * py
        ATA += w[..., None, None] * row[..., :, None] * row[..., None, :]
        ATb += (w * rhs)[..., None] * row
    lam = 0.02
    ATA[..., 0, 0] += lam
    ATA[..., 1, 1] += lam
    ATb[..., 0] += lam * cxs
    ATb[..., 1] += lam * cys
    ATA[..., 2, 2] += 1e-9
    q = np.linalg.solve(ATA, ATb[..., None])[..., 0]
    qx = np.clip(q[..., 0], cxs - res / 2, cxs + res / 2)
    qy = np.clip(q[..., 1], cys - res / 2, cys + res / 2)

    # cells in the shape; the outermost ones put their point on the outline
    inside = outline.sample(cxs, cys)[0] >= 0
    pad = np.pad(inside, 1)
    edge = inside & ~(pad[:-2, 1:-1] & pad[2:, 1:-1] & pad[1:-1, :-2] & pad[1:-1, 2:])
    for _ in range(3):
        d, sx, sy = outline.sample(qx, qy)
        g2 = np.maximum(sx * sx + sy * sy, 1e-9)
        qx = np.where(edge, qx - d * sx / g2, qx)
        qy = np.where(edge, qy - d * sy / g2, qy)
    qz = field(qx, qy, x0 + res / 2, y0 + res / 2)[0]

    idx = -np.ones((ny, nx), int)
    idx[inside] = np.arange(inside.sum())
    top_v = np.stack([qx[inside], qy[inside], qz[inside]], 1)

    # faces round each grid point between four cells (three at an inner corner)
    a, b = idx[:-1, :-1], idx[:-1, 1:]          # cells below-left, below-right of the point
    c, d = idx[1:, 1:], idx[1:, :-1]            # above-right, above-left
    mid = H[1:ny, 1:nx]
    za, zb, zc, zd = (np.where(v >= 0, top_v[v, 2] if len(top_v) else 0, 0) for v in (a, b, c, d))
    all4 = (a >= 0) & (b >= 0) & (c >= 0) & (d >= 0)
    diag = np.abs((za + zc) / 2 - mid) <= np.abs((zb + zd) / 2 - mid)
    tris = [np.stack([a, b, c], -1)[all4 & diag], np.stack([a, c, d], -1)[all4 & diag],
            np.stack([a, b, d], -1)[all4 & ~diag], np.stack([b, c, d], -1)[all4 & ~diag]]
    three = ((a >= 0).astype(int) + (b >= 0) + (c >= 0) + (d >= 0)) == 3
    for keep in ((b, c, d), (a, c, d), (a, b, d), (a, b, c)):     # the missing one left out
        miss = [v for v in (a, b, c, d) if not any(v is kk for kk in keep)][0]
        tris.append(np.stack(keep, -1)[three & (miss < 0)])
    top_t = np.concatenate(tris).reshape(-1, 3)

    # base: the same triangles flat underneath; walls under the edge
    m = len(top_v)
    base = s['plate_thickness']
    bot_v = top_v.copy()
    bot_v[:, 2] = -base
    e = np.concatenate([top_t[:, [0, 1]], top_t[:, [1, 2]], top_t[:, [2, 0]]])
    key = np.minimum(e[:, 0], e[:, 1]) * m + np.maximum(e[:, 0], e[:, 1])
    _, inv, cnt = np.unique(key, return_inverse=True, return_counts=True)
    be = e[cnt[inv] == 1]
    walls = np.concatenate([np.stack([be[:, 0], be[:, 0] + m, be[:, 1] + m], 1),
                            np.stack([be[:, 0], be[:, 1] + m, be[:, 1]], 1)])
    verts = np.concatenate([top_v, bot_v])
    faces = np.concatenate([top_t, top_t[:, [0, 2, 1]] + m, walls])
    if signed_volume(verts, faces) < 0:
        faces = faces[:, [0, 2, 1]]
    return verts, faces


# ---------------------------------------------------------------- weighing
def planar_volume(objects):
    """Volume (mm3) however the parts overlap: rays down from above find the
    top surface, rays up from below the bottom, and the metal is between
    (the piece has no undercuts)."""
    from mathutils.bvhtree import BVHTree
    bm = bmesh.new()
    for o in objects:
        bm.from_mesh(o.data)
    if not bm.verts:
        bm.free()
        return 0.0
    co = np.array([v.co[:] for v in bm.verts])
    tree = BVHTree.FromBMesh(bm)
    bm.free()
    (x0, y0, z0), (x1, y1, z1) = co.min(0), co.max(0)
    step = max(0.02, math.sqrt((x1 - x0) * (y1 - y0) / WEIGH_SAMPLES))
    total = 0.0
    down, up = Vector((0, 0, -1)), Vector((0, 0, 1))
    for x in np.arange(x0 + step / 2, x1, step):
        for y in np.arange(y0 + step / 2, y1, step):
            hit = tree.ray_cast(Vector((x, y, z1 + 1)), down)
            if hit[0] is None:
                continue
            low = tree.ray_cast(Vector((x, y, z0 - 1)), up)
            if low[0] is not None:
                total += hit[0].z - low[0].z
    return total * step * step


def metal_report(volume, alloy):
    name, rho = ALLOYS.get(alloy, ALLOYS['STERLING'])
    others = ", ".join(f"{n} {volume / 1000 * r:.1f} g" for k, (n, r) in ALLOYS.items() if k != alloy)
    line = f"{volume:.0f} mm³: {volume / 1000 * rho:.2f} g of {name}"
    print("Metal: " + line)
    print("  in other materials: " + others)
    return line


def make_metal_material(name, color):
    mat = bpy.data.materials.get(name) or bpy.data.materials.new(name)
    if bpy.app.version < (5, 0, 0):      # always on from Blender 5
        mat.use_nodes = True
    mat.diffuse_color = (*color, 1.0)
    mat.metallic = 1.0
    mat.roughness = 0.15
    bsdf = mat.node_tree.nodes.get("Principled BSDF")
    if bsdf:
        bsdf.inputs["Base Color"].default_value = (*color, 1.0)
        bsdf.inputs["Metallic"].default_value = 1.0
        bsdf.inputs["Roughness"].default_value = 0.12
    return mat


# ---------------------------------------------------------------- build
def build(text):
    """Build the design (the page's export, as text). Returns the objects."""
    design = read_design(text)
    s = design['settings']
    us = bpy.context.scene.unit_settings
    us.system = 'METRIC'
    us.scale_length = 0.001
    us.length_unit = 'MILLIMETERS'
    for obj in [o for o in bpy.data.objects if o.name.startswith("Shape_")]:
        bpy.data.objects.remove(obj, do_unlink=True)
    metal = make_metal_material("Shape_Metal", METAL_COLORS.get(s['metal'], METAL_COLORS['SILVER']))

    strands = [Strand(d, SAMPLES_PER_SPAN, d.get('corners', s['corners']) == 'STRAIGHT', s) for d in design['strands']]
    objects = []
    style = s['style']
    # carved panels are carved into one solid (flat floor in the others,
    # which then stand on it as on a plate); the others' strands are swept
    carved = [st for st in strands if st.style == 'CARVED']
    swept = [st for st in strands if st.style != 'CARVED']
    if carved or style == 'CARVED':
        verts, faces = build_carved(carved, crossings_of(carved), design, s)
        objects.append(mesh_object("Shape_Carved", verts, faces, metal, math.radians(30)))
        volume = abs(signed_volume(verts, faces))
    if swept:
        crossings = crossings_of(swept)
        if carved:
            # (carved panels stand on a base, which the others share as their plate)
            for st in swept:
                st.on_plate, st.woven = True, False
        woven = any(st.woven for st in swept)
        heights = [swept_heights(st, crossings, si, st.height, st.dip, s['cell'], -0.05 if st.on_plate else 0.0, st.woven)
                   for si, st in enumerate(swept)]
        doubles = [si for si, st in enumerate(swept) if st.style == 'DOUBLE']
        parts = [sweep(st.pts, top, bot, st.width * st.scale, PROFILES.get(st.style, PROFILES['ROUNDED']), not st.open)
                 for st, (top, bot) in zip(swept, heights) if st.style != 'DOUBLE']
        if doubles:
            # (crossings are within a panel, so between Double strands only)
            index = {si: k for k, si in enumerate(doubles)}
            dcross = [(index[oi], oj, index[ui], uj, sn) for oi, oj, ui, uj, sn in crossings if oi in index and ui in index]
            parts += double_strands([swept[si] for si in doubles], dcross, [heights[si] for si in doubles], s)
        vs, fs, off = [], [], 0
        for v, f in parts:
            vs.append(v)
            fs.append(f + off)
            off += len(v)
        objects.append(mesh_object("Shape_Knot", np.concatenate(vs), np.concatenate(fs), metal, math.radians(40)))
        if not carved:
            # a piece for each thickness (panels can differ, or have none);
            # older designs have one plate, the outline's
            plates = design.get('plates')
            if plates is None:
                plates = [{'thickness': s['plate_thickness'], 'outline': design['outline']}] if s['plate'] else []
            for i, p in enumerate(plates):
                v, f = extrude(p['outline'], -p['thickness'], 0.0)
                objects.append(mesh_object("Shape_Plate" + (str(i + 1) if len(plates) > 1 else ""), v, f, metal, math.radians(30)))
        if s['rim'] and design.get('rim_outline') and not carved:
            # as high as the highest strands; woven: down to the strands' lowest
            low = max(st.dip * st.height / 2 for st in swept if st.woven) if woven else 0.05
            v, f = extrude(design['rim_outline'], -low, max(st.height for st in swept) + s['rim_height'])
            objects.append(mesh_object("Shape_Rim", v, f, metal, math.radians(30)))
        volume = planar_volume(objects)
    standing = [e for e in design.get('emblems', []) if not e['carved']]
    for i, e in enumerate(standing):
        v, f = extrude(e['outline'], -0.05, e['top'])
        objects.append(mesh_object(f"Shape_Emblem{i + 1}", v, f, metal, math.radians(30)))
    if standing:
        volume = planar_volume(objects)

    if SIMPLIFY_ANGLE > 0:
        for o in objects:
            n = len(o.data.polygons)
            mod = o.modifiers.new("simplify", 'DECIMATE')
            mod.decimate_type = 'DISSOLVE'
            mod.angle_limit = math.radians(SIMPLIFY_ANGLE)
            mod.use_dissolve_boundaries = False
            tri = o.modifiers.new("triangulate", 'TRIANGULATE')
            dg = bpy.context.evaluated_depsgraph_get()
            mesh = bpy.data.meshes.new_from_object(o.evaluated_get(dg))
            o.modifiers.remove(tri)
            o.modifiers.remove(mod)
            old = o.data
            o.data = mesh
            bpy.data.meshes.remove(old)
            print(f"{o.name}: simplified from {n} faces to {len(mesh.polygons)}")

    if MAKE_PRINT_SOLID:
        joined = bpy.data.meshes.new("Shape_Print")
        jbm = bmesh.new()
        for o in objects:
            jbm.from_mesh(o.data)
        jbm.to_mesh(joined)
        jbm.free()
        solid = bpy.data.objects.new("Shape_Print", joined)
        bpy.context.scene.collection.objects.link(solid)
        mod = solid.modifiers.new("fuse", 'REMESH')
        mod.mode = 'VOXEL'
        mod.voxel_size = PRINT_VOXEL
        dg = bpy.context.evaluated_depsgraph_get()
        fused = bpy.data.meshes.new_from_object(solid.evaluated_get(dg))
        solid.modifiers.remove(mod)
        old = solid.data
        solid.data = fused
        bpy.data.meshes.remove(old)
        solid.hide_set(True)
        solid.hide_render = True
        objects.append(solid)

    LAST_REPORT[0] = metal_report(volume, s.get('alloy', 'STERLING'))
    xs = np.concatenate([np.asarray(g['outer'])[:, 0] for g in design['outline']])
    ys = np.concatenate([np.asarray(g['outer'])[:, 1] for g in design['outline']])
    print(f"Celtic knot shape '{design.get('name', '')}': {style.lower()}, {len(strands)} strand(s), "
          f"{len(crossings_of(strands))} crossings, {xs.max() - xs.min():.1f} x {ys.max() - ys.min():.1f} mm")

    if not bpy.app.background:
        for area in [a for screen in bpy.data.screens for a in screen.areas if a.type == 'VIEW_3D']:
            for space in area.spaces:
                if space.type == 'VIEW_3D':
                    space.clip_start = 0.1
                    space.clip_end = 10000
        bpy.ops.object.select_all(action='DESELECT')
        for o in objects:
            if o.visible_get():
                o.select_set(True)
        bpy.context.view_layer.objects.active = objects[0]
    return objects


def clipboard_design():
    text = bpy.context.window_manager.clipboard or ''
    return text if text.lstrip().startswith(HEADER) else None


# ---------------------------------------------------------------- add-on
class WM_OT_knot_shape_done(bpy.types.Operator):
    """Shows what was built"""
    bl_idname = "wm.knot_shape_done"
    bl_label = "Celtic Knot Shape"

    message: bpy.props.StringProperty(name="Message", default="")

    def execute(self, context):
        return {'FINISHED'}

    def invoke(self, context, event):
        return context.window_manager.invoke_props_dialog(self, width=420)

    def draw(self, context):
        for i, line in enumerate(self.message.split("\n")):
            self.layout.label(text=line, icon='CHECKMARK' if i == 0 else 'NONE')


class MESH_OT_add_celtic_knot_shape(bpy.types.Operator):
    """Build a design from the Celtic Knot Shapes page (from the clipboard, or a saved file)"""
    bl_idname = "mesh.add_celtic_knot_shape"
    bl_label = "Celtic Knot Shape"
    bl_options = {'REGISTER', 'UNDO'}

    filepath: bpy.props.StringProperty(subtype='FILE_PATH', default='')
    filter_glob: bpy.props.StringProperty(default='*.json', options={'HIDDEN'})

    def invoke(self, context, event):
        if clipboard_design() or (DESIGN_FILE and os.path.exists(DESIGN_FILE)):
            return self.execute(context)
        context.window_manager.fileselect_add(self)
        return {'RUNNING_MODAL'}

    def execute(self, context):
        text, source = clipboard_design(), "the clipboard"
        if not text:
            path = self.filepath or DESIGN_FILE
            if not path or not os.path.exists(path):
                self.report({'ERROR'}, "No design: press Copy for Blender on the page, or choose a saved design file")
                return {'CANCELLED'}
            with open(path, encoding='utf-8') as f:
                text = f.read()
            source = os.path.basename(path)
        try:
            objects = build(text)
        except Exception as e:
            self.report({'ERROR'}, f"Celtic Knot Shape failed: {e}")
            raise
        if not bpy.app.background:
            bpy.ops.wm.knot_shape_done('INVOKE_DEFAULT',
                                       message=f"Built from {source}: {', '.join(o.name for o in objects)}\n{LAST_REPORT[0]}")
        return {'FINISHED'}


def menu_func(self, context):
    self.layout.operator(MESH_OT_add_celtic_knot_shape.bl_idname, icon='MOD_CURVE')


def register():
    bpy.utils.register_class(WM_OT_knot_shape_done)
    bpy.utils.register_class(MESH_OT_add_celtic_knot_shape)
    bpy.types.VIEW3D_MT_mesh_add.append(menu_func)


def unregister():
    bpy.types.VIEW3D_MT_mesh_add.remove(menu_func)
    bpy.utils.unregister_class(MESH_OT_add_celtic_knot_shape)
    bpy.utils.unregister_class(WM_OT_knot_shape_done)


if __name__ == "__main__":
    args = sys.argv[sys.argv.index('--') + 1:] if '--' in sys.argv else []
    if args:
        # headless: build a saved design and optionally write an STL
        with open(args[0], encoding='utf-8') as f:
            objs = build(f.read())
        if len(args) > 1:
            for o in bpy.context.scene.objects:
                o.select_set(o in objs and o.visible_get())
            if hasattr(bpy.ops.wm, 'stl_export'):
                bpy.ops.wm.stl_export(filepath=args[1], export_selected_objects=True)
            else:
                bpy.ops.export_mesh.stl(filepath=args[1], use_selection=True)
            print("Wrote " + args[1])
    else:
        register()