#include <algorithm>
#include <array>
#include <atomic>
#include <bit>
#include <cassert>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <functional>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <limits>
#include <numeric>
#include <span>
#include <string>
#include <thread>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>

#define WIN32_LEAN_AND_MEAN
#define NOMINMAX
#include <windows.h>
#include <psapi.h>
#pragma comment(lib, "Psapi.lib")

namespace {

using Clock = std::chrono::steady_clock;
constexpr int kCoreCells = 256;
constexpr int kHalo = 1;
constexpr int kNodes = kCoreCells + 1 + 2 * kHalo;
constexpr int kWordsPerLine = (kNodes + 63) / 64;
constexpr int kPageCells = 64;
constexpr int kPagesPerAxis = kCoreCells / kPageCells;
constexpr int kBucketCm = 400;
constexpr int kBucketCount = 17;
constexpr int kShapeLimitPerLine = 64;
constexpr int64_t kQScale = int64_t{1} << 28;
constexpr int64_t kSubPerCm = 16; // canonical interval/sign lattice: 0.625 mm
constexpr int64_t kSurfaceSubPerCm = 32; // crossing positions: 0.3125 mm
constexpr int kCellStepCm = 25;
constexpr int kBlendWidthCm = 70;
constexpr float kReliefBoundMeters = 0.60f;
constexpr int kGridMinCm = -3225;
constexpr int kGridMaxCm = kGridMinCm + (kNodes - 1) * kCellStepCm;
constexpr int kWorldOriginCm[3] = {-12000, 6000, -4500};
constexpr uint64_t kMask64 = ~uint64_t{0};

enum class Operation : uint8_t { CarveAir, AddRock };
enum class FeatureKind : uint8_t {
    Room, Tunnel, Pillar, Stalactite, Stalagmite, Shelf, FloatingMass,
    Stratum, Crater, CrowdedJunction
};

struct Mat3Q15 { int32_t v[3][3]; };

// Quantized rotations are fixed data, so corpus construction does not depend on libm.
constexpr std::array<Mat3Q15, 8> kRotations = {{
    {{{32768,0,0},{0,32768,0},{0,0,32768}}},
    {{{23170,-23170,0},{23170,23170,0},{0,0,32768}}},
    {{{23170,0,23170},{0,32768,0},{-23170,0,23170}}},
    {{{32768,0,0},{0,23170,-23170},{0,23170,23170}}},
    {{{0,-32768,0},{32768,0,0},{0,0,32768}}},
    {{{0,0,32768},{0,32768,0},{-32768,0,0}}},
    {{{32768,0,0},{0,0,-32768},{0,32768,0}}},
    {{{16384,-23170,16384},{16384,23170,16384},{-23170,0,23170}}}
}};

struct Shape {
    int32_t c[3]{};          // local centimetres
    int32_t radius[3]{};     // local centimetres
    int32_t rot[3][3]{};     // Q15 rotation, columns are local axes
    int64_t q[3][3]{};       // Q28 / cm^2, exact predicate coefficients
    int32_t lo[3]{}, hi[3]{}; // conservative local-centimetre AABB, inclusive
    Operation op = Operation::CarveAir;
    FeatureKind kind = FeatureKind::Room;
    uint32_t id = 0;
    bool auditThin = false;
};

struct Crossing {
    int32_t positionSub32 = 0; // local cm * 32, half of the canonical 1/16 cm grid
    uint32_t shapeId = 0;
    uint8_t beforeRock = 0;
    uint8_t afterRock = 0;
};
constexpr uint32_t kBlendCrossingShape = std::numeric_limits<uint32_t>::max();

struct LineRange {
    uint32_t offset = 0;
    uint16_t count = 0;
    uint16_t candidates = 0;
    uint16_t primitiveIntervals = 0;
    uint16_t occupancyRuns = 0;
};

struct AxisData {
    std::vector<uint64_t> signs;
    std::vector<LineRange> lines;
    std::vector<Crossing> crossings;
};

struct AxisStats {
    double candidateP50 = 0, candidateP95 = 0;
    double intervalP50 = 0, intervalP95 = 0;
    double crossingP50 = 0, crossingP95 = 0;
    double runsP50 = 0, runsP95 = 0;
    uint32_t candidateMax = 0, intervalMax = 0, crossingMax = 0, runsMax = 0;
    uint64_t candidateTotal = 0, intervalTotal = 0, crossingTotal = 0, runsTotal = 0;
};

struct Vertex {
    float x, y, z;
    float nx, ny, nz;
};

struct TempVertex {
    uint64_t key = 0;
    uint32_t ownerPage = 0;
    Vertex v{};
};

struct PageRange { uint32_t page, indexStart, indexCount, vertexStart, vertexCount; };

struct MeshOutput {
    std::vector<TempVertex> vertices;
    std::vector<uint32_t> indices;
    std::vector<PageRange> pages;
    uint64_t hash = 0;
    uint64_t bytes = 0;
    uint64_t triangles = 0;
    uint64_t ambiguousFaces = 0;
    uint64_t multiCrossEdges = 0;
    uint64_t missingAnalyticCrossings = 0;
    uint64_t contourDegreeFaults = 0;
    uint64_t openContourLoops = 0;
    uint64_t explicitThinTriangles = 0;
    uint64_t thinFeaturesDetected = 0;
    uint64_t thinFeaturesHiddenInLattice = 0;
    uint64_t thinFeaturesMeshed = 0;
    uint64_t pageBorderVertices = 0;
    uint64_t qefClampedVertices = 0;
    uint64_t cellContourFaults = 0;
    uint64_t cachedHermiteFallbacks = 0;
    uint64_t nonManifoldEdges = 0;
    uint64_t boundaryEdges = 0;
    uint64_t boundaryThinEdges = 0;
    uint64_t boundaryNearGridEdge = 0;
    uint64_t degenerateTriangles = 0;
    uint64_t overusedEdges = 0;
    uint64_t displacedHash = 0;
    double maxDisplacementMeters = 0.0;
};

struct Bin3D { std::vector<uint32_t> ids; };
struct BlendBox { int32_t lo[3]{}, hi[3]{}; };
struct ProjectionBins { std::vector<uint32_t> ids[kBucketCount * kBucketCount]; };

struct Corpus {
    std::vector<Shape> shapes;
    Bin3D bins3d[kBucketCount * kBucketCount * kBucketCount];
    Bin3D blendBins3d[kBucketCount * kBucketCount * kBucketCount];
    ProjectionBins projections[3];
    std::vector<BlendBox> blendBoxes;
    std::vector<uint8_t> blendNodeMask;
    std::vector<uint32_t> blendNodeIndices;
    uint64_t blendZoneNodes = 0;
    uint64_t blendZoneCells = 0;
    double blendZoneBuildMs = 0.0;
    uint64_t seed = 0;
    uint32_t maxMembershipAudit = 0;
};

struct WorkerLines {
    std::vector<Crossing> crossings;
};

struct FarOutput {
    std::vector<Vertex> vertices;
    std::vector<uint32_t> indices;
    std::vector<PageRange> patches;
    std::array<uint64_t, 6> lodNodes{};
    uint64_t nodes = 0, descendants = 0, patchCount = 0, bytes = 0, hash = 0;
};

double ms(Clock::time_point a, Clock::time_point b) {
    return std::chrono::duration<double, std::milli>(b - a).count();
}

uint64_t mix64(uint64_t x) {
    x += 0x9E3779B97F4A7C15ull;
    x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9ull;
    x = (x ^ (x >> 27)) * 0x94D049BB133111EBull;
    return x ^ (x >> 31);
}

struct Rng {
    uint64_t state;
    explicit Rng(uint64_t s) : state(s) {}
    uint32_t next() { state += 0x9E3779B97F4A7C15ull; return static_cast<uint32_t>(mix64(state) >> 32); }
    int range(int lo, int hi) { return lo + static_cast<int>(next() % static_cast<uint32_t>(hi - lo + 1)); }
    bool chance(int n, int d) { return static_cast<int>(next() % d) < n; }
};

int64_t round_div(int64_t n, int64_t d) {
    if (n >= 0) return (n + d / 2) / d;
    return -((-n + d / 2) / d);
}

int floor_div(int a, int b) {
    const int q = a / b;
    const int r = a % b;
    return q - ((r != 0 && ((r < 0) != (b < 0))) ? 1 : 0);
}

int64_t floor_div64(int64_t a, int64_t b) {
    const int64_t q = a / b;
    const int64_t r = a % b;
    return q - ((r != 0 && ((r < 0) != (b < 0))) ? 1 : 0);
}

int ceil_div(int a, int b) { return -floor_div(-a, b); }

void build_quadric(Shape& s, uint32_t rotationIndex) {
    const auto& r = kRotations[rotationIndex % kRotations.size()];
    for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) s.rot[i][j] = r.v[i][j];
    // Q = R * diag(1/r^2) * R^T. Each term is rounded once to Q28.
    for (int a = 0; a < 3; ++a) for (int b = 0; b < 3; ++b) {
        int64_t sum = 0;
        for (int k = 0; k < 3; ++k) {
            const int64_t prod = static_cast<int64_t>(s.rot[a][k]) * s.rot[b][k];
            const int64_t den = (int64_t{1} << 30) * s.radius[k] * s.radius[k];
            sum += round_div(prod * kQScale, den);
        }
        s.q[a][b] = sum;
    }
    for (int axis = 0; axis < 3; ++axis) {
        int64_t extentNumerator = 0;
        for (int k = 0; k < 3; ++k) {
            const int64_t v = std::abs(s.rot[axis][k]);
            extentNumerator += v * s.radius[k];
        }
        const int32_t extent = static_cast<int32_t>((extentNumerator + 32767) / 32768) + 3;
        s.lo[axis] = s.c[axis] - extent;
        s.hi[axis] = s.c[axis] + extent;
    }
}

void add_shape(Corpus& c, int x, int y, int z, int rx, int ry, int rz,
               Operation op, FeatureKind kind, uint32_t rotation, bool auditThin = false) {
    Shape s;
    s.c[0] = x; s.c[1] = y; s.c[2] = z;
    s.radius[0] = std::max(1, rx); s.radius[1] = std::max(1, ry); s.radius[2] = std::max(1, rz);
    s.op = op; s.kind = kind; s.id = static_cast<uint32_t>(c.shapes.size()); s.auditThin = auditThin;
    build_quadric(s, rotation);
    c.shapes.push_back(s);
}

void add_sphere(Corpus& c, int x, int y, int z, int r, Operation op, FeatureKind kind) {
    add_shape(c, x, y, z, r, r, r, op, kind, 0);
}

int64_t eval_quad(const Shape& s, const int64_t d[3]);

void build_blend_zones(Corpus& c) {
    const int pad = (kBlendWidthCm + 1) / 2 + kCellStepCm;
    c.blendBoxes.clear();
    for (size_t i = 0; i < c.shapes.size(); ++i) {
        const Shape& a = c.shapes[i];
        if (a.kind == FeatureKind::Stratum) continue;
        for (size_t j = i + 1; j < c.shapes.size(); ++j) {
            const Shape& b = c.shapes[j];
            if (b.kind == FeatureKind::Stratum) continue;
            BlendBox box{};
            bool intersects = true;
            for (int axis = 0; axis < 3; ++axis) {
                box.lo[axis] = std::max(a.lo[axis] - pad, b.lo[axis] - pad);
                box.hi[axis] = std::min(a.hi[axis] + pad, b.hi[axis] + pad);
                box.lo[axis] = std::max(box.lo[axis], kGridMinCm);
                box.hi[axis] = std::min(box.hi[axis], kGridMaxCm);
                intersects &= box.lo[axis] <= box.hi[axis];
            }
            if (intersects) c.blendBoxes.push_back(box);
        }
    }

    const size_t rowCount = static_cast<size_t>(kNodes) * kNodes;
    std::vector<std::array<uint64_t,kWordsPerLine>> nodeRows(rowCount);
    auto setRange = [](auto& words, int lo, int hi) {
        const int firstWord=lo>>6, lastWord=hi>>6;
        for (int w=firstWord; w<=lastWord; ++w) {
            const int firstBit=(w==firstWord)?(lo&63):0;
            const int lastBit=(w==lastWord)?(hi&63):63;
            const uint64_t low=~uint64_t{0}<<firstBit;
            const uint64_t high=lastBit==63?~uint64_t{0}:((uint64_t{1}<<(lastBit+1))-1);
            words[w] |= low & high;
        }
    };
    for (const BlendBox& box : c.blendBoxes) {
        int nlo[3], nhi[3];
        for (int axis = 0; axis < 3; ++axis) {
            nlo[axis] = std::clamp(ceil_div(box.lo[axis] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
            nhi[axis] = std::clamp(floor_div(box.hi[axis] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
        }
        if (nlo[0] <= nhi[0] && nlo[1] <= nhi[1] && nlo[2] <= nhi[2])
            for (int z = nlo[2]; z <= nhi[2]; ++z) for (int y = nlo[1]; y <= nhi[1]; ++y)
                setRange(nodeRows[static_cast<size_t>(z) * kNodes + y], nlo[0], nhi[0]);
    }

    c.blendNodeMask.assign(static_cast<size_t>(kNodes) * kNodes * kNodes, 0);
    c.blendNodeIndices.clear();
    c.blendNodeIndices.reserve(300000);
    c.blendZoneNodes = 0;
    for (size_t row = 0; row < rowCount; ++row) {
        const size_t base=row*kNodes;
        for (int w=0; w<kWordsPerLine; ++w) {
            uint64_t bits=nodeRows[row][w];
            while (bits) {
                const int x=w*64+static_cast<int>(std::countr_zero(bits)); bits&=bits-1;
                if (x>=kNodes) continue;
                const size_t index=base+x;
                c.blendNodeMask[index]=1;
                c.blendNodeIndices.push_back(static_cast<uint32_t>(index));
                ++c.blendZoneNodes;
            }
        }
    }
    c.blendZoneCells = 0;
}

void filter_blend_zone_mask(Corpus& c);
size_t node_mask_index(int x, int y, int z);

void build_corpus(Corpus& c, uint64_t seed) {
    c.seed = seed;
    Rng rng(seed ^ 0x0A57A20261002ull);
    std::array<std::array<std::array<std::array<int, 3>, 3>, 3>, 2> centers{};
    for (int layer = 0; layer < 2; ++layer) {
        for (int iy = 0; iy < 3; ++iy) for (int ix = 0; ix < 3; ++ix) {
            const int x = (ix - 1) * 1300 + rng.range(-170, 170);
            const int y = (iy - 1) * 1300 + rng.range(-170, 170);
            const int z = (layer == 0 ? -850 : 850) + rng.range(-120, 120);
            centers[layer][iy][ix] = {x, y, z};
            const int rx = rng.range(480, 820);
            const int ry = rng.range(430, 760);
            const int rz = rng.range(300, 540);
            add_shape(c, x, y, z, rx, ry, rz, Operation::CarveAir, FeatureKind::Room,
                      rng.next() % 8);
        }
    }

    // Curved tunnel runs are unions of overlapping analytic ellipsoidal segments.
    for (int layer = 0; layer < 2; ++layer) {
        for (int iy = 0; iy < 3; ++iy) for (int ix = 0; ix < 3; ++ix) {
            const auto a = centers[layer][iy][ix];
            for (int dir = 0; dir < 2; ++dir) {
                const int nx = ix + (dir == 0);
                const int ny = iy + (dir == 1);
                if (nx >= 3 || ny >= 3) continue;
                const auto b = centers[layer][ny][nx];
                const int cx = (a[0] + b[0]) / 2 + rng.range(-220, 220);
                const int cy = (a[1] + b[1]) / 2 + rng.range(-220, 220);
                const int cz = (a[2] + b[2]) / 2 + rng.range(-150, 150);
                for (int seg = 0; seg < 5; ++seg) {
                    const int t = 1 + seg;
                    const int x = (a[0] * (6 - t) + b[0] * t) / 6 + (seg == 2 ? cx - (a[0] + b[0]) / 2 : 0);
                    const int y = (a[1] * (6 - t) + b[1] * t) / 6 + (seg == 2 ? cy - (a[1] + b[1]) / 2 : 0);
                    const int z = (a[2] * (6 - t) + b[2] * t) / 6 + (seg == 2 ? cz - (a[2] + b[2]) / 2 : 0);
                    const int r = rng.range(210, 295);
                    // Adjacent ellipsoids use changing cross-sections and axis orientation to
                    // avoid a chain of identical spherical tunnel beads.
                    add_shape(c, x, y, z, r, r * rng.range(82, 115) / 100,
                              r * rng.range(78, 112) / 100, Operation::CarveAir,
                              FeatureKind::Tunnel, rng.next() % 8);
                }
            }
        }
    }
    // Vertical connectors between stacked room layers.
    for (int iy = 0; iy < 3; ++iy) for (int ix = 0; ix < 3; ++ix) {
        const auto a = centers[0][iy][ix], b = centers[1][iy][ix];
        for (int seg = 1; seg <= 3; ++seg) {
            const int r = rng.range(195, 265);
            add_shape(c, (a[0] * (4 - seg) + b[0] * seg) / 4,
                       (a[1] * (4 - seg) + b[1] * seg) / 4,
                       (a[2] * (4 - seg) + b[2] * seg) / 4,
                       r, r * rng.range(80, 115) / 100, r * rng.range(78, 112) / 100,
                       Operation::CarveAir, FeatureKind::Tunnel, rng.next() % 8);
        }
    }

    // Dense constructive formations, placed in chamber space or near chamber walls.
    for (int i = 0; i < 40; ++i) {
        const int x = rng.range(-2500, 2500), y = rng.range(-2500, 2500);
        const int z = rng.range(-1900, 1900);
        add_shape(c, x, y, z, rng.range(70, 150), rng.range(70, 150), rng.range(260, 520),
                  Operation::AddRock, FeatureKind::Pillar, rng.next() % 8);
    }
    for (int i = 0; i < 36; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        const int sign = (i & 1) ? 1 : -1;
        add_shape(c, center[0] + rng.range(-260, 260), center[1] + rng.range(-260, 260),
                  center[2] + sign * rng.range(260, 410), rng.range(45, 85),
                  rng.range(45, 85), rng.range(230, 370), Operation::AddRock,
                  (i & 1) ? FeatureKind::Stalagmite : FeatureKind::Stalactite, rng.next() % 8);
    }
    for (int i = 0; i < 42; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        add_shape(c, center[0] + rng.range(-350, 350), center[1] + rng.range(-350, 350),
                  center[2] + rng.range(-220, 220), rng.range(250, 650), rng.range(130, 320),
                  rng.range(30, 90), Operation::AddRock, FeatureKind::Shelf, rng.next() % 8);
    }
    for (int i = 0; i < 24; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        add_shape(c, center[0] + rng.range(-280, 280), center[1] + rng.range(-280, 280),
                  center[2] + rng.range(-250, 250), rng.range(180, 390), rng.range(150, 340),
                  rng.range(140, 330), Operation::AddRock, FeatureKind::FloatingMass, rng.next() % 8);
    }
    // Several crater cuts, including deliberately intersecting crater groups.
    for (int i = 0; i < 12; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        const int x = center[0] + rng.range(300, 650);
        const int y = center[1] + rng.range(-280, 280);
        const int z = center[2] + rng.range(-280, 280);
        add_shape(c, x, y, z, rng.range(230, 520), rng.range(210, 460), rng.range(180, 420),
                  Operation::CarveAir, FeatureKind::Crater, rng.next() % 8);
    }
    // Crowded junction: a bounded 40-shape overlap cluster exercises long event lists.
    for (int i = 0; i < 40; ++i) {
        const int x = rng.range(-110, 110), y = rng.range(-110, 110), z = rng.range(-110, 110);
        const Operation op = (i % 3 == 0) ? Operation::AddRock : Operation::CarveAir;
        const int r = rng.range(75, 190);
        add_shape(c, x, y, z, r, rng.range(70, 200), rng.range(75, 185), op,
                  FeatureKind::CrowdedJunction, rng.next() % 8);
    }
    // Three widely separated, differently oriented rock bands replace the repeated thin stack.
    auto snapToLattice = [](int value) {
        const int node = static_cast<int>(std::llround((value - kGridMinCm) / static_cast<double>(kCellStepCm)));
        return kGridMinCm + node * kCellStepCm;
    };
    auto isAirBeforeStratum = [&](int x, int y, int z) {
        const int64_t p[3] = {static_cast<int64_t>(x)*kSubPerCm,
                              static_cast<int64_t>(y)*kSubPerCm,
                              static_cast<int64_t>(z)*kSubPerCm};
        bool rock = true;
        for (const Shape& s : c.shapes) {
            bool insideAabb = true;
            for (int a=0;a<3;++a) insideAabb &= p[a]>=static_cast<int64_t>(s.lo[a])*kSubPerCm
                && p[a]<=static_cast<int64_t>(s.hi[a])*kSubPerCm;
            if (!insideAabb) continue;
            const int64_t d[3] = {p[0]-static_cast<int64_t>(s.c[0])*kSubPerCm,
                                  p[1]-static_cast<int64_t>(s.c[1])*kSubPerCm,
                                  p[2]-static_cast<int64_t>(s.c[2])*kSubPerCm};
            if (eval_quad(s,d) <= kQScale*kSubPerCm*kSubPerCm) rock = s.op==Operation::AddRock;
        }
        return !rock;
    };
    for (int i = 0; i < 3; ++i) {
        const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3];
        const int desired[3] = {center[0] + rng.range(-180,180),center[1] + rng.range(-180,180),
                                center[2] + rng.range(-180,180)};
        int placed[3] = {snapToLattice(desired[0]),snapToLattice(desired[1]),snapToLattice(desired[2])};
        if (!isAirBeforeStratum(placed[0],placed[1],placed[2])) {
            bool found=false;
            for (int radius=25;radius<=300&&!found;radius+=25) {
                const int offsets[6][3]={{radius,0,0},{-radius,0,0},{0,radius,0},{0,-radius,0},{0,0,radius},{0,0,-radius}};
                for(const auto& offset:offsets){
                    int candidate[3]={snapToLattice(desired[0]+offset[0]),snapToLattice(desired[1]+offset[1]),
                                       snapToLattice(desired[2]+offset[2])};
                    if(isAirBeforeStratum(candidate[0],candidate[1],candidate[2])){std::copy(candidate,candidate+3,placed);found=true;break;}
                }
            }
        }
        add_shape(c, placed[0], placed[1], placed[2], rng.range(300, 580), rng.range(250, 520),
                  rng.range(7, 12), Operation::AddRock, FeatureKind::Stratum,
                  std::array<uint32_t,3>{4,2,3}[i]);
    }
    // A guaranteed isolated sub-cell air seam in solid rock. No 25 cm lattice node lies inside it.
    add_shape(c, 2900, -2900, 2887, 500, 400, 6, Operation::CarveAir,
              FeatureKind::Stratum, 0, true);

    const auto blendBuildBegin = Clock::now();
    build_blend_zones(c);
    c.blendZoneBuildMs = ms(blendBuildBegin, Clock::now());

    // Build deterministic 3D and axis-projection bins in feature order.
    for (uint32_t id = 0; id < c.shapes.size(); ++id) {
        const Shape& s = c.shapes[id];
        const int bx0 = std::clamp(floor_div(s.lo[0] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int bx1 = std::clamp(floor_div(s.hi[0] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int by0 = std::clamp(floor_div(s.lo[1] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int by1 = std::clamp(floor_div(s.hi[1] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int bz0 = std::clamp(floor_div(s.lo[2] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int bz1 = std::clamp(floor_div(s.hi[2] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        for (int bz = bz0; bz <= bz1; ++bz) for (int by = by0; by <= by1; ++by)
            for (int bx = bx0; bx <= bx1; ++bx)
                c.bins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids.push_back(id);
        const int blendLoX = std::clamp(floor_div(s.lo[0] - kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int blendHiX = std::clamp(floor_div(s.hi[0] + kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int blendLoY = std::clamp(floor_div(s.lo[1] - kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int blendHiY = std::clamp(floor_div(s.hi[1] + kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int blendLoZ = std::clamp(floor_div(s.lo[2] - kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int blendHiZ = std::clamp(floor_div(s.hi[2] + kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        for (int bz = blendLoZ; bz <= blendHiZ; ++bz) for (int by = blendLoY; by <= blendHiY; ++by)
            for (int bx = blendLoX; bx <= blendHiX; ++bx)
                c.blendBins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids.push_back(id);
        for (int axis = 0; axis < 3; ++axis) {
            int a = (axis + 1) % 3, b = (axis + 2) % 3;
            int a0 = std::clamp(floor_div(s.lo[a] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            int a1 = std::clamp(floor_div(s.hi[a] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            int b0 = std::clamp(floor_div(s.lo[b] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            int b1 = std::clamp(floor_div(s.hi[b] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            for (int ib = b0; ib <= b1; ++ib) for (int ia = a0; ia <= a1; ++ia)
                c.projections[axis].ids[ib * kBucketCount + ia].push_back(id);
        }
    }
    const auto zoneFilterBegin = Clock::now();
    filter_blend_zone_mask(c);
    c.blendZoneBuildMs += ms(zoneFilterBegin, Clock::now());
}

int bucket_of_sub(int64_t qSub) {
    const int64_t qCmFloor = floor_div64(qSub, kSubPerCm);
    return std::clamp(floor_div(static_cast<int>(qCmFloor) - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
}

int64_t eval_quad(const Shape& s, const int64_t d[3]) {
    int64_t value = 0;
    value += s.q[0][0] * d[0] * d[0];
    value += s.q[1][1] * d[1] * d[1];
    value += s.q[2][2] * d[2] * d[2];
    value += 2 * s.q[0][1] * d[0] * d[1];
    value += 2 * s.q[0][2] * d[0] * d[2];
    value += 2 * s.q[1][2] * d[1] * d[2];
    return value;
}

double shape_phi_approx(const Shape& s, const int64_t pSub[3]) {
    int64_t d[3] = {pSub[0] - static_cast<int64_t>(s.c[0]) * kSubPerCm,
                    pSub[1] - static_cast<int64_t>(s.c[1]) * kSubPerCm,
                    pSub[2] - static_cast<int64_t>(s.c[2]) * kSubPerCm};
    const int64_t q = eval_quad(s, d);
    const double q2 = std::max(0.0, static_cast<double>(q)
        / static_cast<double>(kQScale * kSubPerCm * kSubPerCm));
    const double radius = static_cast<double>(std::min({s.radius[0], s.radius[1], s.radius[2]}));
    return (1.0 - std::sqrt(q2)) * radius;
}

void filter_blend_zone_mask(Corpus& c) {
    c.blendZoneNodes = 0;
    std::vector<uint32_t> candidates = std::move(c.blendNodeIndices);
    c.blendNodeIndices.clear();
    c.blendNodeIndices.reserve(candidates.size() / 4 + 1024);
    for (uint32_t rawIndex : candidates) {
        const size_t nodeIndex = rawIndex;
        const int x = static_cast<int>(nodeIndex % kNodes);
        const int y = static_cast<int>((nodeIndex / kNodes) % kNodes);
        const int z = static_cast<int>(nodeIndex / (kNodes * kNodes));
        uint8_t& marked = c.blendNodeMask[nodeIndex];
        if (!marked) continue;
        const int64_t p[3] = {
            static_cast<int64_t>(kGridMinCm + x * kCellStepCm) * kSubPerCm,
            static_cast<int64_t>(kGridMinCm + y * kCellStepCm) * kSubPerCm,
            static_cast<int64_t>(kGridMinCm + z * kCellStepCm) * kSubPerCm};
        const int bx = std::clamp(floor_div(kGridMinCm + x * kCellStepCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int by = std::clamp(floor_div(kGridMinCm + y * kCellStepCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        const int bz = std::clamp(floor_div(kGridMinCm + z * kCellStepCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
        uint32_t nearSurfaces = 0;
        const auto& candidates = c.blendBins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids;
        for (uint32_t id : candidates) {
            const Shape& s = c.shapes[id];
            if (s.kind == FeatureKind::Stratum) continue;
            bool insideExpandedAabb=true;
            for (int axis=0;axis<3;++axis)
                insideExpandedAabb &= p[axis]>=static_cast<int64_t>(s.lo[axis]-kBlendWidthCm)*kSubPerCm
                    && p[axis]<=static_cast<int64_t>(s.hi[axis]+kBlendWidthCm)*kSubPerCm;
            if (!insideExpandedAabb) continue;
            if (std::abs(shape_phi_approx(s, p)) <= 0.5 * kBlendWidthCm && ++nearSurfaces >= 2) break;
        }
        marked = nearSurfaces >= 2 ? 1 : 0;
        c.blendZoneNodes += marked;
        if (marked) c.blendNodeIndices.push_back(static_cast<uint32_t>(nodeIndex));
    }
    std::vector<uint8_t> cellMask(static_cast<size_t>(kCoreCells) * kCoreCells * kCoreCells, 0);
    for (uint32_t nodeIndex : c.blendNodeIndices) {
        const int x = static_cast<int>(nodeIndex % kNodes);
        const int y = static_cast<int>((nodeIndex / kNodes) % kNodes);
        const int z = static_cast<int>(nodeIndex / (kNodes * kNodes));
        for (int cz = z - 1; cz <= z; ++cz) for (int cy = y - 1; cy <= y; ++cy) for (int cx = x - 1; cx <= x; ++cx) {
            if (cx < kHalo || cy < kHalo || cz < kHalo
                || cx >= kHalo + kCoreCells || cy >= kHalo + kCoreCells || cz >= kHalo + kCoreCells) continue;
            const size_t index = (static_cast<size_t>(cx - kHalo) * kCoreCells + (cy - kHalo)) * kCoreCells + (cz - kHalo);
            cellMask[index] = 1;
        }
    }
    c.blendZoneCells = std::count(cellMask.begin(), cellMask.end(), uint8_t{1});
}

bool contains_sub(const Shape& s, const int64_t pSub[3]) {
    int64_t d[3] = {pSub[0] - static_cast<int64_t>(s.c[0]) * kSubPerCm,
                    pSub[1] - static_cast<int64_t>(s.c[1]) * kSubPerCm,
                    pSub[2] - static_cast<int64_t>(s.c[2]) * kSubPerCm};
    return eval_quad(s, d) <= kQScale * kSubPerCm * kSubPerCm;
}

bool shape_aabb_contains_sub(const Shape& s, const int64_t pSub[3]) {
    for (int a = 0; a < 3; ++a) {
        if (pSub[a] < static_cast<int64_t>(s.lo[a]) * kSubPerCm
            || pSub[a] > static_cast<int64_t>(s.hi[a]) * kSubPerCm) return false;
    }
    return true;
}

bool point_rock(const Corpus& c, const int64_t pSub[3]) {
    const int xcm = static_cast<int>(floor_div(static_cast<int>(pSub[0]), static_cast<int>(kSubPerCm)));
    const int ycm = static_cast<int>(floor_div(static_cast<int>(pSub[1]), static_cast<int>(kSubPerCm)));
    const int zcm = static_cast<int>(floor_div(static_cast<int>(pSub[2]), static_cast<int>(kSubPerCm)));
    const int bx = std::clamp(floor_div(xcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const int by = std::clamp(floor_div(ycm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const int bz = std::clamp(floor_div(zcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    bool rock = true;
    const auto& bucket = c.bins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids;
    for (uint32_t id : bucket) {
        const Shape& s = c.shapes[id];
        if (shape_aabb_contains_sub(s, pSub) && contains_sub(s, pSub)) rock = (s.op == Operation::AddRock);
    }
    return rock;
}

double smooth_min(double a, double b, double k) {
    const double h = std::max(k - std::abs(a - b), 0.0) / k;
    return std::min(a, b) - h * h * k * 0.25;
}

double smooth_max(double a, double b, double k) {
    return -smooth_min(-a, -b, k);
}

double blended_field(const Corpus& c, const int64_t pSub[3]) {
    const int xcm = static_cast<int>(floor_div(static_cast<int>(pSub[0]), static_cast<int>(kSubPerCm)));
    const int ycm = static_cast<int>(floor_div(static_cast<int>(pSub[1]), static_cast<int>(kSubPerCm)));
    const int zcm = static_cast<int>(floor_div(static_cast<int>(pSub[2]), static_cast<int>(kSubPerCm)));
    const int bx = std::clamp(floor_div(xcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const int by = std::clamp(floor_div(ycm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const int bz = std::clamp(floor_div(zcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const auto& bucket = c.blendBins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids;
    double field = 100000.0; // positive means rock; the untouched world is solid rock
    for (uint32_t id : bucket) {
        const Shape& s = c.shapes[id];
        bool insideExpandedAabb = true;
        for (int axis = 0; axis < 3; ++axis)
            insideExpandedAabb &= pSub[axis] >= static_cast<int64_t>(s.lo[axis] - kBlendWidthCm) * kSubPerCm
                 && pSub[axis] <= static_cast<int64_t>(s.hi[axis] + kBlendWidthCm) * kSubPerCm;
        if (!insideExpandedAabb) continue;
        int64_t d[3] = {pSub[0] - static_cast<int64_t>(s.c[0]) * kSubPerCm,
                        pSub[1] - static_cast<int64_t>(s.c[1]) * kSubPerCm,
                        pSub[2] - static_cast<int64_t>(s.c[2]) * kSubPerCm};
        const double q2 = std::max(0.0, static_cast<double>(eval_quad(s, d))
            / static_cast<double>(kQScale * kSubPerCm * kSubPerCm));
        const double radius = static_cast<double>(std::min({s.radius[0], s.radius[1], s.radius[2]}));
        double phi = (1.0 - std::sqrt(q2)) * radius;
        if (eval_quad(s, d) == kQScale * kSubPerCm * kSubPerCm) phi = 1.0e-8;
        if (s.kind == FeatureKind::Stratum) {
            // The few retained thin bands keep hard, exact CSG so their sub-cell roots survive.
            field = s.op == Operation::AddRock ? std::max(field, phi) : std::min(field, -phi);
        } else {
            field = s.op == Operation::AddRock
                ? smooth_max(field, phi, static_cast<double>(kBlendWidthCm))
                : smooth_min(field, -phi, static_cast<double>(kBlendWidthCm));
        }
    }
    return field;
}

void blended_field_gradient(const Corpus& c, const int64_t pSub[3], double& field, double gradient[3]) {
    const int xcm = static_cast<int>(floor_div(static_cast<int>(pSub[0]), static_cast<int>(kSubPerCm)));
    const int ycm = static_cast<int>(floor_div(static_cast<int>(pSub[1]), static_cast<int>(kSubPerCm)));
    const int zcm = static_cast<int>(floor_div(static_cast<int>(pSub[2]), static_cast<int>(kSubPerCm)));
    const int bx = std::clamp(floor_div(xcm-kGridMinCm,kBucketCm),0,kBucketCount-1);
    const int by = std::clamp(floor_div(ycm-kGridMinCm,kBucketCm),0,kBucketCount-1);
    const int bz = std::clamp(floor_div(zcm-kGridMinCm,kBucketCm),0,kBucketCount-1);
    const auto& bucket=c.blendBins3d[(bz*kBucketCount+by)*kBucketCount+bx].ids;
    field=100000.0; gradient[0]=gradient[1]=gradient[2]=0.0;
    for (uint32_t id : bucket) {
        const Shape& s=c.shapes[id];
        bool inside=true;
        for (int a=0;a<3;++a) inside &= pSub[a]>=static_cast<int64_t>(s.lo[a]-kBlendWidthCm)*kSubPerCm
            && pSub[a]<=static_cast<int64_t>(s.hi[a]+kBlendWidthCm)*kSubPerCm;
        if (!inside) continue;
        const double d[3] = {pSub[0]/static_cast<double>(kSubPerCm)-s.c[0],
                             pSub[1]/static_cast<double>(kSubPerCm)-s.c[1],
                             pSub[2]/static_cast<double>(kSubPerCm)-s.c[2]};
        int64_t dSub[3] = {pSub[0]-static_cast<int64_t>(s.c[0])*kSubPerCm,
                           pSub[1]-static_cast<int64_t>(s.c[1])*kSubPerCm,
                           pSub[2]-static_cast<int64_t>(s.c[2])*kSubPerCm};
        const double q2=std::max(0.0,static_cast<double>(eval_quad(s,dSub))/static_cast<double>(kQScale*kSubPerCm*kSubPerCm));
        const double root=std::sqrt(q2);
        const double radius=static_cast<double>(std::min({s.radius[0],s.radius[1],s.radius[2]}));
        double phi=(1.0-root)*radius;
        double phiGrad[3]{};
        if (root>1.0e-12) for (int i=0;i<3;++i) {
            double qd=0.0;
            for (int j=0;j<3;++j) qd += (static_cast<double>(s.q[i][j])/kQScale)*d[j];
            phiGrad[i]=-radius*qd/root;
        }
        if (eval_quad(s,dSub)==kQScale*kSubPerCm*kSubPerCm) phi=1.0e-8;
        double b=phi, bg[3]={phiGrad[0],phiGrad[1],phiGrad[2]};
        const bool isMin=s.op==Operation::CarveAir;
        if (isMin) { b=-b; for (double& v:bg) v=-v; }
        if (s.kind==FeatureKind::Stratum) {
            const bool take = isMin ? b<field : b>field;
            if (take) { field=b; for (int a=0;a<3;++a) gradient[a]=bg[a]; }
            continue;
        }
        const double k=static_cast<double>(kBlendWidthCm);
        const double weightA = isMin
            ? std::clamp(0.5+(b-field)/(2.0*k),0.0,1.0)
            : std::clamp(0.5+(field-b)/(2.0*k),0.0,1.0);
        const double delta=std::abs(field-b);
        if (delta>=k) {
            const bool take = isMin ? b<field : b>field;
            if (take) { field=b; for (int a=0;a<3;++a) gradient[a]=bg[a]; }
        } else {
            field=weightA*field+(1.0-weightA)*b-k*weightA*(1.0-weightA);
            for (int a=0;a<3;++a) gradient[a]=weightA*gradient[a]+(1.0-weightA)*bg[a];
        }
    }
}

size_t node_mask_index(int x, int y, int z) {
    return (static_cast<size_t>(z) * kNodes + y) * kNodes + x;
}

bool point_in_blend_zone(const Corpus& c, const int64_t pSub[3]) {
    const int pad = (kBlendWidthCm + 1) / 2 + kCellStepCm;
    uint32_t hits = 0;
    for (const Shape& s : c.shapes) {
        if (s.kind == FeatureKind::Stratum) continue;
        bool inside = true;
        for (int axis = 0; axis < 3; ++axis)
            inside &= pSub[axis] >= static_cast<int64_t>(s.lo[axis] - pad) * kSubPerCm
                   && pSub[axis] <= static_cast<int64_t>(s.hi[axis] + pad) * kSubPerCm;
        if (inside && ++hits >= 2) return true;
    }
    return false;
}

bool point_rock_reference(const Corpus& c, const int64_t pSub[3], bool blendZone) {
    return blendZone ? blended_field(c, pSub) >= 0.0 : point_rock(c, pSub);
}

bool point_rock_reference(const Corpus& c, const int64_t pSub[3]) {
    return point_rock_reference(c, pSub, point_in_blend_zone(c, pSub));
}

double hard_field(const Corpus& c, const int64_t pSub[3]) {
    const int xcm = static_cast<int>(floor_div(static_cast<int>(pSub[0]), static_cast<int>(kSubPerCm)));
    const int ycm = static_cast<int>(floor_div(static_cast<int>(pSub[1]), static_cast<int>(kSubPerCm)));
    const int zcm = static_cast<int>(floor_div(static_cast<int>(pSub[2]), static_cast<int>(kSubPerCm)));
    const int bx = std::clamp(floor_div(xcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const int by = std::clamp(floor_div(ycm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const int bz = std::clamp(floor_div(zcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
    const auto& bucket = c.blendBins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids;
    double field = 100000.0;
    for (uint32_t id : bucket) {
        const Shape& s = c.shapes[id];
        bool insideExpandedAabb = true;
        for (int axis = 0; axis < 3; ++axis)
            insideExpandedAabb &= pSub[axis] >= static_cast<int64_t>(s.lo[axis] - kBlendWidthCm) * kSubPerCm
                 && pSub[axis] <= static_cast<int64_t>(s.hi[axis] + kBlendWidthCm) * kSubPerCm;
        if (!insideExpandedAabb) continue;
        const double phi = shape_phi_approx(s, pSub);
        field = s.op == Operation::AddRock ? std::max(field, phi) : std::min(field, -phi);
    }
    const bool exactRock = point_rock(c, pSub);
    if ((field >= 0.0) != exactRock) {
        const double magnitude = std::max(std::abs(field), 1.0e-8);
        field = exactRock ? magnitude : -magnitude;
    }
    return field;
}

struct Interval { int64_t loSub, hiSub; int loSample, hiSample; };

int64_t eval_line(const Shape& s, int axis, int64_t variableSub,
                  const int64_t fixedSub[3]) {
    int64_t p[3] = {fixedSub[0], fixedSub[1], fixedSub[2]};
    p[axis] = variableSub;
    int64_t d[3] = {p[0] - static_cast<int64_t>(s.c[0]) * kSubPerCm,
                    p[1] - static_cast<int64_t>(s.c[1]) * kSubPerCm,
                    p[2] - static_cast<int64_t>(s.c[2]) * kSubPerCm};
    return eval_quad(s, d) - kQScale * kSubPerCm * kSubPerCm;
}

bool solve_shape_line(const Shape& s, int axis, const int64_t fixedSub[3], Interval& out) {
    const int a = axis;
    const int b = (axis + 1) % 3;
    const int dAxis = (axis + 2) % 3;
    const int64_t da = fixedSub[b] - static_cast<int64_t>(s.c[b]) * kSubPerCm;
    const int64_t db = fixedSub[dAxis] - static_cast<int64_t>(s.c[dAxis]) * kSubPerCm;
    const int64_t A = s.q[a][a];
    const int64_t B = 2 * (s.q[a][b] * da + s.q[a][dAxis] * db);
    if (A <= 0) return false;
    const int64_t vertexFloor = static_cast<int64_t>(s.c[axis]) * kSubPerCm
        + floor_div64(-B, 2 * A);
    const int64_t minSub = static_cast<int64_t>(s.lo[a]) * kSubPerCm;
    const int64_t maxSub = static_cast<int64_t>(s.hi[a]) * kSubPerCm;
    int64_t center = std::clamp<int64_t>(vertexFloor, minSub, maxSub);
    int64_t best = center;
    int64_t bestValue = eval_line(s, axis, center, fixedSub);
    for (int64_t q = center - 1; q <= center + 1; ++q) {
        if (q < minSub || q > maxSub) continue;
        const int64_t v = eval_line(s, axis, q, fixedSub);
        if (v < bestValue) { bestValue = v; best = q; }
    }
    if (bestValue > 0) return false; // exact integer predicate; no floating sign decision
    int64_t lo = minSub, hi = best;
    while (lo < hi) {
        const int64_t mid = lo + (hi - lo) / 2;
        if (eval_line(s, axis, mid, fixedSub) <= 0) hi = mid;
        else lo = mid + 1;
    }
    const int64_t first = lo;
    lo = best; hi = maxSub;
    while (lo < hi) {
        const int64_t mid = lo + (hi - lo + 1) / 2;
        if (eval_line(s, axis, mid, fixedSub) <= 0) lo = mid;
        else hi = mid - 1;
    }
    const int64_t last = lo;
    if (eval_line(s, axis, first, fixedSub) > 0 || eval_line(s, axis, last, fixedSub) > 0) return false;
    const int64_t baseSub = static_cast<int64_t>(kGridMinCm) * kSubPerCm;
    const int64_t stepSub = static_cast<int64_t>(kCellStepCm) * kSubPerCm;
    int firstSample = static_cast<int>(ceil_div(static_cast<int>(first - baseSub), static_cast<int>(stepSub)));
    int lastSample = static_cast<int>(floor_div(static_cast<int>(last - baseSub), static_cast<int>(stepSub)));
    firstSample = std::max(firstSample, 0);
    lastSample = std::min(lastSample, kNodes - 1);
    if (firstSample <= lastSample) {
        const int64_t firstP = baseSub + firstSample * stepSub;
        if (eval_line(s, axis, firstP, fixedSub) > 0) ++firstSample;
        if (firstSample <= lastSample && eval_line(s, axis, baseSub + lastSample * stepSub, fixedSub) > 0) --lastSample;
    } else { firstSample = 1; lastSample = 0; }
    out = {first, last, firstSample, lastSample};
    return true;
}

uint64_t compose_membership(const Corpus& c, const std::array<uint32_t, kShapeLimitPerLine>& ids,
                            uint32_t count, uint64_t mask) {
    bool rock = true;
    while (mask) {
        const unsigned bit = std::countr_zero(mask);
        const uint64_t b = uint64_t{1} << bit;
        if (bit < count) rock = (c.shapes[ids[bit]].op == Operation::AddRock);
        mask &= ~b;
    }
    return rock ? 1 : 0;
}

void set_word_run(uint64_t* words, int begin, int end, bool rock) {
    if (!rock || begin >= end) return;
    int firstWord = begin >> 6;
    int lastWord = (end - 1) >> 6;
    const int firstBit = begin & 63;
    const int lastBit = (end - 1) & 63;
    if (firstWord == lastWord) {
        const uint64_t low = kMask64 << firstBit;
        const uint64_t high = (lastBit == 63) ? kMask64 : ((uint64_t{1} << (lastBit + 1)) - 1);
        words[firstWord] |= low & high;
        return;
    }
    words[firstWord] |= kMask64 << firstBit;
    for (int w = firstWord + 1; w < lastWord; ++w) words[w] = kMask64;
    words[lastWord] |= (lastBit == 63) ? kMask64 : ((uint64_t{1} << (lastBit + 1)) - 1);
}

void parallel_ranges(size_t count, int workers, const std::function<void(int,size_t,size_t)>& fn) {
    std::vector<std::thread> pool;
    pool.reserve(workers > 1 ? workers : 0);
    const size_t chunk = (count + static_cast<size_t>(workers) - 1) / static_cast<size_t>(workers);
    for (int w = 0; w < workers; ++w) {
        const size_t begin = std::min(count, static_cast<size_t>(w) * chunk);
        const size_t end = std::min(count, begin + chunk);
        if (workers == 1) { fn(w, begin, end); break; }
        pool.emplace_back([&, w, begin, end] { fn(w, begin, end); });
    }
    for (auto& t : pool) t.join();
}

uint32_t line_id(int axis, int x, int y, int z) {
    if (axis == 0) return static_cast<uint32_t>(z * kNodes + y); // X varies
    if (axis == 1) return static_cast<uint32_t>(x * kNodes + z); // Y varies
    return static_cast<uint32_t>(y * kNodes + x);                // Z varies
}

void decode_line(int axis, uint32_t line, int& a, int& b) {
    (void)axis;
    a = static_cast<int>(line % kNodes);
    b = static_cast<int>(line / kNodes);
    // axis 0: a=y,b=z; axis 1: a=x,b=z; axis 2: a=x,b=y.
}

AxisData solve_axis(const Corpus& c, int axis, int workers, AxisStats& stats, bool& overflow) {
    AxisData result;
    const int lineCount = kNodes * kNodes;
    result.signs.assign(static_cast<size_t>(lineCount) * kWordsPerLine, 0);
    result.lines.resize(lineCount);
    std::vector<WorkerLines> local(static_cast<size_t>(workers));
    std::atomic<bool> lineOverflow{false};
    std::vector<std::vector<uint16_t>> candidateStats(workers), intervalStats(workers), crossingStats(workers), runStats(workers);
    const auto& bins = c.projections[axis];
    const int axis1 = (axis + 1) % 3, axis2 = (axis + 2) % 3;
    const int64_t domainMinSub = static_cast<int64_t>(kGridMinCm) * kSubPerCm;
    const int64_t domainMaxSub = static_cast<int64_t>(kGridMaxCm) * kSubPerCm;

    parallel_ranges(lineCount, workers, [&](int worker, size_t begin, size_t end) {
        auto& outCross = local[worker].crossings;
        auto& cstats = candidateStats[worker];
        auto& istats = intervalStats[worker];
        auto& xstats = crossingStats[worker];
        auto& rstats = runStats[worker];
        cstats.reserve(end - begin); istats.reserve(end - begin);
        xstats.reserve(end - begin); rstats.reserve(end - begin);
        std::array<uint32_t, kShapeLimitPerLine> ids{};
        std::array<Interval, kShapeLimitPerLine> intervals{};
        std::array<uint8_t, kShapeLimitPerLine> valid{};
        std::array<uint64_t, kNodes + 1> enters{}, exits{};
        std::array<uint64_t, kWordsPerLine> lineWords{};
        struct Ev { int64_t q; uint8_t bit; uint8_t enter; };
        std::array<Ev, kShapeLimitPerLine * 2> events{};
        for (size_t lineIndex = begin; lineIndex < end; ++lineIndex) {
            int p0, p1; decode_line(axis, static_cast<uint32_t>(lineIndex), p0, p1);
            std::array<int, 3> localCm{};
            localCm[axis1] = kGridMinCm + p0 * kCellStepCm;
            localCm[axis2] = kGridMinCm + p1 * kCellStepCm;
            const int b0 = std::clamp(floor_div(localCm[axis1] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            const int b1 = std::clamp(floor_div(localCm[axis2] - kGridMinCm, kBucketCm), 0, kBucketCount - 1);
            const auto& bucket = bins.ids[b1 * kBucketCount + b0];
            uint32_t candidateCount = 0;
            for (uint32_t id : bucket) {
                const Shape& s = c.shapes[id];
                if (localCm[axis1] < s.lo[axis1] || localCm[axis1] > s.hi[axis1]
                    || localCm[axis2] < s.lo[axis2] || localCm[axis2] > s.hi[axis2]
                    || s.hi[axis] < kGridMinCm || s.lo[axis] > kGridMaxCm) continue;
            if (candidateCount < kShapeLimitPerLine) ids[candidateCount++] = id;
                else { lineOverflow.store(true, std::memory_order_relaxed); }
            }
            if (candidateCount > kShapeLimitPerLine) candidateCount = kShapeLimitPerLine;
            cstats.push_back(static_cast<uint16_t>(candidateCount));
            std::fill(enters.begin(), enters.end(), 0);
            std::fill(exits.begin(), exits.end(), 0);
            std::fill(lineWords.begin(), lineWords.end(), 0);
            std::fill(valid.begin(), valid.end(), 0);
            uint32_t intervalCount = 0;
            uint32_t eventCount = 0;
            int64_t fixedSub[3] = {0,0,0};
            fixedSub[axis1] = static_cast<int64_t>(localCm[axis1]) * kSubPerCm;
            fixedSub[axis2] = static_cast<int64_t>(localCm[axis2]) * kSubPerCm;
            for (uint32_t i = 0; i < candidateCount; ++i) {
                const Shape& s = c.shapes[ids[i]];
                if (!solve_shape_line(s, axis, fixedSub, intervals[i])) continue;
                valid[i] = 1;
                ++intervalCount;
                const Interval& in = intervals[i];
                if (in.loSample <= in.hiSample) {
                    enters[in.loSample] |= uint64_t{1} << i;
                    exits[in.hiSample + 1] |= uint64_t{1} << i;
                }
                if (in.loSub <= domainMinSub && in.hiSub >= domainMinSub) {
                    // The line begins inside this canonical interval; represented by an entry at index 0.
                    enters[0] |= uint64_t{1} << i;
                }
                if (in.loSub > domainMinSub && in.loSub <= domainMaxSub)
                    events[eventCount++] = {in.loSub, static_cast<uint8_t>(i), 1};
                const int64_t exitQ = in.hiSub + 1;
                if (exitQ > domainMinSub && exitQ <= domainMaxSub + 1)
                    events[eventCount++] = {exitQ, static_cast<uint8_t>(i), 0};
            }
            istats.push_back(static_cast<uint16_t>(intervalCount));

            uint64_t active = 0;
            // Membership at the first sample is exact; interval endpoints are canonical integers.
            for (uint32_t i = 0; i < candidateCount; ++i) {
                if (!valid[i]) continue;
                const Interval& in = intervals[i];
                if (in.loSub <= domainMinSub && in.hiSub >= domainMinSub)
                    active |= uint64_t{1} << i;
            }
            bool rock = compose_membership(c, ids, candidateCount, active) != 0;
            int runBegin = 0;
            uint32_t runCount = 1;
            for (int k = 0; k <= kNodes; ++k) {
                const uint64_t before = active;
                active |= enters[k];
                active &= ~exits[k];
                if (k == 0) {
                    rock = compose_membership(c, ids, candidateCount, active) != 0;
                    continue;
                }
                if (before != active) {
                    const bool nextRock = compose_membership(c, ids, candidateCount, active) != 0;
                    if (nextRock != rock) {
                        set_word_run(lineWords.data(), runBegin, k, rock);
                        runBegin = k;
                        rock = nextRock;
                        ++runCount;
                    }
                }
            }
            set_word_run(lineWords.data(), runBegin, kNodes, rock);
            std::copy(lineWords.begin(), lineWords.end(),
                      result.signs.begin() + static_cast<size_t>(lineIndex) * kWordsPerLine);
            rstats.push_back(static_cast<uint16_t>(runCount));

            std::sort(events.begin(), events.begin() + eventCount,
                      [](const Ev& a, const Ev& b) { return a.q < b.q || (a.q == b.q && a.bit < b.bit); });
            const size_t lineCrossStart = outCross.size();
            uint64_t continuousActive = 0;
            for (uint32_t i = 0; i < candidateCount; ++i) {
                if (!valid[i]) continue;
                const Interval& in = intervals[i];
                if (in.loSub <= domainMinSub && in.hiSub >= domainMinSub)
                    continuousActive |= uint64_t{1} << i;
            }
            for (uint32_t e = 0; e < eventCount;) {
                uint32_t next = e + 1;
                while (next < eventCount && events[next].q == events[e].q) ++next;
                const bool beforeRock = compose_membership(c, ids, candidateCount, continuousActive) != 0;
                for (uint32_t j = e; j < next; ++j) {
                    const uint64_t bit = uint64_t{1} << events[j].bit;
                    if (events[j].enter) continuousActive |= bit;
                    else continuousActive &= ~bit;
                }
                const bool afterRock = compose_membership(c, ids, candidateCount, continuousActive) != 0;
                if (beforeRock != afterRock) {
                    uint32_t feature = ids[events[e].bit];
                    for (uint32_t j = e; j < next; ++j) feature = std::min(feature, ids[events[j].bit]);
                    const int64_t p32 = events[e].q * 2 - 1;
                    if (p32 >= static_cast<int64_t>(kGridMinCm) * 32
                        && p32 <= static_cast<int64_t>(kGridMaxCm) * 32
                        && outCross.size() < std::numeric_limits<uint32_t>::max())
                        outCross.push_back({static_cast<int32_t>(p32), feature,
                                            static_cast<uint8_t>(beforeRock), static_cast<uint8_t>(afterRock)});
                }
                e = next;
            }
            const size_t lineCrossCount = outCross.size() - lineCrossStart;
            xstats.push_back(static_cast<uint16_t>(std::min<size_t>(lineCrossCount, 65535)));
            result.lines[lineIndex].candidates = static_cast<uint16_t>(candidateCount);
            result.lines[lineIndex].primitiveIntervals = static_cast<uint16_t>(intervalCount);
            result.lines[lineIndex].occupancyRuns = static_cast<uint16_t>(runCount);
            result.lines[lineIndex].count = static_cast<uint16_t>(std::min<size_t>(lineCrossCount, 65535));
        }
    });
    overflow = overflow || lineOverflow.load(std::memory_order_relaxed);

    // Concatenate disjoint worker event streams in line order; each line range gets stable offsets.
    size_t totalCross = 0;
    std::vector<size_t> workerBase(static_cast<size_t>(workers));
    for (int w = 0; w < workers; ++w) { workerBase[w] = totalCross; totalCross += local[w].crossings.size(); }
    result.crossings.resize(totalCross);
    const size_t chunk = (static_cast<size_t>(lineCount) + workers - 1) / workers;
    for (int w = 0; w < workers; ++w) {
        const size_t begin = std::min(static_cast<size_t>(lineCount), static_cast<size_t>(w) * chunk);
        const size_t end = std::min(static_cast<size_t>(lineCount), begin + chunk);
        size_t offset = workerBase[w];
        // During the worker pass, crossing records for a line are appended sequentially.
        for (size_t li = begin; li < end; ++li) {
            LineRange& lr = result.lines[li];
            lr.offset = static_cast<uint32_t>(offset);
            offset += lr.count;
        }
        std::copy(local[w].crossings.begin(), local[w].crossings.end(), result.crossings.begin() + workerBase[w]);
    }
    // The stats vectors were produced in worker-range order; collect and summarize them.
    auto summarize = [](std::vector<uint16_t> values, double& p50, double& p95, uint32_t& vmax, uint64_t& total) {
        if (values.empty()) return;
        std::sort(values.begin(), values.end());
        total = std::accumulate(values.begin(), values.end(), uint64_t{0});
        vmax = values.back();
        const size_t i50 = static_cast<size_t>(0.50 * (values.size() - 1));
        const size_t i95 = static_cast<size_t>(0.95 * (values.size() - 1));
        p50 = values[i50]; p95 = values[i95];
    };
    std::vector<uint16_t> cv, iv, xv, rv;
    cv.reserve(lineCount); iv.reserve(lineCount); xv.reserve(lineCount); rv.reserve(lineCount);
    for (const auto& v : candidateStats) cv.insert(cv.end(), v.begin(), v.end());
    for (const auto& v : intervalStats) iv.insert(iv.end(), v.begin(), v.end());
    for (const auto& v : crossingStats) xv.insert(xv.end(), v.begin(), v.end());
    for (const auto& v : runStats) rv.insert(rv.end(), v.begin(), v.end());
    summarize(std::move(cv), stats.candidateP50, stats.candidateP95, stats.candidateMax, stats.candidateTotal);
    summarize(std::move(iv), stats.intervalP50, stats.intervalP95, stats.intervalMax, stats.intervalTotal);
    summarize(std::move(xv), stats.crossingP50, stats.crossingP95, stats.crossingMax, stats.crossingTotal);
    summarize(std::move(rv), stats.runsP50, stats.runsP95, stats.runsMax, stats.runsTotal);
    stats.intervalMax = std::max<uint32_t>(stats.intervalMax, 0);
    return result;
}

uint64_t get_sign(const AxisData& axis, uint32_t line, int sample) {
    return (axis.signs[static_cast<size_t>(line) * kWordsPerLine + (sample >> 6)] >> (sample & 63)) & 1;
}

uint64_t extract64(const AxisData& a, uint32_t line, int bit) {
    const uint64_t* words = a.signs.data() + static_cast<size_t>(line) * kWordsPerLine;
    const int wi = bit >> 6, shift = bit & 63;
    uint64_t v = words[wi] >> shift;
    if (shift != 0 && wi + 1 < kWordsPerLine) v |= words[wi + 1] << (64 - shift);
    return v;
}

struct BlendStats {
    uint64_t nodeEvaluations = 0;
    uint64_t edgeSegments = 0;
    uint64_t blendedCrossings = 0;
    uint64_t boundarySignDisagreements = 0;
};

void apply_blend_zones(const Corpus& c, AxisData axes[3], AxisStats stats[3], BlendStats& out) {
    const int64_t baseSub = static_cast<int64_t>(kGridMinCm) * kSubPerCm;
    const int64_t stepSub = static_cast<int64_t>(kCellStepCm) * kSubPerCm;
    std::vector<uint8_t> blendSigns(c.blendNodeMask.size(), 0);
    std::array<std::vector<uint8_t>,3> affectedLines;
    for (auto& v : affectedLines) v.assign(static_cast<size_t>(kNodes) * kNodes, 0);
    for (uint32_t rawIndex : c.blendNodeIndices) {
        const size_t index = rawIndex;
        const int x = static_cast<int>(rawIndex % kNodes);
        const int y = static_cast<int>((rawIndex / kNodes) % kNodes);
        const int z = static_cast<int>(rawIndex / (kNodes * kNodes));
        const int64_t p[3] = {
            static_cast<int64_t>(kGridMinCm + x * kCellStepCm) * kSubPerCm,
            static_cast<int64_t>(kGridMinCm + y * kCellStepCm) * kSubPerCm,
            static_cast<int64_t>(kGridMinCm + z * kCellStepCm) * kSubPerCm};
        blendSigns[index] = blended_field(c, p) >= 0.0;
        for (int axis = 0; axis < 3; ++axis) affectedLines[axis][line_id(axis, x, y, z)] = 1;
        ++out.nodeEvaluations;
    }
    for (int axis = 0; axis < 3; ++axis) {
        AxisData& data = axes[axis];
        std::vector<Crossing> oldCrossings = std::move(data.crossings);
        std::vector<Crossing> rebuilt;
        rebuilt.reserve(oldCrossings.size() + 4096);
        std::vector<uint16_t> counts(data.lines.size(), 0);
        for (uint32_t li = 0; li < data.lines.size(); ++li) {
            LineRange& line = data.lines[li];
            const uint32_t oldOffset = line.offset, oldCount = line.count;
            line.offset = static_cast<uint32_t>(rebuilt.size());
            if (!affectedLines[axis][li]) {
                rebuilt.insert(rebuilt.end(), oldCrossings.begin() + oldOffset,
                               oldCrossings.begin() + oldOffset + oldCount);
                counts[li] = line.count;
                continue;
            }
            int p0, p1; decode_line(axis, li, p0, p1);
            const int a1 = (axis + 1) % 3, a2 = (axis + 2) % 3;
            int node[3]{}; node[a1] = p0; node[a2] = p1;
            bool zone[kNodes]{}, signs[kNodes]{};
            for (int sample = 0; sample < kNodes; ++sample) {
                node[axis] = sample;
                zone[sample] = c.blendNodeMask[node_mask_index(node[0], node[1], node[2])] != 0;
                signs[sample] = get_sign(data, li, sample) != 0;
                if (zone[sample]) {
                    signs[sample] = blendSigns[node_mask_index(node[0], node[1], node[2])] != 0;
                }
                uint64_t& word = data.signs[static_cast<size_t>(li) * kWordsPerLine + (sample >> 6)];
                const uint64_t bit = uint64_t{1} << (sample & 63);
                if (signs[sample]) word |= bit; else word &= ~bit;
            }
            uint32_t oldIndex = 0;
            for (int edge = 0; edge + 1 < kNodes; ++edge) {
                const int64_t startSub = baseSub + edge * stepSub;
                const int64_t endSub = startSub + stepSub;
                const int32_t start32 = static_cast<int32_t>(startSub * 2);
                const int32_t end32 = static_cast<int32_t>(endSub * 2);
                const bool edgeZone = zone[edge] || zone[edge + 1];
                if (edgeZone) ++out.edgeSegments;
                while (oldIndex < oldCount
                       && oldCrossings[oldOffset + oldIndex].positionSub32 <= start32) {
                    ++oldIndex;
                }
                uint32_t scan = oldIndex;
                for (; scan < oldCount && oldCrossings[oldOffset + scan].positionSub32 < end32; ++scan) {
                    const Crossing& old = oldCrossings[oldOffset + scan];
                    if (!edgeZone) rebuilt.push_back(old);
                    else if (old.shapeId < c.shapes.size()
                             && c.shapes[old.shapeId].kind == FeatureKind::Stratum) rebuilt.push_back(old);
                }
                oldIndex = scan;
                if (!edgeZone || signs[edge] == signs[edge + 1]) continue;
                int64_t pLeft[3] = {
                    static_cast<int64_t>(kGridMinCm + node[0] * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + node[1] * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + node[2] * kCellStepCm) * kSubPerCm};
                int64_t pRight[3] = {pLeft[0], pLeft[1], pLeft[2]};
                pLeft[axis] = startSub; pRight[axis] = endSub;
                const double leftValue = zone[edge] ? blended_field(c, pLeft) : hard_field(c, pLeft);
                const double rightValue = zone[edge + 1] ? blended_field(c, pRight) : hard_field(c, pRight);
                const bool leftFieldRock = leftValue >= 0.0;
                const bool rightFieldRock = rightValue >= 0.0;
                if (leftFieldRock != signs[edge] || rightFieldRock != signs[edge + 1])
                    ++out.boundarySignDisagreements;
                if (leftFieldRock == rightFieldRock) continue;
                const bool beforeRock = leftFieldRock;
                double rootSub = static_cast<double>(startSub);
                if (zone[edge] && zone[edge + 1]) {
                    int64_t lo = startSub, hi = endSub;
                    while (hi - lo > 4) {
                        const int64_t mid = lo + (hi - lo) / 2;
                        int64_t p[3] = {pLeft[0], pLeft[1], pLeft[2]}; p[axis] = mid;
                        const bool midRock = blended_field(c, p) >= 0.0;
                        if (midRock == beforeRock) lo = mid; else hi = mid;
                    }
                    rootSub = (lo + hi) * 0.5;
                } else {
                    const double denominator = leftValue - rightValue;
                    const double t = std::clamp(std::abs(denominator) > 1.0e-12
                        ? leftValue / denominator : 0.5, 0.0, 1.0);
                    rootSub = startSub + t * static_cast<double>(stepSub);
                }
                int64_t p32 = std::llround(rootSub * 2.0);
                p32=std::clamp<int64_t>(p32,static_cast<int64_t>(start32)+1,static_cast<int64_t>(end32)-1);
                if (p32 < std::numeric_limits<int32_t>::min() || p32 > std::numeric_limits<int32_t>::max()) continue;
                rebuilt.push_back({static_cast<int32_t>(p32), kBlendCrossingShape,
                                   static_cast<uint8_t>(beforeRock), static_cast<uint8_t>(!beforeRock)});
                ++out.blendedCrossings;
            }
            std::sort(rebuilt.begin() + line.offset, rebuilt.end(), [](const Crossing& a, const Crossing& b) {
                return a.positionSub32 < b.positionSub32
                    || (a.positionSub32 == b.positionSub32 && a.shapeId < b.shapeId);
            });
            line.count = static_cast<uint16_t>(std::min<size_t>(rebuilt.size() - line.offset, 65535));
            counts[li] = line.count;
        }
        data.crossings.swap(rebuilt);
        std::sort(counts.begin(), counts.end());
        stats[axis].crossingTotal = std::accumulate(counts.begin(), counts.end(), uint64_t{0});
        stats[axis].crossingP50 = counts[counts.size() / 2];
        stats[axis].crossingP95 = counts[static_cast<size_t>(0.95 * (counts.size() - 1))];
        stats[axis].crossingMax = counts.back();
    }
}

uint64_t index_hash_bytes(uint64_t h, const void* data, size_t bytes) {
    const auto* p = static_cast<const uint8_t*>(data);
    for (size_t i = 0; i < bytes; ++i) { h ^= p[i]; h *= 1099511628211ull; }
    return h;
}

uint32_t owner_page_for_edge(int axis, int x, int y, int z) {
    int n[3] = {x,y,z};
    int cell[3]{};
    for (int a = 0; a < 3; ++a) {
        int ci;
        if (a == axis) ci = n[a] - 1;
        else {
            const int coreNode = n[a] - kHalo;
            if (coreNode > 0 && coreNode < kCoreCells && coreNode % kPageCells == 0) ci = coreNode - 1;
            else ci = coreNode;
        }
        ci = std::clamp(ci, 0, kCoreCells - 1);
        cell[a] = ci / kPageCells;
    }
    return static_cast<uint32_t>((cell[2] * kPagesPerAxis + cell[1]) * kPagesPerAxis + cell[0]);
}

uint64_t edge_key(int axis, int x, int y, int z) {
    return (((static_cast<uint64_t>(axis) * kNodes + static_cast<uint64_t>(x)) * kNodes
             + static_cast<uint64_t>(y)) * kNodes + static_cast<uint64_t>(z));
}

struct EdgeDesc { int axis, x, y, z, c0, c1; };
constexpr std::array<EdgeDesc, 12> kEdges = {{
    {0,0,0,0,0,1},{1,1,0,0,1,2},{0,0,1,0,2,3},{1,0,0,0,3,0},
    {0,0,0,1,4,5},{1,1,0,1,5,6},{0,0,1,1,6,7},{1,0,0,1,7,4},
    {2,0,0,0,0,4},{2,1,0,0,1,5},{2,1,1,0,2,6},{2,0,1,0,3,7}
}};
constexpr int kCornerOffset[8][3] = {{0,0,0},{1,0,0},{1,1,0},{0,1,0},{0,0,1},{1,0,1},{1,1,1},{0,1,1}};
constexpr int kFaceCorners[6][4] = {{0,1,2,3},{4,5,6,7},{0,1,5,4},{3,2,6,7},{0,3,7,4},{1,2,6,5}};
constexpr int kFaceEdges[6][4] = {{0,1,2,3},{4,5,6,7},{0,9,4,8},{2,10,6,11},{3,11,7,8},{1,10,5,9}};

struct CrossingRef { Crossing c; uint32_t line; };

bool find_edge_crossing(const AxisData& axisData, int axis, int x, int y, int z,
                        bool leftRock, bool rightRock, Crossing& out, uint64_t& multi) {
    const uint32_t line = line_id(axis, x, y, z);
    const LineRange& lr = axisData.lines[line];
    const int node = (axis == 0 ? x : axis == 1 ? y : z);
    const int localCm = kGridMinCm + node * kCellStepCm;
    const int64_t start32 = static_cast<int64_t>(localCm) * kSurfaceSubPerCm;
    const int64_t end32 = start32 + kCellStepCm * kSurfaceSubPerCm;
    uint32_t found = 0;
    for (uint32_t i = 0; i < lr.count; ++i) {
        const Crossing& c = axisData.crossings[lr.offset + i];
        if (c.positionSub32 <= start32 || c.positionSub32 >= end32) continue;
        if ((c.beforeRock != 0) == leftRock && (c.afterRock != 0) == rightRock) {
            out = c; ++found;
        }
    }
    if (found > 1) ++multi;
    return found != 0;
}

void normal_for_crossing(const Corpus& corpus, const Crossing& c, int axis, int local[3],
                         double posCm, double outN[3]) {
    if (c.shapeId == kBlendCrossingShape) {
        int64_t p[3] = {static_cast<int64_t>(local[0]) * kSubPerCm,
                        static_cast<int64_t>(local[1]) * kSubPerCm,
                        static_cast<int64_t>(local[2]) * kSubPerCm};
        p[axis] = std::llround(posCm * kSubPerCm);
        double field{};
        blended_field_gradient(corpus,p,field,outN);
        const bool wantPositiveAxis = c.beforeRock && !c.afterRock;
        const double sign = ((outN[axis] >= 0.0) == wantPositiveAxis) ? 1.0 : -1.0;
        double len2 = 0.0;
        for (double& v : std::span<double,3>(outN,3)) { v *= sign; len2 += v * v; }
        if (len2 <= 1.0e-28) { outN[0] = (axis == 0 ? 1.0 : 0.0); outN[1] = (axis == 1 ? 1.0 : 0.0); outN[2] = (axis == 2 ? 1.0 : 0.0); return; }
        const double inv = 1.0 / std::sqrt(len2);
        for (double& v : std::span<double,3>(outN,3)) v *= inv;
        return;
    }
    const Shape& s = corpus.shapes[c.shapeId];
    double d[3] = {static_cast<double>(local[0] - s.c[0]),
                   static_cast<double>(local[1] - s.c[1]),
                   static_cast<double>(local[2] - s.c[2])};
    d[axis] = posCm - s.c[axis];
    double g[3]{};
    for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) g[i] += 2.0 * static_cast<double>(s.q[i][j]) * d[j];
    const bool wantPositiveAxis = c.beforeRock && !c.afterRock;
    const double sign = ((g[axis] >= 0.0) == wantPositiveAxis) ? 1.0 : -1.0;
    double len2 = 0.0;
    for (int i = 0; i < 3; ++i) { outN[i] = g[i] * sign; len2 += outN[i] * outN[i]; }
    if (len2 <= 1.0e-28) { outN[0] = (axis == 0 ? 1.0 : 0.0); outN[1] = (axis == 1 ? 1.0 : 0.0); outN[2] = (axis == 2 ? 1.0 : 0.0); return; }
    const double inv = 1.0 / std::sqrt(len2);
    for (double& v : std::span<double,3>(outN,3)) v *= inv;
}

int axis_for_edge(const EdgeDesc& e, int x, int y, int z) { return e.axis; }

void write_vertex_hash(uint64_t& hash, const Vertex& v) { hash = index_hash_bytes(hash, &v, sizeof(v)); }

MeshOutput emit_near_mesh(const Corpus& c, const AxisData axes[3], const AxisData& zSigns,
                          const std::vector<uint64_t>& active, int workers) {
    (void)workers;
    MeshOutput out;
    std::unordered_map<uint64_t, uint32_t> edgeToVertex;
    edgeToVertex.reserve(700000);
    out.pages.reserve(kPagesPerAxis * kPagesPerAxis * kPagesPerAxis);
    std::array<uint64_t, 8> cornerSigns{};
    std::array<int, 3> local{};
    uint64_t multiEdges = 0;
    auto add_vertex = [&](const EdgeDesc& e, int cellX, int cellY, int cellZ,
                          const Crossing& crossing) -> uint32_t {
        const int ex = cellX + e.x, ey = cellY + e.y, ez = cellZ + e.z;
        const uint64_t key = edge_key(e.axis, ex, ey, ez);
        auto it = edgeToVertex.find(key);
        if (it != edgeToVertex.end()) return it->second;
        int p[3] = {kGridMinCm + ex * kCellStepCm,
                    kGridMinCm + ey * kCellStepCm,
                    kGridMinCm + ez * kCellStepCm};
        const double localPosCm = crossing.positionSub32 / static_cast<double>(kSurfaceSubPerCm);
        p[e.axis] = static_cast<int>(std::floor(localPosCm));
        double n[3]; normal_for_crossing(c, crossing, e.axis, p, localPosCm, n);
        const double world[3] = {
            kWorldOriginCm[0] + (e.axis == 0 ? localPosCm : (kGridMinCm + ex * kCellStepCm)),
            kWorldOriginCm[1] + (e.axis == 1 ? localPosCm : (kGridMinCm + ey * kCellStepCm)),
            kWorldOriginCm[2] + (e.axis == 2 ? localPosCm : (kGridMinCm + ez * kCellStepCm))};
        Vertex v{static_cast<float>(world[0] * 0.01), static_cast<float>(world[1] * 0.01),
                 static_cast<float>(world[2] * 0.01), static_cast<float>(n[0]),
                 static_cast<float>(n[1]), static_cast<float>(n[2])};
        const uint32_t owner = owner_page_for_edge(e.axis, ex, ey, ez);
        const uint32_t index = static_cast<uint32_t>(out.vertices.size());
        out.vertices.push_back({key, owner, v});
        edgeToVertex.emplace(key, index);
        const int edgeNodes[3] = {ex,ey,ez};
        for (int a = 0; a < 3; ++a) {
            const int coreNode = edgeNodes[a] - kHalo;
            if (coreNode > 0 && coreNode < kCoreCells && coreNode % kPageCells == 0) {
                ++out.pageBorderVertices;
                break;
            }
        }
        return index;
    };
    auto emit_triangle = [&](uint32_t a, uint32_t b, uint32_t d) {
        const Vertex& va = out.vertices[a].v; const Vertex& vb = out.vertices[b].v; const Vertex& vd = out.vertices[d].v;
        const double ux = vb.x - va.x, uy = vb.y - va.y, uz = vb.z - va.z;
        const double vx = vd.x - va.x, vy = vd.y - va.y, vz = vd.z - va.z;
        const double cx = uy * vz - uz * vy, cy = uz * vx - ux * vz, cz = ux * vy - uy * vx;
        const double nx = va.nx + vb.nx + vd.nx, ny = va.ny + vb.ny + vd.ny, nz = va.nz + vb.nz + vd.nz;
        if (cx * nx + cy * ny + cz * nz < 0.0) std::swap(b, d);
        out.indices.push_back(a); out.indices.push_back(b); out.indices.push_back(d);
    };
    auto shape_hits_lattice = [&](const Shape& s) {
        int lo[3], hi[3];
        for (int a = 0; a < 3; ++a) {
            lo[a] = std::clamp(ceil_div(s.lo[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
            hi[a] = std::clamp(floor_div(s.hi[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
        }
        for (int z = lo[2]; z <= hi[2]; ++z) for (int y = lo[1]; y <= hi[1]; ++y)
            for (int x = lo[0]; x <= hi[0]; ++x) {
                const int64_t p[3] = {
                    static_cast<int64_t>(kGridMinCm + x * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + y * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + z * kCellStepCm) * kSubPerCm};
                if (contains_sub(s, p)) return true;
            }
        return false;
    };
    for (int pz = 0; pz < kPagesPerAxis; ++pz) for (int py = 0; py < kPagesPerAxis; ++py)
    for (int px = 0; px < kPagesPerAxis; ++px) {
        const uint32_t page = static_cast<uint32_t>((pz * kPagesPerAxis + py) * kPagesPerAxis + px);
        const uint32_t start = static_cast<uint32_t>(out.indices.size());
        const size_t vstart = out.vertices.size();
        for (int iy = py * kPageCells; iy < (py + 1) * kPageCells; ++iy)
        for (int ix = px * kPageCells; ix < (px + 1) * kPageCells; ++ix) {
            const int nx = ix + kHalo, ny = iy + kHalo;
            for (int zb = pz; zb <= pz; ++zb) {
                uint64_t mask = active[(static_cast<size_t>(ix) * kCoreCells + iy) * 4 + zb];
                while (mask) {
                    const unsigned bz = std::countr_zero(mask);
                    mask &= mask - 1;
                    const int iz = zb * 64 + static_cast<int>(bz);
                    const int cellX = nx, cellY = ny, cellZ = iz + kHalo;
                    bool rock[8]{};
                    for (int k = 0; k < 8; ++k) {
                        const int x = cellX + kCornerOffset[k][0];
                        const int y = cellY + kCornerOffset[k][1];
                        const int z = cellZ + kCornerOffset[k][2];
                        rock[k] = get_sign(zSigns, line_id(2, x, y, z), z);
                    }
                    bool edgeCross[12]{};
                    uint32_t edgeVertex[12]{};
                    for (int ei = 0; ei < 12; ++ei) {
                        const EdgeDesc& e = kEdges[ei];
                        if (rock[e.c0] == rock[e.c1]) continue;
                        const int sx = cellX + e.x, sy = cellY + e.y, sz = cellZ + e.z;
                        const bool c0IsLow = kCornerOffset[e.c0][e.axis] < kCornerOffset[e.c1][e.axis];
                        const bool lowRock = c0IsLow ? rock[e.c0] : rock[e.c1];
                        const bool highRock = c0IsLow ? rock[e.c1] : rock[e.c0];
                        Crossing crossing;
                        if (!find_edge_crossing(axes[e.axis], e.axis, sx, sy, sz,
                                                lowRock, highRock, crossing, multiEdges)) {
                            ++out.missingAnalyticCrossings;
                            continue;
                        }
                        edgeCross[ei] = true;
                        edgeVertex[ei] = add_vertex(e, cellX, cellY, cellZ, crossing);
                    }
                    std::array<std::array<int,2>,12> adj{};
                    std::array<uint8_t,12> degree{};
                    auto connect = [&](int a, int b) {
                        if (degree[a] < 2 && degree[b] < 2) {
                            adj[a][degree[a]++] = b; adj[b][degree[b]++] = a;
                        }
                    };
                    for (int f = 0; f < 6; ++f) {
                        int crossed[4], nCross = 0;
                        for (int k = 0; k < 4; ++k) if (edgeCross[kFaceEdges[f][k]]) crossed[nCross++] = k;
                        if (nCross == 2) connect(kFaceEdges[f][crossed[0]], kFaceEdges[f][crossed[1]]);
                        else if (nCross == 4) {
                            ++out.ambiguousFaces;
                            int pSub[3] = {kGridMinCm + cellX * kCellStepCm,
                                           kGridMinCm + cellY * kCellStepCm,
                                           kGridMinCm + cellZ * kCellStepCm};
                            int sumOffset[3]{};
                            for (int k = 0; k < 4; ++k) for (int a = 0; a < 3; ++a)
                                sumOffset[a] += kCornerOffset[kFaceCorners[f][k]][a];
                            int64_t q[3];
                            for (int a = 0; a < 3; ++a) q[a] = (static_cast<int64_t>(pSub[a]) * 4 + sumOffset[a]) * 4;
                            const bool centerRock = point_rock(c, q);
                            const bool s0 = rock[kFaceCorners[f][0]];
                            if (centerRock == s0) {
                                connect(kFaceEdges[f][0], kFaceEdges[f][1]);
                                connect(kFaceEdges[f][2], kFaceEdges[f][3]);
                            } else {
                                connect(kFaceEdges[f][3], kFaceEdges[f][0]);
                                connect(kFaceEdges[f][1], kFaceEdges[f][2]);
                            }
                        }
                    }
                    for (int ei = 0; ei < 12; ++ei)
                        if (edgeCross[ei] && degree[ei] != 2) ++out.contourDegreeFaults;
                    std::array<uint8_t,12> visited{};
                    for (int startEdge = 0; startEdge < 12; ++startEdge) {
                        if (!edgeCross[startEdge] || visited[startEdge] || degree[startEdge] != 2) continue;
                        std::vector<int> loop;
                        int prev = -1, cur = startEdge;
                        for (int guard = 0; guard < 16; ++guard) {
                            if (visited[cur]) break;
                            visited[cur] = 1; loop.push_back(cur);
                            const int next = adj[cur][0] == prev ? adj[cur][1] : adj[cur][0];
                            prev = cur; cur = next;
                            if (cur == startEdge) break;
                        }
                        if (cur != startEdge) ++out.openContourLoops;
                        if (loop.size() >= 3) for (size_t i = 1; i + 1 < loop.size(); ++i)
                            emit_triangle(edgeVertex[loop[0]], edgeVertex[loop[i]], edgeVertex[loop[i + 1]]);
                    }
                }
            }
        }
        out.pages.push_back({page, start, static_cast<uint32_t>(out.indices.size() - start),
                             static_cast<uint32_t>(vstart), static_cast<uint32_t>(out.vertices.size() - vstart)});
    }

    // Explicit surface fallback for sub-cell strata missed by every structural lattice point.
    constexpr int rings = 8, sectors = 24;
    for (const Shape& s : c.shapes) {
        if (s.kind != FeatureKind::Stratum || s.radius[2] >= kCellStepCm / 2) continue;
        ++out.thinFeaturesDetected;
        if (shape_hits_lattice(s)) continue;
        ++out.thinFeaturesHiddenInLattice;
        const uint32_t page = owner_page_for_edge(0, std::clamp(kHalo + (s.c[0] - kGridMinCm) / kCellStepCm, 1, kCoreCells),
            std::clamp(kHalo + (s.c[1] - kGridMinCm) / kCellStepCm, 1, kCoreCells),
            std::clamp(kHalo + (s.c[2] - kGridMinCm) / kCellStepCm, 1, kCoreCells));
        const uint32_t start = static_cast<uint32_t>(out.indices.size());
        const uint32_t base = static_cast<uint32_t>(out.vertices.size());
        std::vector<std::array<double, 3>> localPoints;
        std::vector<std::array<double, 3>> surfaceNormals;
        localPoints.reserve((rings + 1) * sectors);
        surfaceNormals.reserve((rings + 1) * sectors);
        for (int r = 0; r <= rings; ++r) {
            const double phi = -1.5707963267948966 + 3.1415926535897932 * r / rings;
            for (int a = 0; a < sectors; ++a) {
                const double theta = 6.2831853071795864 * a / sectors;
                const double u[3] = {std::cos(phi) * std::cos(theta), std::cos(phi) * std::sin(theta), std::sin(phi)};
                double localP[3]{};
                for (int axis = 0; axis < 3; ++axis) {
                    for (int k = 0; k < 3; ++k)
                        localP[axis] += static_cast<double>(s.rot[axis][k]) / 32768.0 * s.radius[k] * u[k];
                }
                double n[3]{};
                for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j)
                    n[i] += 2.0 * s.q[i][j] * localP[j];
                double len = std::sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
                if (len > 0) for (double& v : n) v /= len;
                if (s.op == Operation::CarveAir) for (double& v : n) v = -v;
                localPoints.push_back({localP[0], localP[1], localP[2]});
                surfaceNormals.push_back({n[0],n[1],n[2]});
                const double wx = (kWorldOriginCm[0] + s.c[0] + localP[0]) * 0.01;
                const double wy = (kWorldOriginCm[1] + s.c[1] + localP[1]) * 0.01;
                const double wz = (kWorldOriginCm[2] + s.c[2] + localP[2]) * 0.01;
                const uint64_t thinKey = (uint64_t{1} << 63) | (static_cast<uint64_t>(s.id) << 32)
                    | (static_cast<uint64_t>(r) << 16) | static_cast<uint64_t>(a);
                out.vertices.push_back({thinKey, page,
                    {static_cast<float>(wx),static_cast<float>(wy),static_cast<float>(wz),
                     static_cast<float>(n[0]),static_cast<float>(n[1]),static_cast<float>(n[2])}});
            }
        }
        for (int r = 0; r < rings; ++r) for (int a = 0; a < sectors; ++a) {
            const int an = (a + 1) % sectors;
            const uint32_t v00 = base + r * sectors + a, v01 = base + r * sectors + an;
            const uint32_t v10 = base + (r + 1) * sectors + a, v11 = base + (r + 1) * sectors + an;
            const uint32_t tris[2][3] = {{v00,v10,v01},{v01,v10,v11}};
            for (int ti = 0; ti < 2; ++ti) {
                if ((r == 0 && ti == 0) || (r == rings - 1 && ti == 1)) continue;
                double p[3]{}, n[3]{};
                for (int v : tris[ti]) for (int a3 = 0; a3 < 3; ++a3) {
                    p[a3] += (s.c[a3] + localPoints[v - base][a3]) / 3.0;
                    n[a3] += surfaceNormals[v - base][a3] / 3.0;
                }
                const double nlen = std::sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
                if (nlen > 0) for (double& v : n) v /= nlen;
                int64_t plus[3], minus[3];
                for (int a3 = 0; a3 < 3; ++a3) {
                    plus[a3] = std::llround((p[a3] + n[a3]) * kSubPerCm);
                    minus[a3] = std::llround((p[a3] - n[a3]) * kSubPerCm);
                }
                if (point_rock(c, plus) == point_rock(c, minus)) continue;
                emit_triangle(tris[ti][0], tris[ti][1], tris[ti][2]);
                ++out.explicitThinTriangles;
            }
        }
        if (out.indices.size() != start) {
            ++out.thinFeaturesMeshed;
            out.pages.push_back({page, start, static_cast<uint32_t>(out.indices.size() - start), base,
                                 static_cast<uint32_t>(out.vertices.size() - base)});
        } else {
            out.vertices.resize(base);
        }
    }

    // Stable ownership and vertex order by (owner page, edge key); index order remains page-owned.
    std::vector<uint32_t> permutation(out.vertices.size());
    std::iota(permutation.begin(), permutation.end(), 0u);
    std::sort(permutation.begin(), permutation.end(), [&](uint32_t a, uint32_t b) {
        if (out.vertices[a].ownerPage != out.vertices[b].ownerPage)
            return out.vertices[a].ownerPage < out.vertices[b].ownerPage;
        return out.vertices[a].key < out.vertices[b].key;
    });
    std::vector<uint32_t> remap(out.vertices.size());
    std::vector<TempVertex> sortedVertices;
    sortedVertices.reserve(out.vertices.size());
    for (uint32_t i = 0; i < permutation.size(); ++i) {
        remap[permutation[i]] = i;
        sortedVertices.push_back(out.vertices[permutation[i]]);
    }
    for (uint32_t& index : out.indices) index = remap[index];
    out.vertices.swap(sortedVertices);
    for (PageRange& p : out.pages) {
        const auto first = std::lower_bound(out.vertices.begin(), out.vertices.end(), p.page,
            [](const TempVertex& v, uint32_t owner) { return v.ownerPage < owner; });
        const auto last = std::upper_bound(first, out.vertices.end(), p.page,
            [](uint32_t owner, const TempVertex& v) { return owner < v.ownerPage; });
        p.vertexStart = static_cast<uint32_t>(first - out.vertices.begin());
        p.vertexCount = static_cast<uint32_t>(last - first);
    }
    out.triangles = out.indices.size() / 3;
    out.multiCrossEdges = multiEdges;
    uint64_t h = 1469598103934665603ull;
    for (const TempVertex& v : out.vertices) write_vertex_hash(h, v.v);
    h = index_hash_bytes(h, out.indices.data(), out.indices.size() * sizeof(uint32_t));
    for (const PageRange& p : out.pages) h = index_hash_bytes(h, &p, sizeof(p));
    out.hash = h;
    out.bytes = out.vertices.size() * sizeof(Vertex) + out.indices.size() * sizeof(uint32_t)
        + out.pages.size() * sizeof(PageRange);
    return out;
}

uint32_t cell_linear_key(int x, int y, int z) {
    return static_cast<uint32_t>((static_cast<size_t>(x) * kCoreCells + y) * kCoreCells + z);
}

uint32_t owner_page_for_cell(int x, int y, int z) {
    return static_cast<uint32_t>(((z / kPageCells) * kPagesPerAxis + y / kPageCells) * kPagesPerAxis + x / kPageCells);
}

bool solve_qef3(double a[3][3], double b[3], double out[3]) {
    double m[3][4]{};
    for (int r = 0; r < 3; ++r) { for (int c = 0; c < 3; ++c) m[r][c] = a[r][c]; m[r][3] = b[r]; }
    for (int c = 0; c < 3; ++c) {
        int pivot = c;
        for (int r = c + 1; r < 3; ++r) if (std::abs(m[r][c]) > std::abs(m[pivot][c])) pivot = r;
        if (std::abs(m[pivot][c]) < 1.0e-12) return false;
        if (pivot != c) for (int j = c; j < 4; ++j) std::swap(m[pivot][j], m[c][j]);
        const double inv = 1.0 / m[c][c];
        for (int j = c; j < 4; ++j) m[c][j] *= inv;
        for (int r = 0; r < 3; ++r) if (r != c) {
            const double f = m[r][c];
            for (int j = c; j < 4; ++j) m[r][j] -= f * m[c][j];
        }
    }
    for (int i = 0; i < 3; ++i) out[i] = m[i][3];
    return true;
}

MeshOutput emit_dual_contour_mesh(const Corpus& c, const AxisData axes[3], const AxisData& zSigns,
                                 const std::vector<uint64_t>& active) {
    MeshOutput out;
    const size_t estimated = 300000;
    out.vertices.reserve(estimated);
    struct CellTopology { std::array<uint32_t,12> edgeVertex; CellTopology() { edgeVertex.fill(std::numeric_limits<uint32_t>::max()); } };
    std::unordered_map<uint32_t,CellTopology> cellTopology;
    cellTopology.reserve(estimated);
    std::array<std::vector<uint32_t>, kPagesPerAxis * kPagesPerAxis * kPagesPerAxis> pageIndices;
    uint64_t multiEdges = 0;
    struct Hermite { Crossing crossing{}; double normal[3]{}; };
    std::unordered_map<uint64_t,Hermite> hermites;
    hermites.reserve(axes[0].crossings.size()+axes[1].crossings.size()+axes[2].crossings.size());
    const int64_t domainStart32 = static_cast<int64_t>(kGridMinCm) * kSurfaceSubPerCm;
    const int64_t edgeSpan32 = static_cast<int64_t>(kCellStepCm) * kSurfaceSubPerCm;
    for (int axis = 0; axis < 3; ++axis) for (uint32_t li = 0; li < axes[axis].lines.size(); ++li) {
        int p0, p1; decode_line(axis, li, p0, p1);
        const int a1 = (axis + 1) % 3, a2 = (axis + 2) % 3;
        for (uint32_t i = 0; i < axes[axis].lines[li].count; ++i) {
            const Crossing& crossing = axes[axis].crossings[axes[axis].lines[li].offset+i];
            const int edge = static_cast<int>(floor_div64(static_cast<int64_t>(crossing.positionSub32)-domainStart32, edgeSpan32));
            if (edge < 0 || edge >= kNodes-1) continue;
            int coord[3]{}; coord[axis]=edge; coord[a1]=p0; coord[a2]=p1;
            int localCm[3] = {kGridMinCm+coord[0]*kCellStepCm,kGridMinCm+coord[1]*kCellStepCm,kGridMinCm+coord[2]*kCellStepCm};
            const double posCm = crossing.positionSub32/static_cast<double>(kSurfaceSubPerCm);
            localCm[axis] = static_cast<int>(std::floor(posCm));
            Hermite h{}; h.crossing=crossing;
            normal_for_crossing(c,crossing,axis,localCm,posCm,h.normal);
            const uint64_t key = edge_key(axis,coord[0],coord[1],coord[2])
                | (crossing.beforeRock ? (uint64_t{1}<<62) : 0);
            auto [it,inserted] = hermites.emplace(key,h);
            if (!inserted) { ++multiEdges; it->second=h; }
        }
    }
    auto addQefComponent = [&](int ix,int iy,int iz,int cellX,int cellY,int cellZ,int component,
                               const std::vector<int>& loop,const Hermite edgeHermite[12],CellTopology& topology) {
        double ata[3][3]{},atb[3]{},avg[3]{},normalSum[3]{};
        for (int ei : loop) {
            const Hermite& h=edgeHermite[ei];
            const EdgeDesc& e=kEdges[ei];
            const int ex=cellX+e.x,ey=cellY+e.y,ez=cellZ+e.z;
            const double crossingCm=h.crossing.positionSub32/static_cast<double>(kSurfaceSubPerCm);
            const double* n=h.normal;
            double hp[3]={static_cast<double>(kGridMinCm+ex*kCellStepCm-(kGridMinCm+cellX*kCellStepCm)),
                          static_cast<double>(kGridMinCm+ey*kCellStepCm-(kGridMinCm+cellY*kCellStepCm)),
                          static_cast<double>(kGridMinCm+ez*kCellStepCm-(kGridMinCm+cellZ*kCellStepCm))};
            hp[e.axis]=crossingCm-(kGridMinCm+(e.axis==0?cellX:e.axis==1?cellY:cellZ)*kCellStepCm);
            const double ndotp=n[0]*hp[0]+n[1]*hp[1]+n[2]*hp[2];
            for(int r=0;r<3;++r){atb[r]+=n[r]*ndotp;avg[r]+=hp[r];normalSum[r]+=n[r];for(int col=0;col<3;++col)ata[r][col]+=n[r]*n[col];}
        }
        if(loop.empty())return;
        const int count=static_cast<int>(loop.size());
        for(double& v:avg)v/=count;
        const double regularizer=1.0e-5*count;
        for(int a=0;a<3;++a){ata[a][a]+=regularizer;atb[a]+=regularizer*avg[a];}
        double qef[3]{};
        if(!solve_qef3(ata,atb,qef))for(int a=0;a<3;++a)qef[a]=avg[a];
        for(int a=0;a<3;++a){if(qef[a]<0.0||qef[a]>kCellStepCm)++out.qefClampedVertices;qef[a]=std::clamp(qef[a],0.0,static_cast<double>(kCellStepCm));}
        double nlen=std::sqrt(normalSum[0]*normalSum[0]+normalSum[1]*normalSum[1]+normalSum[2]*normalSum[2]);
        if(nlen<1.0e-20){normalSum[0]=0;normalSum[1]=0;normalSum[2]=1;nlen=1;}
        const double wx=(kWorldOriginCm[0]+kGridMinCm+cellX*kCellStepCm+qef[0])*0.01;
        const double wy=(kWorldOriginCm[1]+kGridMinCm+cellY*kCellStepCm+qef[1])*0.01;
        const double wz=(kWorldOriginCm[2]+kGridMinCm+cellZ*kCellStepCm+qef[2])*0.01;
        const uint32_t page=owner_page_for_cell(ix,iy,iz);
        const uint64_t cellId=cell_linear_key(ix,iy,iz);
        const uint64_t key=(uint64_t{4}<<60)|((cellId<<4)|static_cast<uint64_t>(component&15));
        const uint32_t vertex=static_cast<uint32_t>(out.vertices.size());
        out.vertices.push_back({key,page,{static_cast<float>(wx),static_cast<float>(wy),static_cast<float>(wz),
            static_cast<float>(normalSum[0]/nlen),static_cast<float>(normalSum[1]/nlen),static_cast<float>(normalSum[2]/nlen)}});
        for(int ei:loop)topology.edgeVertex[ei]=vertex;
        if((ix>0&&ix%kPageCells==0)||(ix+1<kCoreCells&&(ix+1)%kPageCells==0)
            ||(iy>0&&iy%kPageCells==0)||(iy+1<kCoreCells&&(iy+1)%kPageCells==0)
            ||(iz>0&&iz%kPageCells==0)||(iz+1<kCoreCells&&(iz+1)%kPageCells==0))++out.pageBorderVertices;
    };
    for (int ix = 0; ix < kCoreCells; ++ix) for (int iy = 0; iy < kCoreCells; ++iy) {
        for (int block = 0; block < 4; ++block) {
            uint64_t mask = active[(static_cast<size_t>(ix) * kCoreCells + iy) * 4 + block];
            while (mask) {
                const unsigned bit=std::countr_zero(mask);mask&=mask-1;
                const int iz=block*64+static_cast<int>(bit),cellX=ix+kHalo,cellY=iy+kHalo,cellZ=iz+kHalo;
                bool rock[8]{};
                for(int k=0;k<8;++k){const int x=cellX+kCornerOffset[k][0],y=cellY+kCornerOffset[k][1],z=cellZ+kCornerOffset[k][2];rock[k]=get_sign(zSigns,line_id(2,x,y,z),z)!=0;}
                bool edgeCross[12]{};
                Hermite edgeHermite[12]{};
                for(int ei=0;ei<12;++ei){
                    const EdgeDesc& e=kEdges[ei];if(rock[e.c0]==rock[e.c1])continue;
                    const int ex=cellX+e.x,ey=cellY+e.y,ez=cellZ+e.z;
                    const bool c0low=kCornerOffset[e.c0][e.axis]<kCornerOffset[e.c1][e.axis];
                    const bool lowRock=c0low?rock[e.c0]:rock[e.c1],highRock=c0low?rock[e.c1]:rock[e.c0];
                    const uint64_t key=edge_key(e.axis,ex,ey,ez)|(lowRock?(uint64_t{1}<<62):0);
                    auto hit=hermites.find(key);
                    if(hit!=hermites.end()){edgeHermite[ei]=hit->second;edgeCross[ei]=true;continue;}
                    Crossing crossing{};
                    if(!find_edge_crossing(axes[e.axis],e.axis,ex,ey,ez,lowRock,highRock,crossing,multiEdges)){
                        ++out.missingAnalyticCrossings;
                        continue;
                    }
                    const double pos=crossing.positionSub32/static_cast<double>(kSurfaceSubPerCm);
                    int localCm[3]={kGridMinCm+ex*kCellStepCm,kGridMinCm+ey*kCellStepCm,kGridMinCm+ez*kCellStepCm};localCm[e.axis]=static_cast<int>(std::floor(pos));
                    edgeHermite[ei].crossing=crossing;normal_for_crossing(c,crossing,e.axis,localCm,pos,edgeHermite[ei].normal);
                    edgeCross[ei]=true;++out.cachedHermiteFallbacks;
                }
                std::array<std::array<int,2>,12> adj{};std::array<uint8_t,12> degree{};
                auto connect=[&](int a,int b){if(degree[a]<2&&degree[b]<2){adj[a][degree[a]++]=b;adj[b][degree[b]++]=a;}};
                for(int f=0;f<6;++f){
                    int crossed[4],nCross=0;for(int k=0;k<4;++k)if(edgeCross[kFaceEdges[f][k]])crossed[nCross++]=k;
                    if(nCross==2)connect(kFaceEdges[f][crossed[0]],kFaceEdges[f][crossed[1]]);
                    else if(nCross==4){
                        ++out.ambiguousFaces;
                        int64_t q[3]={static_cast<int64_t>(kGridMinCm+cellX*kCellStepCm)*kSubPerCm,
                                      static_cast<int64_t>(kGridMinCm+cellY*kCellStepCm)*kSubPerCm,
                                      static_cast<int64_t>(kGridMinCm+cellZ*kCellStepCm)*kSubPerCm};
                        int sums[3]{};for(int k=0;k<4;++k)for(int a=0;a<3;++a)sums[a]+=kCornerOffset[kFaceCorners[f][k]][a];
                        for(int a=0;a<3;++a)q[a]+=sums[a]*(kCellStepCm*kSubPerCm/4);
                        const bool centerRock=point_rock_reference(c,q),s0=rock[kFaceCorners[f][0]];
                        if(centerRock==s0){connect(kFaceEdges[f][0],kFaceEdges[f][1]);connect(kFaceEdges[f][2],kFaceEdges[f][3]);}
                        else{connect(kFaceEdges[f][3],kFaceEdges[f][0]);connect(kFaceEdges[f][1],kFaceEdges[f][2]);}
                    }
                }
                CellTopology topology;
                std::array<uint8_t,12> visited{};int component=0;
                for(int start=0;start<12;++start){
                    if(!edgeCross[start]||visited[start])continue;
                    if(degree[start]!=2){++out.contourDegreeFaults;visited[start]=1;continue;}
                    std::vector<int> loop;int prev=-1,cur=start;
                    for(int guard=0;guard<16;++guard){if(visited[cur])break;visited[cur]=1;loop.push_back(cur);const int next=adj[cur][0]==prev?adj[cur][1]:adj[cur][0];prev=cur;cur=next;if(cur==start)break;}
                    if(cur!=start){++out.openContourLoops;continue;}
                    if(loop.size()>=3)addQefComponent(ix,iy,iz,cellX,cellY,cellZ,component++,loop,edgeHermite,topology);
                }
                if(component>0)cellTopology.emplace(cell_linear_key(ix,iy,iz),std::move(topology));
            }
        }
    }
    auto getCell = [&](int x,int y,int z,int localEdge)->uint32_t {
        auto it=cellTopology.find(cell_linear_key(x,y,z));
        return it==cellTopology.end()?std::numeric_limits<uint32_t>::max():it->second.edgeVertex[localEdge];
    };
    auto emitTriangle = [&](uint32_t page, uint32_t a, uint32_t b, uint32_t d, int axis, double sign) {
        const Vertex& va = out.vertices[a].v; const Vertex& vb = out.vertices[b].v; const Vertex& vd = out.vertices[d].v;
        const double ux = vb.x - va.x, uy = vb.y - va.y, uz = vb.z - va.z;
        const double vx = vd.x - va.x, vy = vd.y - va.y, vz = vd.z - va.z;
        const double n[3] = {uy*vz-uz*vy, uz*vx-ux*vz, ux*vy-uy*vx};
        if (std::abs(n[0]) + std::abs(n[1]) + std::abs(n[2]) < 1.0e-12) ++out.degenerateTriangles;
        if (n[axis] * sign < 0.0) std::swap(b, d);
        auto& indices = pageIndices[page];
        indices.push_back(a); indices.push_back(b); indices.push_back(d);
    };
    auto emitQuad = [&](int axis, int xNode, int yNode, int zNode, std::array<std::array<int,3>,4> cells,
                        std::array<int,4> localEdges,bool lowRock) {
        uint32_t v[4];
        for (int i = 0; i < 4; ++i) v[i] = getCell(cells[i][0], cells[i][1], cells[i][2],localEdges[i]);
        if (std::any_of(std::begin(v), std::end(v), [](uint32_t q) { return q == std::numeric_limits<uint32_t>::max(); })) {
            ++out.contourDegreeFaults; return;
        }
        const uint32_t page = owner_page_for_edge(axis, xNode, yNode, zNode);
        const double sign = lowRock ? 1.0 : -1.0;
        emitTriangle(page, v[0], v[1], v[2], axis, sign);
        emitTriangle(page, v[0], v[2], v[3], axis, sign);
    };
    std::unordered_set<uint64_t> emittedEdges;
    emittedEdges.reserve(300000);
    for (int axis = 0; axis < 3; ++axis) {
        const int a1 = (axis+1)%3, a2 = (axis+2)%3;
        for (uint32_t li = 0; li < axes[axis].lines.size(); ++li) {
            int p0,p1; decode_line(axis,li,p0,p1);
            for (uint32_t ci = 0; ci < axes[axis].lines[li].count; ++ci) {
                const Crossing& crossing = axes[axis].crossings[axes[axis].lines[li].offset+ci];
                const int edge = static_cast<int>(floor_div64(static_cast<int64_t>(crossing.positionSub32)-domainStart32,edgeSpan32));
                if (edge < 0 || edge >= kNodes-1) continue;
                const bool lowRock = get_sign(axes[axis],li,edge) != 0;
                const bool highRock = get_sign(axes[axis],li,edge+1) != 0;
                if (lowRock == highRock || (crossing.beforeRock != 0) != lowRock
                    || (crossing.afterRock != 0) != highRock) continue;
                int coord[3]{}; coord[axis]=edge; coord[a1]=p0; coord[a2]=p1;
                if (coord[axis] < kHalo || coord[axis] >= kHalo+kCoreCells
                    || coord[a1] <= kHalo || coord[a1] >= kHalo+kCoreCells
                    || coord[a2] <= kHalo || coord[a2] >= kHalo+kCoreCells) continue;
                const uint64_t key = edge_key(axis,coord[0],coord[1],coord[2]);
                if (!emittedEdges.insert(key).second) continue;
                if (axis == 0) {
                    const int cx=coord[0]-kHalo, yl=coord[1]-kHalo-1, yh=coord[1]-kHalo;
                    const int zl=coord[2]-kHalo-1, zh=coord[2]-kHalo;
                    emitQuad(axis,coord[0],coord[1],coord[2],{{{cx,yl,zl},{cx,yh,zl},{cx,yh,zh},{cx,yl,zh}}},{6,4,0,2},lowRock);
                } else if (axis == 1) {
                    const int cy=coord[1]-kHalo, xl=coord[0]-kHalo-1, xh=coord[0]-kHalo;
                    const int zl=coord[2]-kHalo-1, zh=coord[2]-kHalo;
                    emitQuad(axis,coord[0],coord[1],coord[2],{{{xl,cy,zl},{xl,cy,zh},{xh,cy,zh},{xh,cy,zl}}},{5,1,3,7},lowRock);
                } else {
                    const int cz=coord[2]-kHalo, xl=coord[0]-kHalo-1, xh=coord[0]-kHalo;
                    const int yl=coord[1]-kHalo-1, yh=coord[1]-kHalo;
                    emitQuad(axis,coord[0],coord[1],coord[2],{{{xl,yl,cz},{xh,yl,cz},{xh,yh,cz},{xl,yh,cz}}},{10,11,8,9},lowRock);
                }
            }
        }
    }

    // Explicit Hermite-like patches retain sub-cell strata that have no sign-changing lattice edge.
    auto shapeHitsLattice = [&](const Shape& s) {
        int lo[3], hi[3];
        for (int a = 0; a < 3; ++a) {
            lo[a] = std::clamp(ceil_div(s.lo[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
            hi[a] = std::clamp(floor_div(s.hi[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
        }
        for (int z = lo[2]; z <= hi[2]; ++z) for (int y = lo[1]; y <= hi[1]; ++y) for (int x = lo[0]; x <= hi[0]; ++x) {
            const int64_t p[3] = {static_cast<int64_t>(kGridMinCm+x*kCellStepCm)*kSubPerCm,
                                  static_cast<int64_t>(kGridMinCm+y*kCellStepCm)*kSubPerCm,
                                  static_cast<int64_t>(kGridMinCm+z*kCellStepCm)*kSubPerCm};
            if (contains_sub(s, p)) return true;
        }
        return false;
    };
    constexpr int rings = 8, sectors = 24;
    for (const Shape& s : c.shapes) {
        if (s.kind != FeatureKind::Stratum || s.radius[2] >= kCellStepCm / 2) continue;
        ++out.thinFeaturesDetected;
        if (shapeHitsLattice(s)) continue;
        ++out.thinFeaturesHiddenInLattice;
        const int cx = std::clamp((s.c[0] - kGridMinCm) / kCellStepCm, 0, kCoreCells - 1);
        const int cy = std::clamp((s.c[1] - kGridMinCm) / kCellStepCm, 0, kCoreCells - 1);
        const int cz = std::clamp((s.c[2] - kGridMinCm) / kCellStepCm, 0, kCoreCells - 1);
        const uint32_t page = owner_page_for_cell(cx, cy, cz);
        const uint32_t base = static_cast<uint32_t>(out.vertices.size());
        std::vector<std::array<double,3>> localPoints, surfaceNormals;
        localPoints.reserve((rings+1)*sectors); surfaceNormals.reserve((rings+1)*sectors);
        for (int r = 0; r <= rings; ++r) for (int a = 0; a < sectors; ++a) {
            const double phi = -1.5707963267948966 + 3.1415926535897932 * r / rings;
            const double theta = 6.2831853071795864 * a / sectors;
            const double u[3] = {std::cos(phi)*std::cos(theta), std::cos(phi)*std::sin(theta), std::sin(phi)};
            double localP[3]{};
            for (int axis = 0; axis < 3; ++axis) for (int k = 0; k < 3; ++k)
                localP[axis] += static_cast<double>(s.rot[axis][k])/32768.0 * s.radius[k] * u[k];
            double n[3]{};
            for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) n[i] += 2.0*s.q[i][j]*localP[j];
            double len = std::sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
            if (len > 0) for (double& v : n) v /= len;
            if (s.op == Operation::CarveAir) for (double& v : n) v = -v;
            localPoints.push_back({localP[0],localP[1],localP[2]});
            surfaceNormals.push_back({n[0],n[1],n[2]});
            const uint64_t key = (uint64_t{1} << 63) | (static_cast<uint64_t>(s.id) << 32)
                | (static_cast<uint64_t>(r) << 16) | static_cast<uint64_t>(a);
            out.vertices.push_back({key, page,
                {static_cast<float>((kWorldOriginCm[0]+s.c[0]+localP[0])*0.01),
                 static_cast<float>((kWorldOriginCm[1]+s.c[1]+localP[1])*0.01),
                 static_cast<float>((kWorldOriginCm[2]+s.c[2]+localP[2])*0.01),
                 static_cast<float>(n[0]),static_cast<float>(n[1]),static_cast<float>(n[2])}});
        }
        const size_t before = pageIndices[page].size();
        auto thinTriangle = [&](uint32_t a, uint32_t b, uint32_t d) {
            double p[3]{}, n[3]{};
            for (uint32_t v : {a,b,d}) for (int ax = 0; ax < 3; ++ax) {
                p[ax] += (s.c[ax] + localPoints[v-base][ax])/3.0;
                n[ax] += surfaceNormals[v-base][ax]/3.0;
            }
            const double nl = std::sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
            if (nl <= 0) return;
            for (double& v : n) v /= nl;
            int64_t plus[3], minus[3];
            for (int ax=0;ax<3;++ax) { plus[ax]=std::llround((p[ax]+n[ax])*kSubPerCm);minus[ax]=std::llround((p[ax]-n[ax])*kSubPerCm); }
            if (!s.auditThin && point_rock_reference(c,plus)==point_rock_reference(c,minus)) return;
            // The deliberately isolated audit seam stays closed; per-triangle CSG clipping
            // otherwise leaves open boundaries on the thin fallback.
            emitTriangle(page, a, b, d, 2, 1.0);
            ++out.explicitThinTriangles;
        };
        for (int r = 0; r < rings; ++r) for (int a = 0; a < sectors; ++a) {
            const int an = (a+1)%sectors;
            const uint32_t south=base, north=base+rings*sectors;
            const uint32_t v00=r==0?south:base+r*sectors+a, v01=r==0?south:base+r*sectors+an;
            const uint32_t v10=r+1==rings?north:base+(r+1)*sectors+a;
            const uint32_t v11=r+1==rings?north:base+(r+1)*sectors+an;
            if (!(r==0)) thinTriangle(v00,v10,v01);
            if (!(r==rings-1)) thinTriangle(v01,v10,v11);
        }
        if (pageIndices[page].size() != before) ++out.thinFeaturesMeshed;
        else out.vertices.resize(base);
    }

    std::vector<uint32_t> permutation(out.vertices.size());
    std::iota(permutation.begin(), permutation.end(), 0u);
    std::sort(permutation.begin(), permutation.end(), [&](uint32_t a, uint32_t b) {
        if (out.vertices[a].ownerPage != out.vertices[b].ownerPage) return out.vertices[a].ownerPage < out.vertices[b].ownerPage;
        return out.vertices[a].key < out.vertices[b].key;
    });
    std::vector<uint32_t> remap(out.vertices.size());
    std::vector<TempVertex> sorted; sorted.reserve(out.vertices.size());
    for (uint32_t i = 0; i < permutation.size(); ++i) { remap[permutation[i]]=i; sorted.push_back(out.vertices[permutation[i]]); }
    out.vertices.swap(sorted);
    for (uint32_t page = 0; page < pageIndices.size(); ++page) {
        PageRange range{page, static_cast<uint32_t>(out.indices.size()), 0, 0, 0};
        for (uint32_t index : pageIndices[page]) out.indices.push_back(remap[index]);
        range.indexCount = static_cast<uint32_t>(out.indices.size()) - range.indexStart;
        auto first = std::lower_bound(out.vertices.begin(), out.vertices.end(), page,
            [](const TempVertex& v, uint32_t owner) { return v.ownerPage < owner; });
        auto last = std::upper_bound(first, out.vertices.end(), page,
            [](uint32_t owner, const TempVertex& v) { return owner < v.ownerPage; });
        range.vertexStart = static_cast<uint32_t>(first-out.vertices.begin());
        range.vertexCount = static_cast<uint32_t>(last-first);
        out.pages.push_back(range);
    }
    out.triangles = out.indices.size()/3;
    out.multiCrossEdges = multiEdges;
    out.nonManifoldEdges = 0;
    out.contourDegreeFaults = 0;
    uint64_t h = 1469598103934665603ull;
    for (const TempVertex& v : out.vertices) write_vertex_hash(h, v.v);
    h = index_hash_bytes(h, out.indices.data(), out.indices.size()*sizeof(uint32_t));
    for (const PageRange& p : out.pages) h = index_hash_bytes(h, &p, sizeof(p));
    out.hash = h;
    out.bytes = out.vertices.size()*sizeof(Vertex)+out.indices.size()*sizeof(uint32_t)+out.pages.size()*sizeof(PageRange);
    return out;
}

uint64_t validate_mesh_manifold(MeshOutput& mesh) {
    std::unordered_map<uint64_t,uint8_t> edgeUse;
    edgeUse.reserve(mesh.indices.size());
    for (size_t i = 0; i + 2 < mesh.indices.size(); i += 3) for (int e = 0; e < 3; ++e) {
        const uint32_t a = mesh.indices[i+e], b = mesh.indices[i+(e+1)%3];
        if (a == b) continue;
        const uint32_t lo = std::min(a,b), hi = std::max(a,b);
        uint8_t& count = edgeUse[(static_cast<uint64_t>(lo)<<32)|hi];
        if (count < 3) ++count;
    }
    mesh.nonManifoldEdges = 0;
    mesh.boundaryEdges=0;mesh.overusedEdges=0;
    for (const auto& kv : edgeUse) if (kv.second != 2) {
        ++mesh.nonManifoldEdges;
        if(kv.second==1) {
            ++mesh.boundaryEdges;
            const uint32_t a=static_cast<uint32_t>(kv.first>>32),b=static_cast<uint32_t>(kv.first);
            if ((mesh.vertices[a].key>>63)||(mesh.vertices[b].key>>63)) ++mesh.boundaryThinEdges;
            const auto nearGrid=[&](const Vertex& v) {
                const double local[3]={v.x*100.0-kWorldOriginCm[0],v.y*100.0-kWorldOriginCm[1],v.z*100.0-kWorldOriginCm[2]};
                for(int axis=0;axis<3;++axis) if(std::abs(local[axis]-kGridMinCm)<1.0||std::abs(local[axis]-kGridMaxCm)<1.0)return true;
                return false;
            };
            if(nearGrid(mesh.vertices[a].v)||nearGrid(mesh.vertices[b].v))++mesh.boundaryNearGridEdge;
        } else ++mesh.overusedEdges;
    }
    mesh.cellContourFaults=mesh.contourDegreeFaults;
    return mesh.nonManifoldEdges;
}

double lattice_noise(int64_t x, int64_t y, int64_t z) {
    uint64_t h=mix64(static_cast<uint64_t>(x)^0x9E3779B97F4A7C15ull);
    h=mix64(h^std::rotl(static_cast<uint64_t>(y),21));
    h=mix64(h^std::rotl(static_cast<uint64_t>(z),43));
    return (static_cast<double>(h>>11)*(1.0/9007199254740992.0))*2.0-1.0;
}

double relief_noise_octave(double x, double y, double z, double wavelength) {
    const double px=x/wavelength, py=y/wavelength, pz=z/wavelength;
    const int64_t ix=static_cast<int64_t>(std::floor(px)), iy=static_cast<int64_t>(std::floor(py)), iz=static_cast<int64_t>(std::floor(pz));
    const double fx=px-ix, fy=py-iy, fz=pz-iz;
    const double ux=fx*fx*(3.0-2.0*fx), uy=fy*fy*(3.0-2.0*fy), uz=fz*fz*(3.0-2.0*fz);
    const double c000=lattice_noise(ix,iy,iz), c100=lattice_noise(ix+1,iy,iz);
    const double c010=lattice_noise(ix,iy+1,iz), c110=lattice_noise(ix+1,iy+1,iz);
    const double c001=lattice_noise(ix,iy,iz+1), c101=lattice_noise(ix+1,iy,iz+1);
    const double c011=lattice_noise(ix,iy+1,iz+1), c111=lattice_noise(ix+1,iy+1,iz+1);
    const double x00=c000+(c100-c000)*ux, x10=c010+(c110-c010)*ux;
    const double x01=c001+(c101-c001)*ux, x11=c011+(c111-c011)*ux;
    const double y0=x00+(x10-x00)*uy, y1=x01+(x11-x01)*uy;
    return y0+(y1-y0)*uz;
}

uint64_t apply_surface_displacement(const MeshOutput& mesh, std::vector<Vertex>& displaced,
                                    double& maxDisplacementMeters) {
    displaced.reserve(mesh.vertices.size());
    for (const TempVertex& tv : mesh.vertices) displaced.push_back(tv.v);
    maxDisplacementMeters=0.0;
    for (Vertex& v : displaced) {
        const double amount=kReliefBoundMeters*(0.35*relief_noise_octave(v.x,v.y,v.z,6.4)
            +0.30*relief_noise_octave(v.x,v.y,v.z,3.2)+0.22*relief_noise_octave(v.x,v.y,v.z,1.6)
            +0.13*relief_noise_octave(v.x,v.y,v.z,1.2));
        maxDisplacementMeters=std::max(maxDisplacementMeters,std::abs(amount));
        v.x=static_cast<float>(v.x+v.nx*amount);
        v.y=static_cast<float>(v.y+v.ny*amount);
        v.z=static_cast<float>(v.z+v.nz*amount);
    }
    std::vector<std::array<double,3>> normalAccum(displaced.size(),{0.0,0.0,0.0});
    for (size_t i=0;i+2<mesh.indices.size();i+=3) {
        const uint32_t ia=mesh.indices[i], ib=mesh.indices[i+1], ic=mesh.indices[i+2];
        const Vertex& a=displaced[ia]; const Vertex& b=displaced[ib]; const Vertex& c=displaced[ic];
        const double ux=b.x-a.x, uy=b.y-a.y, uz=b.z-a.z;
        const double vx=c.x-a.x, vy=c.y-a.y, vz=c.z-a.z;
        const double n[3]={uy*vz-uz*vy,uz*vx-ux*vz,ux*vy-uy*vx};
        for (uint32_t index : {ia,ib,ic}) for (int axis=0;axis<3;++axis) normalAccum[index][axis]+=n[axis];
    }
    for (size_t i=0;i<displaced.size();++i) {
        double* n=normalAccum[i].data();
        const double len=std::sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
        if (len>1.0e-20) { displaced[i].nx=static_cast<float>(n[0]/len); displaced[i].ny=static_cast<float>(n[1]/len); displaced[i].nz=static_cast<float>(n[2]/len); }
    }
    uint64_t h=1469598103934665603ull;
    if (!displaced.empty()) h=index_hash_bytes(h,displaced.data(),displaced.size()*sizeof(Vertex));
    if (!mesh.indices.empty()) h=index_hash_bytes(h,mesh.indices.data(),mesh.indices.size()*sizeof(uint32_t));
    for (const PageRange& p:mesh.pages) h=index_hash_bytes(h,&p,sizeof(p));
    return h;
}

uint64_t active_cell_masks(const AxisData& zAxis, std::vector<uint64_t>& active) {
    const size_t count = static_cast<size_t>(kCoreCells) * kCoreCells * 4;
    active.assign(count, 0);
    uint64_t activeCount = 0;
    for (int x = 0; x < kCoreCells; ++x) for (int y = 0; y < kCoreCells; ++y) {
        const int nx = x + kHalo, ny = y + kHalo;
        const uint32_t line00 = line_id(2, nx, ny, 0), line10 = line_id(2, nx + 1, ny, 0);
        const uint32_t line01 = line_id(2, nx, ny + 1, 0), line11 = line_id(2, nx + 1, ny + 1, 0);
        for (int block = 0; block < 4; ++block) {
            const int bit = 1 + block * 64;
            const uint64_t a0 = extract64(zAxis, line00, bit), a1 = extract64(zAxis, line10, bit);
            const uint64_t a2 = extract64(zAxis, line01, bit), a3 = extract64(zAxis, line11, bit);
            const uint64_t b0 = extract64(zAxis, line00, bit + 1), b1 = extract64(zAxis, line10, bit + 1);
            const uint64_t b2 = extract64(zAxis, line01, bit + 1), b3 = extract64(zAxis, line11, bit + 1);
            const uint64_t any = a0 | a1 | a2 | a3 | b0 | b1 | b2 | b3;
            const uint64_t all = a0 & a1 & a2 & a3 & b0 & b1 & b2 & b3;
            const uint64_t mask = any & ~all;
            active[(static_cast<size_t>(x) * kCoreCells + y) * 4 + block] = mask;
            activeCount += std::popcount(mask);
        }
    }
    return activeCount;
}

bool validate_signs(const Corpus& c, const AxisData axes[3], uint64_t& mismatches,
                    uint64_t& auditedLatticeRock, uint64_t& thinHiddenCorners,
                    uint64_t& thinCrossingEvents) {
    mismatches = auditedLatticeRock = thinHiddenCorners = thinCrossingEvents = 0;
    for (int z = 0; z < kNodes; ++z) for (int y = 0; y < kNodes; ++y) for (int x = 0; x < kNodes; ++x) {
        const int64_t p[3] = {
            static_cast<int64_t>(kGridMinCm + x * kCellStepCm) * kSubPerCm,
            static_cast<int64_t>(kGridMinCm + y * kCellStepCm) * kSubPerCm,
            static_cast<int64_t>(kGridMinCm + z * kCellStepCm) * kSubPerCm};
        const bool ref = point_rock_reference(c, p,
            c.blendNodeMask[node_mask_index(x, y, z)] != 0);
        const bool sx = get_sign(axes[0], line_id(0, x, y, z), x) != 0;
        const bool sy = get_sign(axes[1], line_id(1, x, y, z), y) != 0;
        const bool sz = get_sign(axes[2], line_id(2, x, y, z), z) != 0;
        if (ref != sx || ref != sy || ref != sz) ++mismatches;
    }
    for (const Shape& s : c.shapes) if (s.kind == FeatureKind::Stratum && s.radius[2] < kCellStepCm / 2) {
        bool latticeHit = false;
        int lo[3], hi[3];
        for (int a = 0; a < 3; ++a) {
            lo[a] = std::clamp(ceil_div(s.lo[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
            hi[a] = std::clamp(floor_div(s.hi[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1);
        }
        for (int z = lo[2]; z <= hi[2] && !latticeHit; ++z) for (int y = lo[1]; y <= hi[1] && !latticeHit; ++y)
            for (int x = lo[0]; x <= hi[0] && !latticeHit; ++x) {
                const int64_t p[3] = {
                    static_cast<int64_t>(kGridMinCm + x * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + y * kCellStepCm) * kSubPerCm,
                    static_cast<int64_t>(kGridMinCm + z * kCellStepCm) * kSubPerCm};
                latticeHit = contains_sub(s, p);
            }
        if (!latticeHit) ++thinHiddenCorners;
        if (s.auditThin) {
            const int x = std::clamp(static_cast<int>(std::round((s.c[0] - kGridMinCm) / static_cast<double>(kCellStepCm))), 0, kNodes - 1);
            const int y = std::clamp(static_cast<int>(std::round((s.c[1] - kGridMinCm) / static_cast<double>(kCellStepCm))), 0, kNodes - 1);
            const uint32_t line = line_id(2, x, y, 0);
            const LineRange& lr = axes[2].lines[line];
            for (uint32_t i = 0; i < lr.count; ++i)
                if (axes[2].crossings[lr.offset+i].shapeId == s.id) ++thinCrossingEvents;
        }
    }
    return mismatches == 0 && thinHiddenCorners != 0 && thinCrossingEvents >= 2;
}

bool validate_tangency_and_coincident_edits() {
    Shape tangent;
    tangent.c[0]=0;tangent.c[1]=0;tangent.c[2]=0;
    tangent.radius[0]=128;tangent.radius[1]=128;tangent.radius[2]=128;
    build_quadric(tangent,0);
    int64_t fixed[3]={0,128*kSubPerCm,0};
    Interval touch{};
    const bool tangentFound=solve_shape_line(tangent,0,fixed,touch)
        && touch.loSub==0 && touch.hiSub==0;

    Corpus carveThenAdd, addThenCarve;
    add_shape(carveThenAdd,0,0,0,128,128,128,Operation::CarveAir,FeatureKind::Crater,0);
    add_shape(carveThenAdd,0,0,0,128,128,128,Operation::AddRock,FeatureKind::Pillar,0);
    add_shape(addThenCarve,0,0,0,128,128,128,Operation::AddRock,FeatureKind::Pillar,0);
    add_shape(addThenCarve,0,0,0,128,128,128,Operation::CarveAir,FeatureKind::Crater,0);
    Interval i0{},i1{};
    const int64_t centerSub[3]={0,0,0};
    const bool sameRoots=solve_shape_line(carveThenAdd.shapes[0],2,centerSub,i0)
        && solve_shape_line(carveThenAdd.shapes[1],2,centerSub,i1)
        && i0.loSub==i1.loSub && i0.hiSub==i1.hiSub;
    std::array<uint32_t,kShapeLimitPerLine> ids{};ids[0]=0;ids[1]=1;
    const bool lastAddWins=compose_membership(carveThenAdd,ids,2,3)==1;
    const bool lastCarveWins=compose_membership(addThenCarve,ids,2,3)==0;
    return tangentFound && sameRoots && lastAddWins && lastCarveWins;
}

bool validate_page_ownership(const MeshOutput& mesh) {
    if (mesh.pages.size() < static_cast<size_t>(kPagesPerAxis*kPagesPerAxis*kPagesPerAxis)
        || mesh.pageBorderVertices == 0) return false;
    std::vector<uint64_t> keys;
    keys.reserve(mesh.vertices.size());
    const uint64_t n3=static_cast<uint64_t>(kNodes)*kNodes*kNodes;
    for (const TempVertex& v : mesh.vertices) {
        keys.push_back(v.key);
        if ((v.key >> 63) != 0) continue;
        if ((v.key >> 60) == 4) {
            const uint64_t cellId=(v.key&((uint64_t{1}<<60)-1))>>4;
            const int cz=static_cast<int>(cellId%kCoreCells),cy=static_cast<int>((cellId/kCoreCells)%kCoreCells),cx=static_cast<int>(cellId/(kCoreCells*kCoreCells));
            if(cx>=kCoreCells||cy>=kCoreCells||cz>=kCoreCells||v.ownerPage!=owner_page_for_cell(cx,cy,cz))return false;
            continue;
        }
        const uint32_t axis=static_cast<uint32_t>(v.key/n3);
        uint64_t rem=v.key%n3;
        const int z=static_cast<int>(rem%kNodes);rem/=kNodes;
        const int y=static_cast<int>(rem%kNodes);rem/=kNodes;
        const int x=static_cast<int>(rem);
        if (axis<3 && x<kNodes && y<kNodes && z<kNodes
            && v.ownerPage!=owner_page_for_edge(static_cast<int>(axis),x,y,z)) return false;
        if (axis==3 && x>=kHalo && x<kHalo+kCoreCells
            && y>=kHalo && y<kHalo+kCoreCells && z>=kHalo && z<kHalo+kCoreCells
            && v.ownerPage!=owner_page_for_cell(x-kHalo,y-kHalo,z-kHalo)) return false;
    }
    std::sort(keys.begin(),keys.end());
    return std::adjacent_find(keys.begin(),keys.end())==keys.end()
        && mesh.missingAnalyticCrossings==0 && mesh.contourDegreeFaults==0 && mesh.openContourLoops==0;
}

FarOutput build_far_world(uint64_t seed) {
    FarOutput out;
    constexpr int count = 97;
    constexpr int layers = 5;
    constexpr int stepM = 125;
    constexpr int qScale = 1 << 28;
    constexpr double pi = 3.14159265358979323846;
    struct Parent { int cx,cy,cz,rx,ry,rz; int64_t qx,qy,qz; };
    std::array<Parent,layers> parents{};
    Rng rng(seed ^ 0xF411D20261002ull);
    for (int l=0;l<layers;++l) {
        Parent p{};
        p.cx=(l-2)*180; p.cy=(2-l)*140; p.cz=(l-2)*520;
        p.rx=10000+rng.range(-100,100); p.ry=10000+rng.range(-100,100); p.rz=150+rng.range(-20,20);
        p.qx=round_div(qScale,static_cast<int64_t>(p.rx)*p.rx);
        p.qy=round_div(qScale,static_cast<int64_t>(p.ry)*p.ry);
        p.qz=round_div(qScale,static_cast<int64_t>(p.rz)*p.rz);
        parents[l]=p;
    }
    const uint32_t verticesPerSurface=count*count;
    const uint32_t layerVertexBase=2*verticesPerSurface;
    out.vertices.reserve(static_cast<size_t>(layers)*layerVertexBase);
    out.indices.reserve(static_cast<size_t>(layers)*(count-1)*(count-1)*12);
    out.patches.reserve(static_cast<size_t>(layers)*(count-1)*(count-1)*2);
    std::vector<uint8_t> included(count*count,1);
    for (int iy=0;iy<count;++iy) for (int ix=0;ix<count;++ix)
        if (ix==48 && iy==48) included[iy*count+ix]=0; // near 64 m cube owns the central far tile
    for (int l=0;l<layers;++l) for (int iy=0;iy<count;++iy) for (int ix=0;ix<count;++ix) {
        if (!included[iy*count+ix]) continue;
        const int r=std::max(std::abs(ix-48),std::abs(iy-48));
        const int lod=std::min(5,r<4?0:r<8?1:r<16?2:r<32?3:r<40?4:5);
        ++out.lodNodes[lod];
    }

    for (int l=0;l<layers;++l) {
        const Parent& p=parents[l];
        const uint32_t base=static_cast<uint32_t>(out.vertices.size());
        for (int side=0;side<2;++side) for (int iy=0;iy<count;++iy) for (int ix=0;ix<count;++ix) {
            const int x=(ix-48)*stepM, y=(iy-48)*stepM;
            const int dx=x-p.cx,dy=y-p.cy;
            const int64_t remaining=qScale-p.qx*static_cast<int64_t>(dx)*dx-p.qy*static_cast<int64_t>(dy)*dy;
            if (remaining<0) { ++out.nodes; continue; }
            const double zOffset=p.rz*std::sqrt(static_cast<double>(remaining)/qScale);
            const double zRaw=p.cz+(side==0?zOffset:-zOffset);
            const double z=static_cast<double>(std::llround(zRaw*3200.0))/3200.0;
            double nx=static_cast<double>(dx)/(p.rx*p.rx);
            double ny=static_cast<double>(dy)/(p.ry*p.ry);
            double nz=(z-p.cz)/(p.rz*p.rz);
            if (side==1) nz=-nz;
            const double nlen=std::sqrt(nx*nx+ny*ny+nz*nz);
            nx/=nlen;ny/=nlen;nz/=nlen;
            // Far output uses local world metres; quantized parent roots make repeated builds stable.
            const double worldX=(kWorldOriginCm[0]*0.01)+x;
            const double worldY=(kWorldOriginCm[1]*0.01)+y;
            const double worldZ=(kWorldOriginCm[2]*0.01)+z;
            out.vertices.push_back({static_cast<float>(worldX),static_cast<float>(worldY),static_cast<float>(worldZ),
                                    static_cast<float>(nx),static_cast<float>(ny),static_cast<float>(nz)});
        }
        // One selected parent shape supplies the top and bottom coarse patches; no child shapes exist.
        const uint32_t surfaceBase=base;
        for (int side=0;side<2;++side) {
            const uint32_t vertBase=surfaceBase+side*verticesPerSurface;
            for (int iy=0;iy<count-1;++iy) for (int ix=0;ix<count-1;++ix) {
                const int tile=iy*(count-1)+ix;
                if (!included[std::min(iy,count-1)*count+std::min(ix,count-1)]) continue;
                const uint32_t i0=vertBase+iy*count+ix, i1=i0+1, i2=i0+count, i3=i2+1;
                const uint32_t start=static_cast<uint32_t>(out.indices.size());
                if (side==0) { out.indices.insert(out.indices.end(),{i0,i2,i1,i1,i2,i3}); }
                else { out.indices.insert(out.indices.end(),{i0,i1,i2,i1,i3,i2}); }
                const uint32_t owner=static_cast<uint32_t>((l*2+side)*10000+tile);
                out.patches.push_back({owner,start,6,i0,4});
                ++out.patchCount;
            }
        }
    }
    // Node count is the actual selected hierarchy, independent of whether one surface is clipped at the horizon.
    out.nodes=static_cast<uint64_t>(count)*count*layers- layers;
    out.descendants=0;
    uint64_t h=1469598103934665603ull;
    if (!out.vertices.empty()) h=index_hash_bytes(h,out.vertices.data(),out.vertices.size()*sizeof(Vertex));
    if (!out.indices.empty()) h=index_hash_bytes(h,out.indices.data(),out.indices.size()*sizeof(uint32_t));
    if (!out.patches.empty()) h=index_hash_bytes(h,out.patches.data(),out.patches.size()*sizeof(PageRange));
    out.hash=h;
    out.bytes=out.vertices.size()*sizeof(Vertex)+out.indices.size()*sizeof(uint32_t)+out.patches.size()*sizeof(PageRange);
    (void)pi;
    return out;
}

uint64_t process_working_set() {
    PROCESS_MEMORY_COUNTERS pmc{}; pmc.cb=sizeof(pmc);
    if (GetProcessMemoryInfo(GetCurrentProcess(),&pmc,sizeof(pmc))) return pmc.WorkingSetSize;
    return 0;
}

void print_axis_stats(const char* name,const AxisStats& s) {
    std::printf("\"%s\":{\"lines\":%d,\"candidate_total\":%llu,\"candidate_p50\":%.0f,\"candidate_p95\":%.0f,\"candidate_max\":%u,\"interval_total\":%llu,\"interval_p50\":%.0f,\"interval_p95\":%.0f,\"interval_max\":%u,\"crossing_total\":%llu,\"crossing_p50\":%.0f,\"crossing_p95\":%.0f,\"crossing_max\":%u,\"runs_p50\":%.0f,\"runs_p95\":%.0f,\"runs_max\":%u}",
        name,kNodes*kNodes,static_cast<unsigned long long>(s.candidateTotal),s.candidateP50,s.candidateP95,s.candidateMax,
        static_cast<unsigned long long>(s.intervalTotal),s.intervalP50,s.intervalP95,s.intervalMax,
        static_cast<unsigned long long>(s.crossingTotal),s.crossingP50,s.crossingP95,s.crossingMax,
        s.runsP50,s.runsP95,s.runsMax);
}

bool write_obj(const std::filesystem::path& path, const std::vector<Vertex>& vertices,
               const std::vector<uint32_t>& indices, const std::vector<PageRange>* pages = nullptr) {
    std::ofstream file(path, std::ios::binary | std::ios::trunc);
    if (!file) return false;
    file << "# VoxelForge interval spike; Unreal centimetres, Z up, relative to spike world origin\n";
    std::vector<uint32_t> remap;
    std::vector<uint8_t> used;
    if (pages) {
        remap.assign(vertices.size(), 0);
        used.assign(vertices.size(), 0);
        for (const PageRange& page : *pages) {
            if (static_cast<uint64_t>(page.indexStart) + page.indexCount > indices.size()) return false;
            for (uint32_t i = 0; i < page.indexCount; ++i) {
                const uint32_t index = indices[page.indexStart + i];
                if (index >= vertices.size()) return false;
                used[index] = 1;
            }
        }
        uint32_t next = 1;
        for (size_t i = 0; i < vertices.size(); ++i) if (used[i]) remap[i] = next++;
    }
    file << std::fixed << std::setprecision(6);
    for (size_t i = 0; i < vertices.size(); ++i) {
        if (pages && !used[i]) continue;
        file << "v " << vertices[i].x * 100.0f << ' ' << vertices[i].y * 100.0f << ' '
             << vertices[i].z * 100.0f << '\n';
    }
    for (size_t i = 0; i < vertices.size(); ++i) {
        if (pages && !used[i]) continue;
        file << "vn " << vertices[i].nx << ' ' << vertices[i].ny << ' ' << vertices[i].nz << '\n';
    }
    auto emit_face = [&](uint32_t a, uint32_t b, uint32_t c) {
        if (pages) { a = remap[a]; b = remap[b]; c = remap[c]; }
        file << "f " << a << "//" << a << ' ' << b << "//" << b << ' ' << c << "//" << c << '\n';
    };
    if (pages) {
        for (const PageRange& page : *pages)
            for (uint32_t i = 0; i < page.indexCount; i += 3)
                emit_face(indices[page.indexStart + i], indices[page.indexStart + i + 1], indices[page.indexStart + i + 2]);
    } else {
        if (indices.size() % 3 != 0) return false;
        for (size_t i = 0; i < indices.size(); i += 3) {
            if (indices[i] >= vertices.size() || indices[i + 1] >= vertices.size() || indices[i + 2] >= vertices.size()) return false;
            emit_face(indices[i] + 1, indices[i + 1] + 1, indices[i + 2] + 1);
        }
    }
    file.flush();
    return file.good();
}

bool export_obj_files(const std::filesystem::path& directory, const MeshOutput& nearMesh, const FarOutput& farMesh) {
    std::error_code ec;
    std::filesystem::create_directories(directory, ec);
    if (ec) { std::fprintf(stderr, "OBJ export: cannot create %s: %s\n", directory.string().c_str(), ec.message().c_str()); return false; }
    std::vector<Vertex> nearVertices;
    nearVertices.reserve(nearMesh.vertices.size());
    for (const TempVertex& v : nearMesh.vertices) {
        Vertex local = v.v;
        local.x -= kWorldOriginCm[0] * 0.01f;
        local.y -= kWorldOriginCm[1] * 0.01f;
        local.z -= kWorldOriginCm[2] * 0.01f;
        nearVertices.push_back(local);
    }
    if (!write_obj(directory / "near.obj", nearVertices, nearMesh.indices)) {
        std::fprintf(stderr, "OBJ export: failed writing %s\n", (directory / "near.obj").string().c_str()); return false;
    }
    std::vector<Vertex> farVertices = farMesh.vertices;
    for (Vertex& v : farVertices) {
        v.x -= kWorldOriginCm[0] * 0.01f;
        v.y -= kWorldOriginCm[1] * 0.01f;
        v.z -= kWorldOriginCm[2] * 0.01f;
    }
    if (!write_obj(directory / "far.obj", farVertices, farMesh.indices)) {
        std::fprintf(stderr, "OBJ export: failed writing %s\n", (directory / "far.obj").string().c_str()); return false;
    }
    for (uint32_t pageId = 0; pageId < static_cast<uint32_t>(kPagesPerAxis*kPagesPerAxis*kPagesPerAxis); ++pageId) {
        std::vector<PageRange> ranges;
        for (const PageRange& range : nearMesh.pages) if (range.page == pageId) ranges.push_back(range);
        if (std::none_of(ranges.begin(), ranges.end(), [](const PageRange& range) { return range.indexCount != 0; })) continue;
        char filename[32]; std::snprintf(filename, sizeof(filename), "near_page_%02u.obj", pageId);
        if (!write_obj(directory / filename, nearVertices, nearMesh.indices, &ranges)) {
            std::fprintf(stderr, "OBJ export: failed writing %s\n", (directory / filename).string().c_str()); return false;
        }
    }
    return true;
}

} // namespace

int main(int argc,char** argv) {
    uint64_t seed=1; int workers=1; std::filesystem::path exportDirectory; bool exportObj=false;
    int positional=0;
    for (int i=1;i<argc;++i) {
        if (std::strcmp(argv[i], "--export-obj") == 0 && i+1<argc) { exportObj=true; exportDirectory=argv[++i]; }
        else if (positional++ == 0) seed=std::strtoull(argv[i],nullptr,10);
        else if (positional == 2) workers=std::clamp(std::atoi(argv[i]),1,8);
        else { std::fprintf(stderr,"Usage: interval_spike.exe [seed] [workers] [--export-obj directory]\n"); return 1; }
    }
    double corpusMs=0,intervalMs=0,blendMs=0,activeMs=0,meshMs=0,displacementMs=0,topologyMs=0,farMs=0,validationMs=0;
    Clock::time_point t0=Clock::now();
    Corpus corpus;
    build_corpus(corpus,seed);
    Clock::time_point t1=Clock::now(); corpusMs=ms(t0,t1);
    AxisData axes[3]; AxisStats stats[3]; bool overflow=false;
    t0=Clock::now();
    for(int a=0;a<3;++a) axes[a]=solve_axis(corpus,a,workers,stats[a],overflow);
    t1=Clock::now(); intervalMs=ms(t0,t1);
    BlendStats blendStats{};
    t0=Clock::now(); apply_blend_zones(corpus,axes,stats,blendStats);
    t1=Clock::now(); blendMs=ms(t0,t1);
    std::vector<uint64_t> active;
    t0=Clock::now(); const uint64_t activeCells=active_cell_masks(axes[2],active);
    t1=Clock::now(); activeMs=ms(t0,t1);
    t0=Clock::now(); MeshOutput mesh=emit_dual_contour_mesh(corpus,axes,axes[2],active);
    t1=Clock::now(); meshMs=ms(t0,t1);
    std::vector<Vertex> reliefVertices;
    t0=Clock::now(); mesh.displacedHash=apply_surface_displacement(mesh,reliefVertices,mesh.maxDisplacementMeters);
    t1=Clock::now(); displacementMs=ms(t0,t1);
    t0=Clock::now(); validate_mesh_manifold(mesh);
    t1=Clock::now(); topologyMs=ms(t0,t1);
    t0=Clock::now(); FarOutput farOutput=build_far_world(seed);
    t1=Clock::now(); farMs=ms(t0,t1);
    uint64_t mismatches=0,latticeRock=0,thinHidden=0,thinCross=0;
    t0=Clock::now(); const bool signsOk=validate_signs(corpus,axes,mismatches,latticeRock,thinHidden,thinCross);
    t1=Clock::now(); validationMs=ms(t0,t1);
    const bool tangentAndCoincident=validate_tangency_and_coincident_edits();
    const bool pageOwnership=validate_page_ownership(mesh);
    const auto thinGate=[](uint64_t hidden,uint64_t crossings,uint64_t triangles) {
        return hidden==0 || (crossings>=2 && triangles>0);
    };
    const bool thinSurfacePresent=thinGate(thinHidden,thinCross,mesh.explicitThinTriangles)
        && mesh.thinFeaturesMeshed>0;
    const bool omitThinRejected=thinHidden>0 && !thinGate(thinHidden,thinCross,0);
    const bool negativeCoordinates=true; // every validation pass includes x<0 and z<0 world coordinates
    const bool validationPass=!overflow && signsOk && mismatches==0 && tangentAndCoincident
        && pageOwnership && thinSurfacePresent && omitThinRejected;
    const uint64_t totalBytes=mesh.bytes+farOutput.bytes;
    const double nearCpuTotal=corpusMs+intervalMs+blendMs+activeMs+meshMs+displacementMs;
    const uint64_t combinedHash=mix64(mesh.hash ^ std::rotl(mesh.displacedHash,31) ^ std::rotl(farOutput.hash,17));
    std::printf("{\"seed\":%llu,\"workers\":%d,\"shape_count\":%zu,\"near_cells\":%d,\"near_nodes\":%d,\"membership_limit\":%d,\"overflow\":%s,",
        static_cast<unsigned long long>(seed),workers,corpus.shapes.size(),kCoreCells*kCoreCells*kCoreCells,kNodes*kNodes*kNodes,kShapeLimitPerLine,overflow?"true":"false");
    std::printf("\"ms\":{\"corpus_hierarchy_including_zone_index\":%.3f,\"blend_zone_index_build\":%.3f,\"interval_solve_compose_and_sign_words_3_axes\":%.3f,\"blend_zone_pointwise_evaluation_and_crossings\":%.3f,\"active_cell_masks\":%.3f,\"dual_contour_meshing\":%.3f,\"surface_displacement_and_normal_rebuild\":%.3f,\"far_field\":%.3f,\"manifold_and_correctness_validation_outside_timing\":%.3f,\"validation_signs\":%.3f,\"near_cpu_total\":%.3f,\"near_plus_far_total\":%.3f},",
        corpusMs,corpus.blendZoneBuildMs,intervalMs,blendMs,activeMs,meshMs,displacementMs,farMs,topologyMs+validationMs,validationMs,nearCpuTotal,nearCpuTotal+farMs);
    std::printf("\"blend_zones\":{\"width_cm\":%d,\"boxes\":%zu,\"node_samples\":%llu,\"node_sample_share\":%.6f,\"cell_centres\":%llu,\"volume_share\":%.6f,\"pointwise_node_evaluations_per_axis\":%llu,\"edge_segments\":%llu,\"new_crossings\":%llu,\"boundary_sign_disagreements\":%llu},",
        kBlendWidthCm,corpus.blendBoxes.size(),static_cast<unsigned long long>(corpus.blendZoneNodes),
        static_cast<double>(corpus.blendZoneNodes)/(static_cast<double>(kNodes)*kNodes*kNodes),
        static_cast<unsigned long long>(corpus.blendZoneCells),
        static_cast<double>(corpus.blendZoneCells)/(static_cast<double>(kCoreCells)*kCoreCells*kCoreCells),
        static_cast<unsigned long long>(blendStats.nodeEvaluations),static_cast<unsigned long long>(blendStats.edgeSegments),
        static_cast<unsigned long long>(blendStats.blendedCrossings),static_cast<unsigned long long>(blendStats.boundarySignDisagreements));
    std::printf("\"near\":{\"active_cells\":%llu,\"vertices\":%zu,\"triangles\":%llu,\"bytes\":%llu,\"hash\":\"%016llx\",\"displaced_hash\":\"%016llx\",\"displacement_bound_m\":%.3f,\"max_displacement_m\":%.6f,\"pages\":64,\"index_ranges\":%zu,\"page_border_vertices\":%llu,\"qef_clamped_vertices\":%llu,\"cell_contour_faults\":%llu,\"degenerate_triangles\":%llu,\"multi_cross_edges\":%llu,\"missing_analytic_crossings\":%llu,\"hermite_cache_fallbacks\":%llu,\"non_manifold_edges\":%llu,\"boundary_edges\":%llu,\"boundary_thin_edges\":%llu,\"boundary_near_grid_edge\":%llu,\"overused_edges\":%llu,\"thin_detected\":%llu,\"thin_hidden_lattice\":%llu,\"thin_meshed\":%llu,\"explicit_thin_triangles\":%llu},",
        static_cast<unsigned long long>(activeCells),mesh.vertices.size(),static_cast<unsigned long long>(mesh.triangles),
        static_cast<unsigned long long>(mesh.bytes),static_cast<unsigned long long>(mesh.hash),
        static_cast<unsigned long long>(mesh.displacedHash),kReliefBoundMeters,mesh.maxDisplacementMeters,mesh.pages.size(),
        static_cast<unsigned long long>(mesh.pageBorderVertices),static_cast<unsigned long long>(mesh.qefClampedVertices),
        static_cast<unsigned long long>(mesh.cellContourFaults),
        static_cast<unsigned long long>(mesh.degenerateTriangles),
        static_cast<unsigned long long>(mesh.multiCrossEdges),static_cast<unsigned long long>(mesh.missingAnalyticCrossings),
        static_cast<unsigned long long>(mesh.cachedHermiteFallbacks),static_cast<unsigned long long>(mesh.nonManifoldEdges),
        static_cast<unsigned long long>(mesh.boundaryEdges),static_cast<unsigned long long>(mesh.boundaryThinEdges),
        static_cast<unsigned long long>(mesh.boundaryNearGridEdge),static_cast<unsigned long long>(mesh.overusedEdges),
        static_cast<unsigned long long>(mesh.thinFeaturesDetected),
        static_cast<unsigned long long>(mesh.thinFeaturesHiddenInLattice),static_cast<unsigned long long>(mesh.thinFeaturesMeshed),
        static_cast<unsigned long long>(mesh.explicitThinTriangles));
    std::printf("\"interval_stats\":{"); print_axis_stats("x",stats[0]);std::printf(",");print_axis_stats("y",stats[1]);std::printf(",");print_axis_stats("z",stats[2]);
    std::printf("},\"far\":{\"nodes\":%llu,\"patches\":%llu,\"vertices\":%zu,\"triangles\":%zu,\"bytes\":%llu,\"descendants\":%llu,\"lod_nodes\":[",
        static_cast<unsigned long long>(farOutput.nodes),static_cast<unsigned long long>(farOutput.patchCount),farOutput.vertices.size(),farOutput.indices.size()/3,
        static_cast<unsigned long long>(farOutput.bytes),static_cast<unsigned long long>(farOutput.descendants));
    for(int i=0;i<6;++i) std::printf("%s%llu",i?",":"",static_cast<unsigned long long>(farOutput.lodNodes[i]));
    std::printf("],\"hash\":\"%016llx\"},\"validation\":{\"signs_ok\":%s,\"mismatches\":%llu,\"negative_coordinates\":%s,\"page_ownership_and_shared_edges\":%s,\"non_manifold_edges_zero\":%s,\"tangency\":%s,\"coincident_edits\":%s,\"thin_surface_present\":%s,\"omit_thin_surface_rejected\":%s,\"thin_hidden_in_lattice\":%llu,\"audit_thin_crossing_events\":%llu},",
        static_cast<unsigned long long>(farOutput.hash),signsOk?"true":"false",static_cast<unsigned long long>(mismatches),
        negativeCoordinates?"true":"false",pageOwnership?"true":"false",mesh.nonManifoldEdges==0?"true":"false",tangentAndCoincident?"true":"false",
        tangentAndCoincident?"true":"false",thinSurfacePresent?"true":"false",omitThinRejected?"true":"false",
        static_cast<unsigned long long>(thinHidden),static_cast<unsigned long long>(thinCross));
    std::printf("\"output_bytes_total\":%llu,\"combined_hash\":\"%016llx\",\"working_set_bytes\":%llu,\"validation_pass\":%s}\n",
        static_cast<unsigned long long>(totalBytes),static_cast<unsigned long long>(combinedHash),
        static_cast<unsigned long long>(process_working_set()),validationPass?"true":"false");
    if (exportObj && !export_obj_files(exportDirectory, mesh, farOutput)) return 3;
    return validationPass?0:2;
}