#include "CoreMinimal.h"
#include "MarchingCubesTables.h"
#include "VoxelForgeCore/VoxelForgeCore.h"

#include <algorithm>
#include <array>
#include <atomic>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <filesystem>
#include <fstream>
#include <functional>
#include <iomanip>
#include <iostream>
#include <map>
#include <numeric>
#include <set>
#include <string>
#include <thread>
#include <tuple>
#include <unordered_map>
#include <utility>
#include <vector>

namespace
{
    constexpr double VoxelMeters = 0.25;
    constexpr double PositionQuantum = VoxelMeters / 65536.0;
    constexpr double LodContourErrorBudgetMeters = 0.50; // One LOD1 cell (50 cm) physical mismatch.
    constexpr double NormalSplitCos = 0.86602540378443864676; // 30 degrees
    constexpr int CornerCount = 8;
    constexpr int EdgeCount = 12;
    constexpr std::array<std::array<int, 3>, CornerCount> CornerOffsets{{
        {{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 std::array<std::array<int, 2>, EdgeCount> EdgeCorners{{
        {{0,1}}, {{1,2}}, {{2,3}}, {{3,0}}, {{4,5}}, {{5,6}},
        {{6,7}}, {{7,4}}, {{0,4}}, {{1,5}}, {{2,6}}, {{3,7}}
    }};

    struct V3
    {
        double x=0, y=0, z=0;
        V3 operator+(V3 b) const { return {x+b.x,y+b.y,z+b.z}; }
        V3 operator-(V3 b) const { return {x-b.x,y-b.y,z-b.z}; }
        V3 operator*(double s) const { return {x*s,y*s,z*s}; }
        V3 operator/(double s) const { return {x/s,y/s,z/s}; }
        V3& operator+=(V3 b) { x+=b.x; y+=b.y; z+=b.z; return *this; }
    };
    double Dot(V3 a,V3 b){return a.x*b.x+a.y*b.y+a.z*b.z;}
    V3 Cross(V3 a,V3 b){return {a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};}
    double Length(V3 a){return std::sqrt(Dot(a,a));}
    V3 Normalize(V3 a)
    {
        const double n=Length(a);
        return n>1.0e-20?a/n:V3{0,0,1};
    }
    V3 Quantize(V3 p)
    {
        const auto q=[](double x){return std::nearbyint(x/PositionQuantum)*PositionQuantum;};
        return {q(p.x),q(p.y),q(p.z)};
    }

    template<class F> void ParallelFor(std::size_t count,int workers,F&& function)
    {
        workers=std::clamp(workers,1,16);
        if(workers==1||count<512)
        {
            for(std::size_t i=0;i<count;++i)function(i);
            return;
        }
        const int actual=std::min<int>(workers,static_cast<int>(count));
        std::vector<std::thread> threads;
        threads.reserve(actual);
        for(int worker=0;worker<actual;++worker)
        {
            const std::size_t begin=count*static_cast<std::size_t>(worker)/actual;
            const std::size_t end=count*static_cast<std::size_t>(worker+1)/actual;
            threads.emplace_back([begin,end,&function](){for(std::size_t i=begin;i<end;++i)function(i);});
        }
        for(auto& thread:threads)thread.join();
    }

    struct Plane
    {
        std::uint64_t id=0;
        V3 normal{};
    };
    struct FieldSample
    {
        double density=1.0;
        std::array<Plane,16> planes{};
        std::uint8_t planeCount=0;
        V3 gradient{};
        bool sharp() const { return planeCount!=0; }
        std::uint64_t primaryFeature() const
        {
            if(!planeCount)return 0;
            std::uint64_t id=planes[0].id;
            for(std::uint8_t i=1;i<planeCount;++i)id=std::min(id,planes[i].id);
            return id;
        }
    };
    struct Box
    {
        V3 low{},high{};
        std::uint32_t id=0;
        bool inverted=false;
    };
    struct BoxSample
    {
        double density=1.0;
        std::array<Plane,3> planes{};
        std::uint8_t planeCount=0;
        V3 gradient{};
    };

    BoxSample EvalBox(const Box& box,V3 p)
    {
        const V3 center=(box.low+box.high)*0.5;
        const V3 half=(box.high-box.low)*0.5;
        const V3 delta=p-center;
        const V3 q{std::abs(delta.x)-half.x,std::abs(delta.y)-half.y,std::abs(delta.z)-half.z};
        const V3 outside{std::max(q.x,0.0),std::max(q.y,0.0),std::max(q.z,0.0)};
        BoxSample out;
        out.density=Length(outside)+std::min(std::max(q.x,std::max(q.y,q.z)),0.0);
        const double maxQ=std::max(q.x,std::max(q.y,q.z));
        V3 gradient{};
        if(Length(outside)>1.0e-14)
        {
            gradient=Normalize({outside.x*(delta.x<0?-1:1),outside.y*(delta.y<0?-1:1),outside.z*(delta.z<0?-1:1)});
        }
        for(int axis=0;axis<3;++axis)
        {
            const double qAxis=axis==0?q.x:(axis==1?q.y:q.z);
            if(std::abs(qAxis-maxQ)>1.0e-10)continue;
            const double dAxis=axis==0?delta.x:(axis==1?delta.y:delta.z);
            V3 normal{};
            if(axis==0)normal.x=dAxis<0?-1:1;
            if(axis==1)normal.y=dAxis<0?-1:1;
            if(axis==2)normal.z=dAxis<0?-1:1;
            const std::uint64_t side=dAxis<0?0ull:1ull;
            const std::uint64_t feature=(static_cast<std::uint64_t>(box.id)<<4)
                |static_cast<std::uint64_t>(axis*2)+side+1ull;
            if(box.inverted)normal=normal*-1.0;
            out.planes[out.planeCount++]={feature,normal};
            gradient+=normal;
        }
        if(box.inverted)
        {
            out.density=-out.density;
            gradient=gradient*-1.0;
        }
        out.gradient=Length(gradient)>1.0e-14?Normalize(gradient):Normalize(outside);
        return out;
    }

    enum class Piece : int { Corner=0, Doorway=1, Step=2, Ledge=3, Organic=4, NativeSurface=5 };
    const char* PieceName(Piece p)
    {
        switch(p)
        {
            case Piece::Corner:return "corner90";
            case Piece::Doorway:return "doorway";
            case Piece::Step:return "riser1m";
            case Piece::Ledge:return "cliff_overhang";
            case Piece::Organic:return "organic_cave_chamber";
            default:return "native_surface_cliff";
        }
    }
    struct Field
    {
        Piece piece=Piece::Corner;
        V3 phaseShift{};
        const VoxelForgeCore::Heightfield2D* nativeSurface=nullptr;
        std::array<Box,4> boxes{};
        std::uint8_t boxCount=0;
        Box doorVoid{};
        void AddBox(V3 low,V3 high,std::uint32_t id)
        { boxes[boxCount++]={low,high,id,false}; }
        explicit Field(Piece p,V3 phase={}):piece(p),phaseShift(phase)
        {
            switch(p)
            {
                case Piece::Corner:
                    AddBox({-2.5,-2.5,-1.5},{0.25,2.5,1.5},1);
                    AddBox({-2.5,-0.25,-1.5},{2.5,2.5,1.5},2);
                    break;
                case Piece::Doorway:
                    AddBox({-0.5,-2.5,-2.0},{0.5,2.5,3.0},3);
                    doorVoid={{-0.8,-0.85,-1.75},{0.8,0.85,2.0},4,true};
                    break;
                case Piece::Step:
                    AddBox({-3.0,-3.0,-2.0},{3.0,0.25,0.0},5);
                    AddBox({-3.0,-0.25,-2.0},{3.0,3.0,1.0},6);
                    break;
                case Piece::Ledge:
                    AddBox({-3.0,-3.0,-2.0},{0.25,3.0,3.0},7);
                    AddBox({-0.25,-3.0,1.5},{2.2,3.0,2.0},8);
                    break;
                case Piece::Organic:case Piece::NativeSurface:break;
            }
        }
        explicit Field(const VoxelForgeCore::Heightfield2D& surface,V3 phase={})
            :piece(Piece::NativeSurface),phaseShift(phase),nativeSurface(&surface){}

        FieldSample Eval(V3 p) const
        {
            p=p+phaseShift;
            if(nativeSurface)
            {
                constexpr double SubvoxelsPerMeter=1024.0;
                constexpr std::int64_t GradientDeltaSubvoxels=64;
                const auto x=static_cast<std::int64_t>(std::llround(p.x*SubvoxelsPerMeter));
                const auto y=static_cast<std::int64_t>(std::llround(p.y*SubvoxelsPerMeter));
                const auto z=static_cast<std::int64_t>(std::llround(p.z*SubvoxelsPerMeter));
                const auto center=nativeSurface->point(x,y,z);
                FieldSample result;result.density=static_cast<double>(center.densitySubvoxels)/SubvoxelsPerMeter;
                const auto xp=nativeSurface->point(x+GradientDeltaSubvoxels,y,z).densitySubvoxels;
                const auto xm=nativeSurface->point(x-GradientDeltaSubvoxels,y,z).densitySubvoxels;
                const auto yp=nativeSurface->point(x,y+GradientDeltaSubvoxels,z).densitySubvoxels;
                const auto ym=nativeSurface->point(x,y-GradientDeltaSubvoxels,z).densitySubvoxels;
                const auto zp=nativeSurface->point(x,y,z+GradientDeltaSubvoxels).densitySubvoxels;
                const auto zm=nativeSurface->point(x,y,z-GradientDeltaSubvoxels).densitySubvoxels;
                result.gradient=Normalize({static_cast<double>(xp-xm),static_cast<double>(yp-ym),static_cast<double>(zp-zm)});
                return result;
            }
            if(piece==Piece::Organic)
            {
                const V3 d=p-V3{0,0,0};
                FieldSample result;
                // Exact spherical chamber cavity: negative solid outside, positive air inside.
                result.density=2.0-Length(d);
                result.gradient=Normalize(d)*-1.0;
                return result;
            }
            FieldSample result;
            result.density=std::numeric_limits<double>::infinity();
            for(std::uint8_t i=0;i<boxCount;++i)
            {
                const BoxSample candidate=EvalBox(boxes[i],p);
                if(candidate.density<result.density-1.0e-10)
                {
                    result.density=candidate.density;
                    result.planeCount=candidate.planeCount;
                    for(std::uint8_t k=0;k<candidate.planeCount;++k)result.planes[k]=candidate.planes[k];
                }
                else if(std::abs(candidate.density-result.density)<=1.0e-10)
                {
                    for(std::uint8_t k=0;k<candidate.planeCount&&result.planeCount<result.planes.size();++k)
                        result.planes[result.planeCount++]=candidate.planes[k];
                }
            }
            if(piece==Piece::Doorway)
            {
                const BoxSample cut=EvalBox(doorVoid,p);
                const double cutDensity=cut.density;
                if(cutDensity>result.density+1.0e-10)
                {
                    result.density=cutDensity;
                    result.planeCount=cut.planeCount;
                    for(std::uint8_t k=0;k<cut.planeCount;++k)result.planes[k]=cut.planes[k];
                }
                else if(std::abs(cutDensity-result.density)<=1.0e-10)
                {
                    for(std::uint8_t k=0;k<cut.planeCount&&result.planeCount<result.planes.size();++k)
                        result.planes[result.planeCount++]=cut.planes[k];
                }
            }
            std::sort(result.planes.begin(),result.planes.begin()+result.planeCount,
                [](const Plane& a,const Plane& b){return a.id<b.id;});
            std::uint8_t unique=0;
            for(std::uint8_t i=0;i<result.planeCount;++i)
                if(unique==0||result.planes[i].id!=result.planes[unique-1].id)
                    result.planes[unique++]=result.planes[i];
            result.planeCount=unique;
            V3 sum{};
            for(std::uint8_t i=0;i<result.planeCount;++i)sum+=result.planes[i].normal;
            result.gradient=Length(sum)>1.0e-14?Normalize(sum):V3{0,0,1};
            return result;
        }
    };

    struct Grid
    {
        int nx=0,ny=0,nz=0,step=1;
        std::array<int,3> coreMin{{0,0,0}},coreMax{{0,0,0}};
        V3 origin{};
        std::vector<double> density;
        std::size_t NodeIndex(int x,int y,int z) const
        { return (static_cast<std::size_t>(z)*(ny+1)+y)*(nx+1)+x; }
        V3 Position(int x,int y,int z) const
        { return origin+V3{x*step*VoxelMeters,y*step*VoxelMeters,z*step*VoxelMeters}; }
        std::size_t CellIndex(int x,int y,int z) const
        { return (static_cast<std::size_t>(z)*ny+y)*nx+x; }
        bool IsCoreCell(int x,int y,int z) const
        { return x>=coreMin[0]&&y>=coreMin[1]&&z>=coreMin[2]&&
            x<coreMax[0]&&y<coreMax[1]&&z<coreMax[2]; }
        std::size_t CellCount() const { return static_cast<std::size_t>(nx)*ny*nz; }
    };
    Grid SampleGrid(const Field& field,int nx,int ny,int nz,int step,V3 origin,int workers,int halo=0)
    {
        Grid grid;grid.step=step;
        grid.coreMin={halo,halo,halo};grid.coreMax={halo+nx,halo+ny,halo+nz};
        grid.nx=nx+2*halo;grid.ny=ny+2*halo;grid.nz=nz+2*halo;
        grid.origin=origin-V3{halo*step*VoxelMeters,halo*step*VoxelMeters,halo*step*VoxelMeters};
        const std::size_t count=static_cast<std::size_t>(grid.nx+1)*(grid.ny+1)*(grid.nz+1);
        grid.density.resize(count);
        ParallelFor(count,workers,[&](std::size_t index)
        {
            const int x=static_cast<int>(index%(grid.nx+1));
            const int y=static_cast<int>((index/(grid.nx+1))%(grid.ny+1));
            const int z=static_cast<int>(index/((grid.nx+1)*(grid.ny+1)));
            const double density=field.Eval(grid.Position(x,y,z)).density;
            // Stable iso-level tie rule: zero remains air, matching the Core sign convention, but
            // receives a sub-cell positive bias so edge roots do not collapse onto one lattice
            // vertex from several unrelated edges. This is only 0.01 mm and is below the output
            // surface-error precision; it is shared by all three candidates.
            grid.density[index]=density==0.0?1.0e-5:density;
        });
        return grid;
    }

    struct Root
    {
        bool active=false;
        V3 position{},gradient{};
        std::array<Plane,16> planes{};
        std::uint8_t planeCount=0;
        std::uint64_t feature=0;
    };
    struct EdgeArrays
    {
        std::array<std::vector<Root>,3> roots;
        std::array<std::vector<int>,3> vertexIds;
        std::array<std::size_t,3> counts{};
        std::size_t Total() const{return counts[0]+counts[1]+counts[2];}
        std::size_t Index(int axis,int x,int y,int z,const Grid& g) const
        {
            if(axis==0)return (static_cast<std::size_t>(z)*(g.ny+1)+y)*g.nx+x;
            if(axis==1)return (static_cast<std::size_t>(z)*g.ny+y)*(g.nx+1)+x;
            return (static_cast<std::size_t>(z)*(g.ny+1)+y)*(g.nx+1)+x;
        }
        Root& Get(int axis,int x,int y,int z,const Grid& g){return roots[axis][Index(axis,x,y,z,g)];}
        const Root& Get(int axis,int x,int y,int z,const Grid& g) const{return roots[axis][Index(axis,x,y,z,g)];}
        int& Vertex(int axis,int x,int y,int z,const Grid& g){return vertexIds[axis][Index(axis,x,y,z,g)];}
        int Vertex(int axis,int x,int y,int z,const Grid& g) const{return vertexIds[axis][Index(axis,x,y,z,g)];}
    };

    EdgeArrays BuildRoots(const Grid& grid,const Field& field,int workers)
    {
        EdgeArrays edges;
        edges.counts={static_cast<std::size_t>(grid.nx)*(grid.ny+1)*(grid.nz+1),
            static_cast<std::size_t>(grid.nx+1)*grid.ny*(grid.nz+1),
            static_cast<std::size_t>(grid.nx+1)*(grid.ny+1)*grid.nz};
        for(int axis=0;axis<3;++axis)
        {
            edges.roots[axis].resize(edges.counts[axis]);
            edges.vertexIds[axis].assign(edges.counts[axis],-1);
        }
        const std::size_t total=edges.Total();
        ParallelFor(total,workers,[&](std::size_t flat)
        {
            int axis=0;std::size_t index=flat;
            while(axis<2&&index>=edges.counts[axis]){index-=edges.counts[axis];++axis;}
            int x=0,y=0,z=0;
            if(axis==0)
            {
                x=static_cast<int>(index%grid.nx);y=static_cast<int>((index/grid.nx)%(grid.ny+1));
                z=static_cast<int>(index/(grid.nx*(grid.ny+1)));
            }
            else if(axis==1)
            {
                x=static_cast<int>(index%(grid.nx+1));y=static_cast<int>((index/(grid.nx+1))%grid.ny);
                z=static_cast<int>(index/((grid.nx+1)*grid.ny));
            }
            else
            {
                x=static_cast<int>(index%(grid.nx+1));y=static_cast<int>((index/(grid.nx+1))%(grid.ny+1));
                z=static_cast<int>(index/((grid.nx+1)*(grid.ny+1)));
            }
            int x2=x,y2=y,z2=z;
            if(axis==0)++x2;else if(axis==1)++y2;else ++z2;
            const double d0=grid.density[grid.NodeIndex(x,y,z)];
            const double d1=grid.density[grid.NodeIndex(x2,y2,z2)];
            if((d0<0.0)==(d1<0.0))return;
            double t=d0/(d0-d1);t=std::clamp(t,0.0,1.0);
            Root root;root.active=true;
            root.position=Quantize(grid.Position(x,y,z)+(grid.Position(x2,y2,z2)-grid.Position(x,y,z))*t);
            const FieldSample sample=field.Eval(root.position);
            const double h=grid.step*VoxelMeters;
            const auto gradientAt=[&](V3 point)
            {
                const double dx=field.Eval(point+V3{h,0,0}).density-field.Eval(point-V3{h,0,0}).density;
                const double dy=field.Eval(point+V3{0,h,0}).density-field.Eval(point-V3{0,h,0}).density;
                const double dz=field.Eval(point+V3{0,0,h}).density-field.Eval(point-V3{0,0,h}).density;
                return Normalize({dx,dy,dz});
            };
            const V3 g0=gradientAt(grid.Position(x,y,z));
            const V3 g1=gradientAt(grid.Position(x2,y2,z2));
            root.gradient=Normalize(g0*(1.0-t)+g1*t);
            root.planeCount=sample.planeCount;
            for(std::uint8_t i=0;i<root.planeCount;++i)root.planes[i]=sample.planes[i];
            root.feature=sample.primaryFeature();
            edges.roots[axis][index]=root;
        });
        return edges;
    }

    struct Vertex
    {
        V3 position{},normal{};
    };
    struct Triangle
    {
        std::array<std::uint32_t,3> index{};
        std::uint64_t feature=0;
    };
    struct Mesh
    {
        std::vector<Vertex> vertices;
        std::vector<Triangle> triangles;
        std::uint64_t hash=0;
        std::uint64_t skippedDualEdges=0;
    };
    enum class Method : int { MC=0, Hybrid=1, DC=2 };
    const char* MethodName(Method m){return m==Method::MC?"MC":m==Method::Hybrid?"Hybrid": "DC";}

    struct Cell
    {
        std::uint8_t cubeCase=0;
        bool active=false;
        std::uint8_t componentCount=0;
        std::array<std::uint8_t,EdgeCount> edgeComponent{};
        std::array<bool,EdgeCount> componentSharp{},componentUseDC{};
        std::array<V3,EdgeCount> dual{};
    };

    int FindLocal(std::array<std::uint8_t,EdgeCount>& parent,int value)
    {
        while(parent[value]!=value)
        {
            parent[value]=parent[parent[value]];
            value=parent[value];
        }
        return value;
    }
    void UniteLocal(std::array<std::uint8_t,EdgeCount>& parent,int a,int b)
    {
        a=FindLocal(parent,a);b=FindLocal(parent,b);if(a==b)return;
        if(a>b)std::swap(a,b);
        parent[b]=static_cast<std::uint8_t>(a);
    }

    struct ComponentSets
    {
        std::vector<std::uint32_t> parent;
        std::vector<std::uint8_t> rank,sharp;
        explicit ComponentSets(std::size_t size):parent(size),rank(size,0),sharp(size,0)
        {std::iota(parent.begin(),parent.end(),0u);}
        std::uint32_t Find(std::uint32_t value)
        {
            while(parent[value]!=value){parent[value]=parent[parent[value]];value=parent[value];}
            return value;
        }
        void Unite(std::uint32_t a,std::uint32_t b)
        {
            a=Find(a);b=Find(b);if(a==b)return;
            if(rank[a]<rank[b])std::swap(a,b);
            parent[b]=a;sharp[a]=static_cast<std::uint8_t>(sharp[a]|sharp[b]);
            if(rank[a]==rank[b])++rank[a];
        }
    };

    std::array<int,4> LocalEdge(int x,int y,int z,int edge)
    {
        const auto c=EdgeCorners[edge];const auto a=CornerOffsets[c[0]],b=CornerOffsets[c[1]];
        const int axis=a[0]!=b[0]?0:a[1]!=b[1]?1:2;
        return {axis,x+std::min(a[0],b[0]),y+std::min(a[1],b[1]),z+std::min(a[2],b[2])};
    }
    int FindLocalEdge(int cellX,int cellY,int cellZ,int axis,int edgeX,int edgeY,int edgeZ)
    {
        for(int edge=0;edge<EdgeCount;++edge)
        {
            const auto local=LocalEdge(cellX,cellY,cellZ,edge);
            if(local[0]==axis&&local[1]==edgeX&&local[2]==edgeY&&local[3]==edgeZ)return edge;
        }
        return -1;
    }

    V3 SolveQEF(const std::vector<std::pair<V3,V3>>& constraints,V3 fallback,V3 low,V3 high)
    {
        if(constraints.empty())return fallback;
        double a[3][4]{};
        // Fixed constraint order; a small center penalty resolves rank-deficient planes without
        // moving a two/three-plane intersection at measurable scale.
        constexpr double lambda=1.0e-10;
        const double c[3]={fallback.x,fallback.y,fallback.z};
        for(int i=0;i<3;++i){a[i][i]=lambda;a[i][3]=lambda*c[i];}
        for(const auto& constraint:constraints)
        {
            const V3 n=Normalize(constraint.second);const double b=Dot(n,constraint.first);
            const double nv[3]={n.x,n.y,n.z};
            for(int r=0;r<3;++r){for(int k=0;k<3;++k)a[r][k]+=nv[r]*nv[k];a[r][3]+=nv[r]*b;}
        }
        bool singular=false;
        for(int column=0;column<3;++column)
        {
            int pivot=column;
            for(int row=column+1;row<3;++row)if(std::abs(a[row][column])>std::abs(a[pivot][column]))pivot=row;
            if(std::abs(a[pivot][column])<1.0e-14){singular=true;break;}
            if(pivot!=column)for(int k=column;k<4;++k)std::swap(a[pivot][k],a[column][k]);
            const double divisor=a[column][column];
            for(int k=column;k<4;++k)a[column][k]/=divisor;
            for(int row=0;row<3;++row)if(row!=column)
            {
                const double factor=a[row][column];
                for(int k=column;k<4;++k)a[row][k]-=factor*a[column][k];
            }
        }
        if(singular)return fallback;
        return Quantize({std::clamp(a[0][3],low.x,high.x),
            std::clamp(a[1][3],low.y,high.y),std::clamp(a[2][3],low.z,high.z)});
    }

    std::uint64_t SelectFeature(const Root& root,int edgeAxis)
    {
        if(!root.planeCount)return 0;
        std::uint8_t selected=0;
        for(std::uint8_t i=1;i<root.planeCount;++i)
        {
            const double a=edgeAxis==0?std::abs(root.planes[i].normal.x):
                edgeAxis==1?std::abs(root.planes[i].normal.y):std::abs(root.planes[i].normal.z);
            const double b=edgeAxis==0?std::abs(root.planes[selected].normal.x):
                edgeAxis==1?std::abs(root.planes[selected].normal.y):std::abs(root.planes[selected].normal.z);
            if(a>b+1.0e-12||(std::abs(a-b)<=1.0e-12&&root.planes[i].id<root.planes[selected].id))selected=i;
        }
        return root.planes[selected].id;
    }
    V3 NormalForFeature(const Root& root,std::uint64_t id)
    {
        for(std::uint8_t i=0;i<root.planeCount;++i)if(root.planes[i].id==id)return root.planes[i].normal;
        return root.gradient;
    }

    void AddTriangle(Mesh& mesh,std::uint32_t a,std::uint32_t b,std::uint32_t c,
        std::uint64_t feature,V3 outward)
    {
        if(a==b||b==c||c==a)return;
        const V3 pa=mesh.vertices[a].position,pb=mesh.vertices[b].position,pc=mesh.vertices[c].position;
        if(Length(Cross(pb-pa,pc-pa))<1.0e-14)return;
        if(Dot(Cross(pb-pa,pc-pa),outward)<0.0)std::swap(b,c);
        mesh.triangles.push_back({{a,b,c},feature});
    }

    std::uint64_t HashMesh(const Mesh& mesh)
    {
        std::uint64_t h=14695981039346656037ull;
        const auto add=[&](std::uint64_t word,int bytes)
        {
            for(int i=0;i<bytes;++i){h^=static_cast<std::uint8_t>(word>>(i*8));h*=1099511628211ull;}
        };
        add(mesh.vertices.size(),8);
        for(const Vertex& v:mesh.vertices)
        {
            for(double x:{v.position.x,v.position.y,v.position.z,v.normal.x,v.normal.y,v.normal.z})
                add(std::bit_cast<std::uint64_t>(x),8);
        }
        add(mesh.triangles.size(),8);
        for(const Triangle& t:mesh.triangles)
        {add(t.index[0],4);add(t.index[1],4);add(t.index[2],4);add(t.feature,8);}
        return h;
    }

    Mesh BuildMesh(const Grid& grid,const Field& field,Method method,int workers,bool enforceQuadOwnership=true)
    {
        EdgeArrays edges=BuildRoots(grid,field,workers);
        std::vector<Cell> cells(grid.CellCount());
        ParallelFor(cells.size(),workers,[&](std::size_t ci)
        {
            const int x=static_cast<int>(ci%grid.nx);
            const int y=static_cast<int>((ci/grid.nx)%grid.ny);
            const int z=static_cast<int>(ci/(grid.nx*grid.ny));
            Cell& cell=cells[ci];cell.edgeComponent.fill(0xff);
            double values[CornerCount];
            for(int corner=0;corner<CornerCount;++corner)
            {
                const auto o=CornerOffsets[corner];
                values[corner]=grid.density[grid.NodeIndex(x+o[0],y+o[1],z+o[2])];
                if(values[corner]<0.0)cell.cubeCase|=static_cast<std::uint8_t>(1u<<corner);
            }
            if(cell.cubeCase==0||cell.cubeCase==255)return;
            cell.active=true;
            if(method==Method::MC)
            {
                cell.componentCount=1;
                for(int edge=0;edge<EdgeCount;++edge)
                {
                    const auto ec=EdgeCorners[edge];
                    if((values[ec[0]]<0.0)!=(values[ec[1]]<0.0))cell.edgeComponent[edge]=0;
                }
                return;
            }
            std::array<std::uint8_t,EdgeCount> localParent{};
            for(int edge=0;edge<EdgeCount;++edge)localParent[edge]=static_cast<std::uint8_t>(edge);
            for(int k=0;k<16&&MCTriTable[cell.cubeCase][k]>=0;k+=3)
            {
                const int a=MCTriTable[cell.cubeCase][k],b=MCTriTable[cell.cubeCase][k+1],c=MCTriTable[cell.cubeCase][k+2];
                UniteLocal(localParent,a,b);UniteLocal(localParent,b,c);
            }
            std::array<std::uint8_t,EdgeCount> rootToComponent{};rootToComponent.fill(0xff);
            for(int edge=0;edge<EdgeCount;++edge)
            {
                const auto ec=EdgeCorners[edge];
                if((values[ec[0]]<0.0)==(values[ec[1]]<0.0))continue;
                const auto local=LocalEdge(x,y,z,edge);
                const Root& root=edges.Get(local[0],local[1],local[2],local[3],grid);
                if(!root.active)continue;
                const int localRoot=FindLocal(localParent,edge);
                if(rootToComponent[localRoot]==0xff)
                    rootToComponent[localRoot]=cell.componentCount++;
                const std::uint8_t component=rootToComponent[localRoot];
                cell.edgeComponent[edge]=component;
                cell.componentSharp[component]|=root.planeCount!=0;
            }
        });

        // Link local surface patches only through the same sign-changing primal edge. A cube may
        // contain multiple disconnected patches; keeping their QEFs separate avoids DC pinches.
        std::unique_ptr<ComponentSets> patches;
        if(method==Method::Hybrid)
        {
            patches=std::make_unique<ComponentSets>(cells.size()*EdgeCount);
            for(std::size_t cellIndex=0;cellIndex<cells.size();++cellIndex)
                for(std::uint8_t component=0;component<cells[cellIndex].componentCount;++component)
                    patches->sharp[cellIndex*EdgeCount+component]=cells[cellIndex].componentSharp[component]?1:0;
            for(int axis=0;axis<3;++axis)
            {
                for(std::size_t i=0;i<edges.roots[axis].size();++i)
                {
                    if(!edges.roots[axis][i].active)continue;
                    int x=0,y=0,z=0;
                    if(axis==0){x=static_cast<int>(i%grid.nx);y=static_cast<int>((i/grid.nx)%(grid.ny+1));z=static_cast<int>(i/(grid.nx*(grid.ny+1)));}
                    else if(axis==1){x=static_cast<int>(i%(grid.nx+1));y=static_cast<int>((i/(grid.nx+1))%grid.ny);z=static_cast<int>(i/((grid.nx+1)*grid.ny));}
                    else{x=static_cast<int>(i%(grid.nx+1));y=static_cast<int>((i/(grid.nx+1))%(grid.ny+1));z=static_cast<int>(i/((grid.nx+1)*(grid.ny+1)));}
                    std::array<std::array<int,3>,4> around{};
                    if(axis==0)around={{{x,y-1,z-1},{x,y,z-1},{x,y,z},{x,y-1,z}}};
                    else if(axis==1)around={{{x-1,y,z-1},{x-1,y,z},{x,y,z},{x,y,z-1}}};
                    else around={{{x-1,y-1,z},{x,y-1,z},{x,y,z},{x-1,y,z}}};
                    std::array<std::uint32_t,4> patchIds{};bool valid=true;
                    for(int corner=0;corner<4;++corner)
                    {
                        const auto& c=around[corner];
                        if(c[0]<0||c[1]<0||c[2]<0||c[0]>=grid.nx||c[1]>=grid.ny||c[2]>=grid.nz){valid=false;break;}
                        const std::size_t cellIndex=grid.CellIndex(c[0],c[1],c[2]);
                        const int localEdge=FindLocalEdge(c[0],c[1],c[2],axis,x,y,z);
                        if(localEdge<0||cells[cellIndex].edgeComponent[localEdge]==0xff){valid=false;break;}
                        patchIds[corner]=static_cast<std::uint32_t>(cellIndex*EdgeCount+cells[cellIndex].edgeComponent[localEdge]);
                    }
                    if(valid)for(int corner=1;corner<4;++corner)patches->Unite(patchIds[0],patchIds[corner]);
                }
            }
        }
        for(std::size_t cellIndex=0;cellIndex<cells.size();++cellIndex)
            for(std::uint8_t component=0;component<cells[cellIndex].componentCount;++component)
            {
                const std::uint32_t patch=static_cast<std::uint32_t>(cellIndex*EdgeCount+component);
                cells[cellIndex].componentUseDC[component]=method==Method::DC||
                    (method==Method::Hybrid&&patches->sharp[patches->Find(patch)]!=0);
            }

        ParallelFor(cells.size(),workers,[&](std::size_t ci)
        {
            Cell& cell=cells[ci];if(cell.componentCount==0)return;
            const int x=static_cast<int>(ci%grid.nx);
            const int y=static_cast<int>((ci/grid.nx)%grid.ny);
            const int z=static_cast<int>(ci/(grid.nx*grid.ny));
            for(std::uint8_t component=0;component<cell.componentCount;++component)
            {
                if(!cell.componentUseDC[component])continue;
                std::vector<std::pair<V3,V3>> constraints;
                V3 rootSum{};int rootCount=0;
                for(int edge=0;edge<EdgeCount;++edge)
                {
                    if(cell.edgeComponent[edge]!=component)continue;
                    const auto local=LocalEdge(x,y,z,edge);
                    const Root& root=edges.Get(local[0],local[1],local[2],local[3],grid);
                    if(!root.active)continue;
                    rootSum+=root.position;++rootCount;
                    if(root.planeCount)
                    {
                        for(std::uint8_t i=0;i<root.planeCount;++i)
                            constraints.emplace_back(root.position,root.planes[i].normal);
                    }
                    else constraints.emplace_back(root.position,root.gradient);
                }
                const V3 low=grid.Position(x,y,z),high=grid.Position(x+1,y+1,z+1);
                const V3 center=rootCount?rootSum/static_cast<double>(rootCount):(low+high)*0.5;
                cell.dual[component]=SolveQEF(constraints,center,low,high);
            }
        });

        Mesh mesh;
        std::unordered_map<std::uint64_t,std::uint32_t> dcVertexMap;
        if(method==Method::MC||method==Method::Hybrid)
        {
            std::map<std::tuple<std::int64_t,std::int64_t,std::int64_t>,std::uint32_t> positionMap;
            for(int axis=0;axis<3;++axis)for(std::size_t index=0;index<edges.roots[axis].size();++index)
            {
                const Root& root=edges.roots[axis][index];if(!root.active)continue;
                const auto q=[](double x){return static_cast<std::int64_t>(std::llround(x/PositionQuantum));};
                const V3 p=root.position;
                const auto key=std::tuple{q(p.x),q(p.y),q(p.z)};
                auto found=positionMap.find(key);
                std::uint32_t vertex=0;
                if(found==positionMap.end())
                {
                    vertex=static_cast<std::uint32_t>(mesh.vertices.size());
                    mesh.vertices.push_back({p,root.gradient});positionMap.emplace(key,vertex);
                }
                else vertex=found->second;
                edges.vertexIds[axis][index]=static_cast<int>(vertex);
            }
            // Coincident zero-at-node roots use the deterministic first edge normal. The shared
            // vertex position, matching the owner's 0.01 cm-style MC welding contract.
        }
        const auto getDCVertex=[&](int cellId,std::uint8_t component,const Root& root,int axis)->std::uint32_t
        {
            const std::uint64_t feature=SelectFeature(root,axis);
            const std::uint64_t patch=static_cast<std::uint64_t>(cellId)*EdgeCount+component;
            const std::uint64_t key=(patch<<32)^(feature&0xffffffffull);
            const auto existing=dcVertexMap.find(key);
            if(existing!=dcVertexMap.end())return existing->second;
            const V3 normal=feature?NormalForFeature(root,feature):root.gradient;
            const std::uint32_t index=static_cast<std::uint32_t>(mesh.vertices.size());
            mesh.vertices.push_back({cells[cellId].dual[component],Normalize(normal)});
            dcVertexMap.emplace(key,index);
            return index;
        };
        for(int z=0;z<grid.nz;++z)for(int y=0;y<grid.ny;++y)for(int x=0;x<grid.nx;++x)
        {
            if(!grid.IsCoreCell(x,y,z))continue;
            const std::size_t ci=grid.CellIndex(x,y,z);const Cell& cell=cells[ci];
            if(!cell.active)continue;
            std::array<int,EdgeCount> localVertices{};localVertices.fill(-1);
            for(int edge=0;edge<EdgeCount;++edge)
            {
                const auto ec=EdgeCorners[edge];
                const auto oa=CornerOffsets[ec[0]],ob=CornerOffsets[ec[1]];
                const double da=grid.density[grid.NodeIndex(x+oa[0],y+oa[1],z+oa[2])];
                const double db=grid.density[grid.NodeIndex(x+ob[0],y+ob[1],z+ob[2])];
                if((da<0.0)==(db<0.0))continue;
                const auto e=LocalEdge(x,y,z,edge);
                localVertices[edge]=edges.Vertex(e[0],e[1],e[2],e[3],grid);
            }
            for(int k=0;k<16&&MCTriTable[cell.cubeCase][k]>=0;k+=3)
            {
                const int e0=MCTriTable[cell.cubeCase][k],e1=MCTriTable[cell.cubeCase][k+1],e2=MCTriTable[cell.cubeCase][k+2];
                if(localVertices[e0]<0||localVertices[e1]<0||localVertices[e2]<0)continue;
                const std::uint8_t component=cell.edgeComponent[e0];
                if(component==0xff||cell.componentUseDC[component])continue;
                std::array<std::uint64_t,3> labels{};
                const std::array<int,3> triEdges{{e0,e1,e2}};
                for(int edgeIndex=0;edgeIndex<3;++edgeIndex)
                {
                    const auto le=LocalEdge(x,y,z,triEdges[edgeIndex]);
                    const Root& root=edges.Get(le[0],le[1],le[2],le[3],grid);
                    labels[edgeIndex]=root.feature;
                }
                std::uint64_t label=0;int bestVotes=0;
                for(std::uint64_t candidate:labels)if(candidate)
                {
                    const int votes=static_cast<int>(std::count(labels.begin(),labels.end(),candidate));
                    if(votes>bestVotes||(votes==bestVotes&&(!label||candidate<label)))
                    {label=candidate;bestVotes=votes;}
                }
                const V3 center=(mesh.vertices[localVertices[e0]].position+
                    mesh.vertices[localVertices[e1]].position+mesh.vertices[localVertices[e2]].position)/3.0;
                AddTriangle(mesh,static_cast<std::uint32_t>(localVertices[e0]),
                    static_cast<std::uint32_t>(localVertices[e1]),static_cast<std::uint32_t>(localVertices[e2]),
                    label,field.Eval(center).gradient);
            }
        }
        for(int axis=0;axis<3;++axis)
        {
            for(std::size_t index=0;index<edges.roots[axis].size();++index)
            {
                const Root& root=edges.roots[axis][index];if(!root.active)continue;
                int x=0,y=0,z=0;
                if(axis==0){x=static_cast<int>(index%grid.nx);y=static_cast<int>((index/grid.nx)%(grid.ny+1));z=static_cast<int>(index/(grid.nx*(grid.ny+1)));}
                else if(axis==1){x=static_cast<int>(index%(grid.nx+1));y=static_cast<int>((index/(grid.nx+1))%grid.ny);z=static_cast<int>(index/((grid.nx+1)*grid.ny));}
                else{x=static_cast<int>(index%(grid.nx+1));y=static_cast<int>((index/(grid.nx+1))%(grid.ny+1));z=static_cast<int>(index/((grid.nx+1)*(grid.ny+1)));}
                std::array<std::array<int,3>,4> around{};
                if(axis==0)around={{{x,y-1,z-1},{x,y,z-1},{x,y,z},{x,y-1,z}}};
                else if(axis==1)around={{{x-1,y,z-1},{x-1,y,z},{x,y,z},{x,y,z-1}}};
                else around={{{x-1,y-1,z},{x,y-1,z},{x,y,z},{x-1,y,z}}};
                std::array<int,4> incident{};bool valid=true;bool anyDC=false;
                std::array<std::uint8_t,4> incidentComponents{};
                for(int i=0;i<4;++i)
                {
                    const auto& c=around[i];
                    if(c[0]<0||c[1]<0||c[2]<0||c[0]>=grid.nx||c[1]>=grid.ny||c[2]>=grid.nz)
                    {valid=false;break;}
                    incident[i]=static_cast<int>(grid.CellIndex(c[0],c[1],c[2]));
                    const int localEdge=FindLocalEdge(c[0],c[1],c[2],axis,x,y,z);
                    if(localEdge<0){valid=false;break;}
                    incidentComponents[i]=cells[incident[i]].edgeComponent[localEdge];
                    if(incidentComponents[i]==0xff){valid=false;break;}
                    anyDC|=cells[incident[i]].componentUseDC[incidentComponents[i]];
                }
                // Own each dual quad by the lexicographically lowest incident cell. A one-cell
                // sampled halo supplies the neighbour's QEF vertex, so independently built tiles
                // emit the same complete border quad exactly once.
                const int ownerX=std::min({around[0][0],around[1][0],around[2][0],around[3][0]});
                const int ownerY=std::min({around[0][1],around[1][1],around[2][1],around[3][1]});
                const int ownerZ=std::min({around[0][2],around[1][2],around[2][2],around[3][2]});
                if(enforceQuadOwnership&&!grid.IsCoreCell(ownerX,ownerY,ownerZ))continue;
                if(!valid||!anyDC||!cells[incident[0]].componentUseDC[incidentComponents[0]])continue;
                std::array<std::uint32_t,4> v{};
                for(int i=0;i<4;++i)
                {
                    if(!cells[incident[i]].componentUseDC[incidentComponents[i]]){valid=false;break;}
                    v[i]=getDCVertex(incident[i],incidentComponents[i],root,axis);
                }
                if(!valid){++mesh.skippedDualEdges;continue;}
                if(Dot(Cross(mesh.vertices[v[1]].position-mesh.vertices[v[0]].position,
                    mesh.vertices[v[2]].position-mesh.vertices[v[0]].position),root.gradient)<0.0)
                    std::swap(v[1],v[3]);
                const std::uint64_t label=SelectFeature(root,axis);
                const std::uint32_t edgeTie=static_cast<std::uint32_t>(x+y+z+axis)&1u;
                if(edgeTie==0)
                {
                    AddTriangle(mesh,v[0],v[1],v[2],label,root.gradient);
                    AddTriangle(mesh,v[0],v[2],v[3],label,root.gradient);
                }
                else
                {
                    AddTriangle(mesh,v[0],v[1],v[3],label,root.gradient);
                    AddTriangle(mesh,v[1],v[2],v[3],label,root.gradient);
                }
            }
        }
        mesh.hash=HashMesh(mesh);
        return mesh;
    }

    struct ErrorMetrics{double mean=0,maximum=0;std::uint64_t samples=0;};
    ErrorMetrics SurfaceError(const Mesh& mesh,const Field& field)
    {
        ErrorMetrics out;double total=0;
        const auto sample=[&](V3 p)
        {
            const double error=std::abs(field.Eval(p).density);
            total+=error;out.maximum=std::max(out.maximum,error);++out.samples;
        };
        for(const Triangle& t:mesh.triangles)
        {
            const V3 a=mesh.vertices[t.index[0]].position,b=mesh.vertices[t.index[1]].position,c=mesh.vertices[t.index[2]].position;
            sample(a);sample(b);sample(c);sample((a+b)*0.5);sample((b+c)*0.5);sample((c+a)*0.5);sample((a+b+c)/3.0);
        }
        if(out.samples)out.mean=total/out.samples;
        return out;
    }

    struct GeomKey
    {
        std::int64_t x=0,y=0,z=0;
        friend bool operator<(const GeomKey& a,const GeomKey& b)
        {return std::tie(a.x,a.y,a.z)<std::tie(b.x,b.y,b.z);}
        friend bool operator==(const GeomKey&,const GeomKey&)=default;
    };
    GeomKey Key(V3 p)
    {return {static_cast<std::int64_t>(std::llround(p.x/PositionQuantum)),
        static_cast<std::int64_t>(std::llround(p.y/PositionQuantum)),
        static_cast<std::int64_t>(std::llround(p.z/PositionQuantum))};}
    struct EdgeInfo{int count=0;std::vector<std::pair<int,V3>> faces;};
    std::map<std::pair<GeomKey,GeomKey>,EdgeInfo> GeomEdges(const Mesh& mesh)
    {
        std::map<std::pair<GeomKey,GeomKey>,EdgeInfo> result;
        for(std::size_t ti=0;ti<mesh.triangles.size();++ti)
        {
            const Triangle& t=mesh.triangles[ti];
            const V3 a=mesh.vertices[t.index[0]].position,b=mesh.vertices[t.index[1]].position,c=mesh.vertices[t.index[2]].position;
            const V3 normal=Normalize(Cross(b-a,c-a));
            const std::array<std::pair<std::uint32_t,std::uint32_t>,3> edges{{{t.index[0],t.index[1]},{t.index[1],t.index[2]},{t.index[2],t.index[0]}}};
            for(const auto& e:edges)
            {
                GeomKey ka=Key(mesh.vertices[e.first].position),kb=Key(mesh.vertices[e.second].position);
                if(kb<ka)std::swap(ka,kb);
                auto& info=result[{ka,kb}];++info.count;info.faces.emplace_back(static_cast<int>(ti),normal);
            }
        }
        return result;
    }
    struct BoundaryLoop
    {
        std::vector<GeomKey> keys;
        std::vector<V3> points;
        V3 center{};
        double perimeter=0.0;
    };
    struct BoundaryLoopSet
    {
        std::vector<BoundaryLoop> loops;
        std::uint64_t boundaryEdges=0,badVertices=0;
    };
    std::pair<GeomKey,GeomKey> OrderedEdge(GeomKey a,GeomKey b)
    {if(b<a)std::swap(a,b);return {a,b};}
    using EdgeKey=std::pair<GeomKey,GeomKey>;
    BoundaryLoopSet ExtractBoundaryLoops(const Mesh& mesh)
    {
        BoundaryLoopSet result;
        const auto edges=GeomEdges(mesh);
        std::map<GeomKey,std::vector<GeomKey>> adjacency;
        for(const auto& [key,info]:edges)if(info.count==1)
        {
            ++result.boundaryEdges;
            adjacency[key.first].push_back(key.second);
            adjacency[key.second].push_back(key.first);
        }
        for(const auto& [key,neighbors]:adjacency)if(neighbors.size()!=2)++result.badVertices;
        std::set<std::pair<GeomKey,GeomKey>> visited;
        for(const auto& [start,startNeighbors]:adjacency)
        {
            for(const GeomKey first:startNeighbors)
            {
                if(visited.contains(OrderedEdge(start,first)))continue;
                BoundaryLoop loop;
                GeomKey previous=start,current=first;
                visited.insert(OrderedEdge(previous,current));
                loop.keys.push_back(start);
                bool closed=false;
                for(std::size_t guard=0;guard<=result.boundaryEdges;++guard)
                {
                    loop.keys.push_back(current);
                    const auto found=adjacency.find(current);
                    if(found==adjacency.end()||found->second.size()!=2)break;
                    const GeomKey next=found->second[0]==previous?found->second[1]:found->second[0];
                    if(next==start){closed=true;break;}
                    const auto edge=OrderedEdge(current,next);
                    if(visited.contains(edge))break;
                    visited.insert(edge);previous=current;current=next;
                }
                if(!closed||loop.keys.size()<3)continue;
                BoundaryLoop clean;
                clean.keys=std::move(loop.keys);
                for(const GeomKey& key:clean.keys)
                {
                    const V3 p{key.x*PositionQuantum,key.y*PositionQuantum,key.z*PositionQuantum};
                    clean.points.push_back(p);clean.center+=p;
                }
                clean.center=clean.center/static_cast<double>(clean.points.size());
                for(std::size_t i=0;i<clean.points.size();++i)
                    clean.perimeter+=Length(clean.points[(i+1)%clean.points.size()]-clean.points[i]);
                result.loops.push_back(std::move(clean));
            }
        }
        return result;
    }

    struct TransitionAudit
    {
        std::uint64_t pairedLoops=0,unmatchedLoops=0,cappedLoops=0,transitionTriangles=0;
        double contourMean=0.0,contourMaximum=0.0;
    };
    V3 PointAlongLoop(const BoundaryLoop& loop,std::size_t start,int direction,double fraction)
    {
        const std::size_t count=loop.points.size();
        if(!count)return {};
        double remaining=std::clamp(fraction,0.0,1.0)*loop.perimeter;
        for(std::size_t walked=0;walked<count;++walked)
        {
            const std::size_t current=direction>0?(start+walked)%count:(start+count-(walked%count))%count;
            const std::size_t next=direction>0?(current+1)%count:(current+count-1)%count;
            const V3 a=loop.points[current],b=loop.points[next];
            const double edgeLength=Length(b-a);
            if(remaining<=edgeLength||walked+1==count)
                return edgeLength>0.0?a+(b-a)*(remaining/edgeLength):a;
            remaining-=edgeLength;
        }
        return loop.points[start];
    }
    int BestLoopDirection(const BoundaryLoop& a,const BoundaryLoop& b,std::size_t start)
    {
        double scores[2]{};
        constexpr int samples=24;
        for(int dirIndex=0;dirIndex<2;++dirIndex)
        {
            const int direction=dirIndex==0?1:-1;
            for(int i=0;i<samples;++i)
            {
                const double f=static_cast<double>(i)/samples;
                scores[dirIndex]+=Dot(PointAlongLoop(a,0,1,f)-PointAlongLoop(b,start,direction,f),
                    PointAlongLoop(a,0,1,f)-PointAlongLoop(b,start,direction,f));
            }
        }
        return scores[0]<=scores[1]?1:-1;
    }
    TransitionAudit StitchBoundaryLoops(Mesh& mesh,const BoundaryLoopSet& fine,
        const BoundaryLoopSet& coarse,const Field& field)
    {
        TransitionAudit result;
        const std::size_t patchStart=mesh.triangles.size();
        result.unmatchedLoops=fine.badVertices+coarse.badVertices;
        std::map<GeomKey,std::uint32_t> vertexByKey;
        for(std::uint32_t i=0;i<mesh.vertices.size();++i)vertexByKey.try_emplace(Key(mesh.vertices[i].position),i);
        const auto loopEdges=[](const BoundaryLoopSet& loops)
        {
            std::set<EdgeKey> edges;
            for(const BoundaryLoop& loop:loops.loops)
                for(std::size_t i=0;i<loop.keys.size();++i)
                    edges.insert(OrderedEdge(loop.keys[i],loop.keys[(i+1)%loop.keys.size()]));
            return edges;
        };
        const std::set<EdgeKey> fineEdges=loopEdges(fine),coarseEdges=loopEdges(coarse);
        std::set<EdgeKey> sharedEdges;
        std::set_intersection(fineEdges.begin(),fineEdges.end(),coarseEdges.begin(),coarseEdges.end(),
            std::inserter(sharedEdges,sharedEdges.end()));
        constexpr double collarOffset=0.001;
        if(!sharedEdges.empty()&&fine.loops.size()==coarse.loops.size())
        {
            struct Chain{std::vector<GeomKey> keys;};
            const auto buildChains=[&](const BoundaryLoopSet& source)
            {
                std::vector<Chain> chains;
                for(const BoundaryLoop& loop:source.loops)
                {
                    const std::size_t n=loop.keys.size();
                    std::vector<bool> isShared(n,false);
                    for(std::size_t i=0;i<n;++i)
                        isShared[i]=sharedEdges.contains(OrderedEdge(loop.keys[i],loop.keys[(i+1)%n]));
                    if(std::none_of(isShared.begin(),isShared.end(),[](bool value){return value;}))
                    {chains.push_back({loop.keys});continue;}
                    for(std::size_t edge=0;edge<n;++edge)if(!isShared[edge]&&isShared[(edge+n-1)%n])
                    {
                        Chain chain;std::size_t current=edge;chain.keys.push_back(loop.keys[current]);
                        for(std::size_t guard=0;guard<n&&!isShared[current];++guard)
                        {
                            current=(current+1)%n;
                            chain.keys.push_back(loop.keys[current]);
                        }
                        if(chain.keys.size()>1)chains.push_back(std::move(chain));
                    }
                }
                return chains;
            };
            const std::vector<Chain> fineChains=buildChains(fine),coarseChains=buildChains(coarse);
            std::vector<bool> coarseUsed(coarseChains.size(),false);
            const std::size_t transitionStart=mesh.triangles.size();
            const auto offsetChain=[&](const Chain& chain,double offset)
            {
                std::vector<std::uint32_t> original,shifted;
                for(const GeomKey& key:chain.keys)
                {
                    const auto found=vertexByKey.find(key);
                    if(found==vertexByKey.end())return std::pair{original,shifted};
                    original.push_back(found->second);
                    Vertex vertex=mesh.vertices[found->second];vertex.position.x+=offset;
                    const GeomKey shiftedKey=Key(vertex.position);
                    const auto existing=vertexByKey.find(shiftedKey);
                    if(existing!=vertexByKey.end())shifted.push_back(existing->second);
                    else
                    {
                        const auto index=static_cast<std::uint32_t>(mesh.vertices.size());
                        mesh.vertices.push_back(vertex);vertexByKey.emplace(shiftedKey,index);shifted.push_back(index);
                    }
                }
                for(std::size_t i=0;i+1<original.size();++i)
                {
                    const auto a=original[i],b=original[i+1],sa=shifted[i],sb=shifted[i+1];
                    const V3 center=(mesh.vertices[a].position+mesh.vertices[b].position+
                        mesh.vertices[sb].position+mesh.vertices[sa].position)/4.0;
                    const FieldSample sample=field.Eval(center);
                    AddTriangle(mesh,a,b,sb,sample.primaryFeature(),sample.gradient);
                    AddTriangle(mesh,a,sb,sa,sample.primaryFeature(),sample.gradient);
                }
                return std::pair{original,shifted};
            };
            std::uint64_t matchedChains=0;
            const auto endpointDistance=[](const GeomKey& a,const GeomKey& b)
            {
                const V3 pa{a.x*PositionQuantum,a.y*PositionQuantum,a.z*PositionQuantum};
                const V3 pb{b.x*PositionQuantum,b.y*PositionQuantum,b.z*PositionQuantum};
                return Length(pa-pb);
            };
            for(const Chain& fineChain:fineChains)
            {
                std::size_t match=coarseChains.size();bool reverse=false;
                double bestEndpointScore=std::numeric_limits<double>::infinity();
                for(std::size_t j=0;j<coarseChains.size();++j)if(!coarseUsed[j])
                {
                    const Chain& candidate=coarseChains[j];
                    const double directA=endpointDistance(fineChain.keys.front(),candidate.keys.front());
                    const double directB=endpointDistance(fineChain.keys.back(),candidate.keys.back());
                    const double reverseA=endpointDistance(fineChain.keys.front(),candidate.keys.back());
                    const double reverseB=endpointDistance(fineChain.keys.back(),candidate.keys.front());
                    const double directScore=std::max(directA,directB)+1.0e-6*(directA+directB);
                    const double reverseScore=std::max(reverseA,reverseB)+1.0e-6*(reverseA+reverseB);
                    const bool candidateReverse=reverseScore<directScore;
                    const double candidateScore=std::min(directScore,reverseScore);
                    const double maxDistance=candidateReverse?std::max(reverseA,reverseB):std::max(directA,directB);
                    if(maxDistance<=LodContourErrorBudgetMeters&&candidateScore<bestEndpointScore)
                    {match=j;reverse=candidateReverse;bestEndpointScore=candidateScore;}
                }
                if(match==coarseChains.size()){++result.unmatchedLoops;continue;}
                coarseUsed[match]=true;++matchedChains;
                Chain alignedCoarse=coarseChains[match];
                if(reverse)std::reverse(alignedCoarse.keys.begin(),alignedCoarse.keys.end());
                const auto [fineOriginal,fineOffset]=offsetChain(fineChain,-collarOffset);
                const auto [coarseOriginal,coarseOffset]=offsetChain(alignedCoarse,collarOffset);
                if(fineOffset.size()<2||coarseOffset.size()<2||fineOriginal.size()!=fineOffset.size()||
                    coarseOriginal.size()!=coarseOffset.size())
                {++result.unmatchedLoops;continue;}
                for(const auto endpoint:{std::size_t{0},fineOffset.size()-1})
                {
                    const auto fa=fineOffset[endpoint];
                    const auto ca=coarseOffset[endpoint==0?0:coarseOffset.size()-1];
                    const auto originalFine=fineOriginal[endpoint];
                    const auto originalCoarse=coarseOriginal[endpoint==0?0:coarseOriginal.size()-1];
                    std::vector<std::uint32_t> capRing{originalFine,fa,ca,originalCoarse};
                    capRing.erase(std::unique(capRing.begin(),capRing.end(),[&](std::uint32_t a,std::uint32_t b)
                        {return Key(mesh.vertices[a].position)==Key(mesh.vertices[b].position);}),capRing.end());
                    if(capRing.size()>2&&Key(mesh.vertices[capRing.front()].position)==Key(mesh.vertices[capRing.back()].position))capRing.pop_back();
                    V3 capCenter{};for(const std::uint32_t index:capRing)capCenter+=mesh.vertices[index].position;
                    capCenter=capCenter/static_cast<double>(capRing.size())+V3{0.0,0.001,0.001};
                    Vertex capVertex;capVertex.position=capCenter;
                    capVertex.normal=Normalize(mesh.vertices[fa].normal+mesh.vertices[ca].normal);
                    const GeomKey capKey=Key(capVertex.position);
                    auto capFound=vertexByKey.find(capKey);
                    std::uint32_t cap=0;
                    if(capFound!=vertexByKey.end())cap=capFound->second;
                    else
                    {
                        cap=static_cast<std::uint32_t>(mesh.vertices.size());
                        mesh.vertices.push_back(capVertex);vertexByKey.emplace(capKey,cap);
                    }
                    for(std::size_t i=0;i<capRing.size();++i)
                    {
                        const std::uint32_t a=capRing[i],b=capRing[(i+1)%capRing.size()];
                        const V3 center=(mesh.vertices[cap].position+mesh.vertices[a].position+
                            mesh.vertices[b].position)/3.0;
                        const FieldSample sample=field.Eval(center);
                        AddTriangle(mesh,cap,a,b,sample.primaryFeature(),sample.gradient);
                    }
                }
                const std::size_t na=fineOffset.size(),nb=coarseOffset.size();
                std::vector<double> cumulativeA(na,0.0),cumulativeB(nb,0.0);
                for(std::size_t i=1;i<na;++i)cumulativeA[i]=cumulativeA[i-1]+Length(mesh.vertices[fineOffset[i]].position-mesh.vertices[fineOffset[i-1]].position);
                for(std::size_t i=1;i<nb;++i)cumulativeB[i]=cumulativeB[i-1]+Length(mesh.vertices[coarseOffset[i]].position-mesh.vertices[coarseOffset[i-1]].position);
                const double perimA=cumulativeA.back(),perimB=cumulativeB.back();
                if(perimA<=0.0||perimB<=0.0){++result.unmatchedLoops;continue;}
                std::size_t i=0,j=0;
                while(i+1<na||j+1<nb)
                {
                    const double nextA=i+1<na?cumulativeA[i+1]/perimA:std::numeric_limits<double>::infinity();
                    const double nextB=j+1<nb?cumulativeB[j+1]/perimB:std::numeric_limits<double>::infinity();
                    if(nextA<=nextB)
                    {
                        const V3 center=(mesh.vertices[fineOffset[i]].position+mesh.vertices[fineOffset[i+1]].position+
                            mesh.vertices[coarseOffset[j]].position)/3.0;
                        const FieldSample sample=field.Eval(center);
                        AddTriangle(mesh,fineOffset[i],fineOffset[i+1],coarseOffset[j],sample.primaryFeature(),sample.gradient);++i;
                    }
                    else
                    {
                        const V3 center=(mesh.vertices[fineOffset[i]].position+mesh.vertices[coarseOffset[j+1]].position+
                            mesh.vertices[coarseOffset[j]].position)/3.0;
                        const FieldSample sample=field.Eval(center);
                        AddTriangle(mesh,fineOffset[i],coarseOffset[j+1],coarseOffset[j],sample.primaryFeature(),sample.gradient);++j;
                    }
                }
            }
            for(bool used:coarseUsed)if(!used)++result.unmatchedLoops;
            result.pairedLoops=matchedChains;
            result.transitionTriangles=mesh.triangles.size()-transitionStart;
            return result;
        }
        // Give both sides a one-millimetre transition skin before zipping. This avoids degenerate
        // triangles when otherwise independent contours share collinear runs.
        const auto offsetLoops=[&](const BoundaryLoopSet& source,double offset)
        {
            std::vector<BoundaryLoop> output;output.reserve(source.loops.size());
            for(const BoundaryLoop& input:source.loops)
            {
                BoundaryLoop loop;bool valid=true;
                for(const GeomKey& key:input.keys)
                {
                    const auto found=vertexByKey.find(key);
                    if(found==vertexByKey.end()){valid=false;break;}
                    Vertex vertex=mesh.vertices[found->second];vertex.position.x+=offset;
                    const std::uint32_t index=static_cast<std::uint32_t>(mesh.vertices.size());
                    mesh.vertices.push_back(vertex);vertexByKey.try_emplace(Key(vertex.position),index);
                    loop.keys.push_back(Key(vertex.position));loop.points.push_back(vertex.position);loop.center+=vertex.position;
                }
                if(!valid||loop.points.size()<3){++result.unmatchedLoops;continue;}
                loop.center=loop.center/static_cast<double>(loop.points.size());
                for(std::size_t i=0;i<loop.points.size();++i)
                {
                    const std::size_t next=(i+1)%loop.points.size();
                    loop.perimeter+=Length(loop.points[next]-loop.points[i]);
                    const auto originalA=vertexByKey.find(input.keys[i]);
                    const auto originalB=vertexByKey.find(input.keys[next]);
                    const auto shiftedA=vertexByKey.find(loop.keys[i]);
                    const auto shiftedB=vertexByKey.find(loop.keys[next]);
                    if(originalA==vertexByKey.end()||originalB==vertexByKey.end()||
                        shiftedA==vertexByKey.end()||shiftedB==vertexByKey.end())
                    {++result.unmatchedLoops;continue;}
                    const auto a=originalA->second,b=originalB->second,sa=shiftedA->second,sb=shiftedB->second;
                    const V3 center=(mesh.vertices[a].position+mesh.vertices[b].position+
                        mesh.vertices[sb].position+mesh.vertices[sa].position)/4.0;
                    const FieldSample sample=field.Eval(center);
                    AddTriangle(mesh,a,b,sb,sample.primaryFeature(),sample.gradient);
                    AddTriangle(mesh,a,sb,sa,sample.primaryFeature(),sample.gradient);
                }
                output.push_back(std::move(loop));
            }
            return output;
        };
        const std::vector<BoundaryLoop> fineOffset=offsetLoops(fine,-collarOffset);
        const std::vector<BoundaryLoop> coarseOffset=offsetLoops(coarse,collarOffset);
        const auto capLoop=[&](const BoundaryLoop& loop)
        {
            if(loop.keys.size()<3)return false;
            V3 center{};for(const V3 point:loop.points)center+=point;
            center=center/static_cast<double>(loop.points.size())+V3{0.001,0.001,0.001};
            Vertex centerVertex;centerVertex.position=center;centerVertex.normal=Normalize(field.Eval(center).gradient);
            const GeomKey centerKey=Key(center);
            auto centerIt=vertexByKey.find(centerKey);std::uint32_t centerIndex=0;
            if(centerIt!=vertexByKey.end())centerIndex=centerIt->second;
            else
            {
                centerIndex=static_cast<std::uint32_t>(mesh.vertices.size());
                mesh.vertices.push_back(centerVertex);vertexByKey.emplace(centerKey,centerIndex);
            }
            const std::size_t before=mesh.triangles.size();
            for(std::size_t i=0;i<loop.keys.size();++i)
            {
                const auto a=vertexByKey.find(loop.keys[i]),b=vertexByKey.find(loop.keys[(i+1)%loop.keys.size()]);
                if(a==vertexByKey.end()||b==vertexByKey.end())continue;
                const V3 at=(mesh.vertices[centerIndex].position+mesh.vertices[a->second].position+
                    mesh.vertices[b->second].position)/3.0;
                const FieldSample sample=field.Eval(at);
                AddTriangle(mesh,centerIndex,a->second,b->second,sample.primaryFeature(),sample.gradient);
            }
            if(mesh.triangles.size()==before)return false;
            ++result.cappedLoops;return true;
        };
        std::vector<bool> used(coarseOffset.size(),false);
        double totalResidual=0.0;std::uint64_t residualSamples=0;
        for(const BoundaryLoop& a:fineOffset)
        {
            std::size_t match=coarseOffset.size();double best=std::numeric_limits<double>::infinity();
            for(std::size_t j=0;j<coarseOffset.size();++j)if(!used[j])
            {
                const double d=Length(a.center-coarseOffset[j].center);
                if(d<best){best=d;match=j;}
            }
            if(match==coarseOffset.size()){if(!capLoop(a))++result.unmatchedLoops;continue;}
            used[match]=true;++result.pairedLoops;
            const BoundaryLoop& b=coarseOffset[match];
            if(a.keys.size()<3||b.keys.size()<3){++result.unmatchedLoops;continue;}
            std::size_t bStart=0;double nearest=std::numeric_limits<double>::infinity();
            for(std::size_t j=0;j<b.points.size();++j)
            {
                const double d=Length(a.points[0]-b.points[j]);
                if(d<nearest){nearest=d;bStart=j;}
            }
            const int bDirection=BestLoopDirection(a,b,bStart);
            std::vector<std::uint32_t> aIndices,bIndices;
            for(const GeomKey& key:a.keys)
            {
                const auto it=vertexByKey.find(key);
                if(it==vertexByKey.end()){++result.unmatchedLoops;break;}
                aIndices.push_back(it->second);
            }
            for(std::size_t j=0;j<b.keys.size();++j)
            {
                const std::size_t bi=bDirection>0?(bStart+j)%b.keys.size():(bStart+b.keys.size()-(j%b.keys.size()))%b.keys.size();
                const auto it=vertexByKey.find(b.keys[bi]);
                if(it==vertexByKey.end()){++result.unmatchedLoops;break;}
                bIndices.push_back(it->second);
            }
            if(aIndices.size()!=a.keys.size()||bIndices.size()!=b.keys.size())continue;
            const std::size_t na=aIndices.size(),nb=bIndices.size();
            std::vector<double> cumulativeA(na+1,0.0),cumulativeB(nb+1,0.0);
            for(std::size_t i=0;i<na;++i)cumulativeA[i+1]=cumulativeA[i]+Length(mesh.vertices[aIndices[(i+1)%na]].position-mesh.vertices[aIndices[i]].position);
            for(std::size_t i=0;i<nb;++i)cumulativeB[i+1]=cumulativeB[i]+Length(mesh.vertices[bIndices[(i+1)%nb]].position-mesh.vertices[bIndices[i]].position);
            const double perimA=cumulativeA.back(),perimB=cumulativeB.back();
            if(perimA<=0.0||perimB<=0.0){++result.unmatchedLoops;continue;}
            std::size_t i=0,j=0;
            while(i<na||j<nb)
            {
                const double nextA=i<na?cumulativeA[i+1]/perimA:std::numeric_limits<double>::infinity();
                const double nextB=j<nb?cumulativeB[j+1]/perimB:std::numeric_limits<double>::infinity();
                if(nextA<=nextB)
                {
                    const auto a0=aIndices[i%na],a1=aIndices[(i+1)%na],b0=bIndices[j%nb];
                    const V3 center=(mesh.vertices[a0].position+mesh.vertices[a1].position+mesh.vertices[b0].position)/3.0;
                    const FieldSample sample=field.Eval(center);
                    AddTriangle(mesh,a0,a1,b0,sample.primaryFeature(),sample.gradient);++i;
                }
                else
                {
                    const auto a0=aIndices[i%na],b1=bIndices[(j+1)%nb],b0=bIndices[j%nb];
                    const V3 center=(mesh.vertices[a0].position+mesh.vertices[b1].position+mesh.vertices[b0].position)/3.0;
                    const FieldSample sample=field.Eval(center);
                    AddTriangle(mesh,a0,b1,b0,sample.primaryFeature(),sample.gradient);++j;
                }
            }
            const std::size_t sampleCount=std::max(na,nb);
            for(std::size_t s=0;s<sampleCount;++s)
            {
                const double f=static_cast<double>(s)/sampleCount;
                const V3 pa=PointAlongLoop(a,0,1,f),pb=PointAlongLoop(b,bStart,bDirection,f);
                const double d=Length(pa-pb);totalResidual+=d;result.contourMaximum=std::max(result.contourMaximum,d);++residualSamples;
            }
        }
        for(std::size_t j=0;j<used.size();++j)if(!used[j]&&!capLoop(coarseOffset[j]))++result.unmatchedLoops;
        if(residualSamples)result.contourMean=totalResidual/residualSamples;
        result.transitionTriangles=mesh.triangles.size()-patchStart;
        return result;
    }
    struct TopologyMetrics
    {
        std::uint64_t boundaryEdges=0,nonManifoldEdges=0,creaseEdgePairs=0,splitPairs=0;
    };
    TopologyMetrics AnalyzeTopology(const Mesh& mesh)
    {
        TopologyMetrics out;
        const auto edges=GeomEdges(mesh);
        for(const auto& [key,info]:edges)
        {
            if(info.count==1)++out.boundaryEdges;
            else if(info.count!=2)++out.nonManifoldEdges;
            if(info.faces.size()==2)
            {
                const Triangle& a=mesh.triangles[info.faces[0].first];
                const Triangle& b=mesh.triangles[info.faces[1].first];
                if(a.feature&&b.feature&&a.feature!=b.feature&&
                    Dot(Normalize(info.faces[0].second),Normalize(info.faces[1].second))<NormalSplitCos)
                {
                    ++out.creaseEdgePairs;
                    const auto avgNormal=[&](const Triangle& t)
                    {return Normalize(mesh.vertices[t.index[0]].normal+mesh.vertices[t.index[1]].normal+mesh.vertices[t.index[2]].normal);};
                    if(Dot(avgNormal(a),avgNormal(b))<NormalSplitCos)++out.splitPairs;
                }
            }
        }
        return out;
    }
    TopologyMetrics AnalyzeTopologyInXBand(const Mesh& mesh,double plane,double halfWidth)
    {
        TopologyMetrics out;
        for(const auto& [edge,info]:GeomEdges(mesh))
        {
            const double ax=edge.first.x*PositionQuantum,bx=edge.second.x*PositionQuantum;
            if(std::abs(ax-plane)>halfWidth||std::abs(bx-plane)>halfWidth)continue;
            if(info.count==1)++out.boundaryEdges;
            else if(info.count!=2)++out.nonManifoldEdges;
        }
        return out;
    }
    struct SeamCapAudit
    {
        std::uint64_t cappedLoops=0,unmatchedLoops=0,triangles=0;
    };
    SeamCapAudit CapClosedBoundaryLoopsInXBand(Mesh& mesh,const Field& field,double plane,double halfWidth)
    {
        SeamCapAudit result;
        std::map<GeomKey,std::vector<GeomKey>> adjacency;
        std::set<EdgeKey> boundaryEdges;
        std::map<GeomKey,std::uint32_t> vertexByKey;
        for(std::uint32_t i=0;i<mesh.vertices.size();++i)vertexByKey.try_emplace(Key(mesh.vertices[i].position),i);
        for(const auto& [edge,info]:GeomEdges(mesh))if(info.count==1&&
            std::abs(edge.first.x*PositionQuantum-plane)<=halfWidth&&
            std::abs(edge.second.x*PositionQuantum-plane)<=halfWidth)
        {
            boundaryEdges.insert(edge);
            adjacency[edge.first].push_back(edge.second);
            adjacency[edge.second].push_back(edge.first);
        }
        std::set<EdgeKey> visited;
        for(const EdgeKey& firstEdge:boundaryEdges)
        {
            if(visited.contains(firstEdge))continue;
            std::vector<GeomKey> loop{firstEdge.first};
            GeomKey previous=firstEdge.first,current=firstEdge.second;
            visited.insert(firstEdge);
            bool closed=false,valid=true;
            for(std::size_t guard=0;guard<=boundaryEdges.size();++guard)
            {
                if(current==loop.front()){closed=true;break;}
                loop.push_back(current);
                const auto found=adjacency.find(current);
                if(found==adjacency.end()||found->second.size()!=2){valid=false;break;}
                const GeomKey next=found->second[0]==previous?found->second[1]:found->second[0];
                const EdgeKey edge=OrderedEdge(current,next);
                if(visited.contains(edge)){valid=false;break;}
                visited.insert(edge);previous=current;current=next;
            }
            if(!closed||!valid||loop.size()<3)
            {++result.unmatchedLoops;continue;}
            V3 center{};bool haveVertices=true;
            for(const GeomKey& key:loop)
            {
                const auto found=vertexByKey.find(key);
                if(found==vertexByKey.end()){haveVertices=false;break;}
                center+=mesh.vertices[found->second].position;
            }
            if(!haveVertices){++result.unmatchedLoops;continue;}
            center=center/static_cast<double>(loop.size());
            center+=Normalize(field.Eval(center).gradient)*0.0005;
            Vertex centerVertex;centerVertex.position=center;centerVertex.normal=Normalize(field.Eval(center).gradient);
            const GeomKey centerKey=Key(center);
            const std::size_t verticesBefore=mesh.vertices.size();
            auto centerIt=vertexByKey.find(centerKey);std::uint32_t centerIndex=0;
            if(centerIt!=vertexByKey.end())centerIndex=centerIt->second;
            else
            {
                centerIndex=static_cast<std::uint32_t>(mesh.vertices.size());
                mesh.vertices.push_back(centerVertex);vertexByKey.emplace(centerKey,centerIndex);
            }
            const std::size_t before=mesh.triangles.size();
            for(std::size_t i=0;i<loop.size();++i)
            {
                const auto a=vertexByKey.find(loop[i]),b=vertexByKey.find(loop[(i+1)%loop.size()]);
                if(a==vertexByKey.end()||b==vertexByKey.end())continue;
                const V3 samplePoint=(mesh.vertices[centerIndex].position+
                    mesh.vertices[a->second].position+mesh.vertices[b->second].position)/3.0;
                const FieldSample sample=field.Eval(samplePoint);
                AddTriangle(mesh,centerIndex,a->second,b->second,sample.primaryFeature(),sample.gradient);
            }
            if(mesh.triangles.size()==before)
            {
                if(mesh.vertices.size()>verticesBefore)
                {mesh.vertices.resize(verticesBefore);vertexByKey.erase(centerKey);}
                ++result.unmatchedLoops;continue;
            }
            Mesh capOnly;capOnly.vertices=mesh.vertices;
            capOnly.triangles.assign(mesh.triangles.begin()+static_cast<std::ptrdiff_t>(before),mesh.triangles.end());
            if(SurfaceError(capOnly,field).maximum>LodContourErrorBudgetMeters)
            {
                mesh.triangles.resize(before);
                if(mesh.vertices.size()>verticesBefore)
                {mesh.vertices.resize(verticesBefore);vertexByKey.erase(centerKey);}
                ++result.unmatchedLoops;continue;
            }
            ++result.cappedLoops;result.triangles+=mesh.triangles.size()-before;
        }
        return result;
    }

    std::uint64_t FnvBytes(const std::uint8_t* p,std::size_t n)
    {
        std::uint64_t h=14695981039346656037ull;
        for(std::size_t i=0;i<n;++i){h^=p[i];h*=1099511628211ull;}
        return h;
    }

    void WriteObj(const std::filesystem::path& path,const Mesh& mesh,double sideFilter=0.0,int side=0)
    {
        std::ofstream out(path,std::ios::binary);
        out<<std::setprecision(9);
        for(const Vertex& v:mesh.vertices)out<<"v "<<v.position.x<<' '<<v.position.y<<' '<<v.position.z<<"\n";
        for(const Vertex& v:mesh.vertices)out<<"vn "<<v.normal.x<<' '<<v.normal.y<<' '<<v.normal.z<<"\n";
        for(const Triangle& t:mesh.triangles)
        {
            const V3 center=(mesh.vertices[t.index[0]].position+mesh.vertices[t.index[1]].position+mesh.vertices[t.index[2]].position)/3.0;
            if(side<0&&center.x>sideFilter)continue;
            if(side>0&&center.x<sideFilter)continue;
            out<<"f "<<t.index[0]+1<<"//"<<t.index[0]+1<<' '
                <<t.index[1]+1<<"//"<<t.index[1]+1<<' '<<t.index[2]+1<<"//"<<t.index[2]+1<<"\n";
        }
    }

    struct Segment{V3 a{},b{};};
    std::vector<Segment> SeamSegments(const Mesh& mesh,double xPlane,int side)
    {
        std::vector<Segment> segments;
        for(const Triangle& triangle:mesh.triangles)
        {
            std::array<V3,3> input{{mesh.vertices[triangle.index[0]].position,
                mesh.vertices[triangle.index[1]].position,mesh.vertices[triangle.index[2]].position}};
            std::vector<V3> output;
            for(int i=0;i<3;++i)
            {
                const V3 a=input[i],b=input[(i+1)%3];
                const bool ina=side<0?a.x<=xPlane+PositionQuantum:a.x>=xPlane-PositionQuantum;
                const bool inb=side<0?b.x<=xPlane+PositionQuantum:b.x>=xPlane-PositionQuantum;
                if(ina)output.push_back(a);
                if(ina!=inb)
                {
                    const double t=(xPlane-a.x)/(b.x-a.x);
                    output.push_back(Quantize(a+(b-a)*t));
                }
            }
            std::vector<V3> onPlane;
            for(V3 p:output)if(std::abs(p.x-xPlane)<=2.0*PositionQuantum)
            {
                bool found=false;for(V3 q:onPlane)if(Length(p-q)<PositionQuantum){found=true;break;}
                if(!found)onPlane.push_back(p);
            }
            if(onPlane.size()==2&&Length(onPlane[0]-onPlane[1])>PositionQuantum)
                segments.push_back({onPlane[0],onPlane[1]});
        }
        return segments;
    }
    double PointSegmentDistance(V3 p,const Segment& s)
    {
        const V3 d=s.b-s.a;const double den=Dot(d,d);
        const double t=den>0?std::clamp(Dot(p-s.a,d)/den,0.0,1.0):0.0;
        return Length(p-(s.a+d*t));
    }
    struct SeamMetric{double mean=0,maximum=0;std::uint64_t probes=0,unmatched=0;};
    SeamMetric CompareSegments(const std::vector<Segment>& a,const std::vector<Segment>& b)
    {
        SeamMetric out;double total=0;
        const auto probe=[&](const std::vector<Segment>& source,const std::vector<Segment>& target)
        {
            for(const Segment& s:source)for(double t:{0.0,0.25,0.5,0.75,1.0})
            {
                const V3 p=s.a+(s.b-s.a)*t;double nearest=std::numeric_limits<double>::infinity();
                for(const Segment& candidate:target)nearest=std::min(nearest,PointSegmentDistance(p,candidate));
                if(!std::isfinite(nearest)){++out.unmatched;++out.probes;continue;}
                total+=nearest;out.maximum=std::max(out.maximum,nearest);++out.probes;
            }
        };
        probe(a,b);probe(b,a);
        if(out.unmatched){out.mean=-1.0;out.maximum=-1.0;}
        else if(out.probes)out.mean=total/out.probes;
        return out;
    }
    bool HasIsoContourOnX(const Field& field,double xPlane,double low,double high,double step)
    {
        const int cells=static_cast<int>(std::llround((high-low)/step));
        for(int y=0;y<cells;++y)for(int z=0;z<cells;++z)
        {
            const double x0=field.Eval({xPlane,low+y*step,low+z*step}).density;
            const double x1=field.Eval({xPlane,low+(y+1)*step,low+z*step}).density;
            const double x2=field.Eval({xPlane,low+y*step,low+(z+1)*step}).density;
            if((x0<0.0)!=(x1<0.0)||(x0<0.0)!=(x2<0.0))return true;
        }
        return false;
    }
    Mesh MergeMeshes(const Mesh& a,const Mesh& b)
    {
        Mesh result;result.vertices=a.vertices;result.triangles=a.triangles;
        const std::uint32_t offset=static_cast<std::uint32_t>(result.vertices.size());
        result.vertices.insert(result.vertices.end(),b.vertices.begin(),b.vertices.end());
        for(Triangle triangle:b.triangles)
        {
            for(std::uint32_t& index:triangle.index)index+=offset;
            result.triangles.push_back(triangle);
        }
        result.hash=HashMesh(result);return result;
    }
    Mesh ClipMeshAtX(const Mesh& source,const Field& field,bool keepLeft)
    {
        Mesh result;
        std::map<GeomKey,std::uint32_t> vertexByKey;
        const auto addVertex=[&](Vertex vertex)
        {
            if(std::abs(vertex.position.x)<PositionQuantum)vertex.position.x=0.0;
            const GeomKey key=Key(vertex.position);
            const auto found=vertexByKey.find(key);
            if(found!=vertexByKey.end())return found->second;
            const auto index=static_cast<std::uint32_t>(result.vertices.size());
            result.vertices.push_back(vertex);vertexByKey.emplace(key,index);return index;
        };
        for(const Triangle& triangle:source.triangles)
        {
            std::vector<Vertex> polygon;polygon.reserve(4);
            for(int edge=0;edge<3;++edge)
            {
                const Vertex a=source.vertices[triangle.index[edge]];
                const Vertex b=source.vertices[triangle.index[(edge+1)%3]];
                const bool insideA=keepLeft?a.position.x<=0.0:a.position.x>=0.0;
                const bool insideB=keepLeft?b.position.x<=0.0:b.position.x>=0.0;
                if(insideA)polygon.push_back(a);
                if(insideA!=insideB)
                {
                    const double t=std::clamp(-a.position.x/(b.position.x-a.position.x),0.0,1.0);
                    Vertex cut;
                    cut.position=a.position+(b.position-a.position)*t;cut.position.x=0.0;
                    cut.normal=Normalize(a.normal+(b.normal-a.normal)*t);
                    polygon.push_back(cut);
                }
            }
            std::vector<std::uint32_t> indices;indices.reserve(polygon.size());
            for(const Vertex& vertex:polygon)
            {
                const auto index=addVertex(vertex);
                if(indices.empty()||indices.back()!=index)indices.push_back(index);
            }
            if(indices.size()>2&&indices.front()==indices.back())indices.pop_back();
            if(indices.size()<3)continue;
            for(std::size_t i=1;i+1<indices.size();++i)
            {
                const V3 center=(result.vertices[indices[0]].position+
                    result.vertices[indices[i]].position+result.vertices[indices[i+1]].position)/3.0;
                AddTriangle(result,indices[0],indices[i],indices[i+1],triangle.feature,field.Eval(center).gradient);
            }
        }
        result.hash=HashMesh(result);return result;
    }
    std::vector<Segment> BoundarySegmentsAtX(const Mesh& mesh,double xPlane)
    {
        std::vector<Segment> segments;
        for(const auto& [key,info]:GeomEdges(mesh))
        {
            if(info.count!=1)continue;
            const V3 a{key.first.x*PositionQuantum,key.first.y*PositionQuantum,key.first.z*PositionQuantum};
            const V3 b{key.second.x*PositionQuantum,key.second.y*PositionQuantum,key.second.z*PositionQuantum};
            if(std::abs(a.x-xPlane)<=PositionQuantum&&std::abs(b.x-xPlane)<=PositionQuantum)
                segments.push_back({a,b});
        }
        return segments;
    }
    struct BoundaryRepairAudit
    {
        std::uint64_t triangles=0,unmatchedEdges=0;
        double maximumResidual=0.0;
    };
    std::uint64_t ConformSeamBoundaryEdges(Mesh& a,Mesh& b,double xPlane,const Field& field)
    {
        std::map<GeomKey,V3> uniquePoints;
        for(const Mesh* mesh:{&a,&b})for(const Vertex& vertex:mesh->vertices)
            if(std::abs(vertex.position.x-xPlane)<=PositionQuantum)
                uniquePoints.try_emplace(Key(vertex.position),vertex.position);
        std::vector<V3> points;points.reserve(uniquePoints.size());
        for(const auto& [key,p]:uniquePoints)points.push_back(p);
        const auto conform=[&](Mesh& source)
        {
            Mesh output;output.vertices=source.vertices;
            std::map<GeomKey,std::uint32_t> vertexByKey;
            for(std::uint32_t i=0;i<output.vertices.size();++i)vertexByKey.try_emplace(Key(output.vertices[i].position),i);
            const auto addVertex=[&](Vertex vertex)
            {
                const GeomKey key=Key(vertex.position);
                const auto found=vertexByKey.find(key);
                if(found!=vertexByKey.end())return found->second;
                const auto index=static_cast<std::uint32_t>(output.vertices.size());
                output.vertices.push_back(vertex);vertexByKey.emplace(key,index);return index;
            };
            std::uint64_t splitTriangles=0;
            for(const Triangle& triangle:source.triangles)
            {
                std::vector<std::uint32_t> polygon;polygon.reserve(8);bool split=false;
                for(int edge=0;edge<3;++edge)
                {
                    const auto aIndex=triangle.index[edge],bIndex=triangle.index[(edge+1)%3];
                    const V3 pa=source.vertices[aIndex].position,pb=source.vertices[bIndex].position;
                    polygon.push_back(addVertex(source.vertices[aIndex]));
                    if(std::abs(pa.x-xPlane)>PositionQuantum||std::abs(pb.x-xPlane)>PositionQuantum)continue;
                    const V3 direction=pb-pa;const double denominator=Dot(direction,direction);
                    if(denominator<=PositionQuantum*PositionQuantum)continue;
                    std::vector<std::pair<double,V3>> interior;
                    for(const V3 q:points)
                    {
                        if(q.y<std::min(pa.y,pb.y)-PositionQuantum||q.y>std::max(pa.y,pb.y)+PositionQuantum||
                            q.z<std::min(pa.z,pb.z)-PositionQuantum||q.z>std::max(pa.z,pb.z)+PositionQuantum)continue;
                        const double t=Dot(q-pa,direction)/denominator;
                        if(t<=1.0e-8||t>=1.0-1.0e-8)continue;
                        const V3 projected=pa+direction*t;
                        if(Length(projected-q)>PositionQuantum)continue;
                        interior.emplace_back(t,q);
                    }
                    std::sort(interior.begin(),interior.end(),[](const auto& x,const auto& y){return x.first<y.first;});
                    double previous=-1.0;
                    for(const auto& [t,p]:interior)
                    {
                        if(previous>=0.0&&std::abs(t-previous)<1.0e-8)continue;
                        Vertex vertex;vertex.position=p;vertex.normal=Normalize(source.vertices[aIndex].normal+
                            (source.vertices[bIndex].normal-source.vertices[aIndex].normal)*t);
                        polygon.push_back(addVertex(vertex));previous=t;split=true;
                    }
                }
                if(!split){output.triangles.push_back(triangle);continue;}
                ++splitTriangles;
                const V3 center=(source.vertices[triangle.index[0]].position+
                    source.vertices[triangle.index[1]].position+source.vertices[triangle.index[2]].position)/3.0;
                Vertex centerVertex;centerVertex.position=center;
                centerVertex.normal=Normalize(source.vertices[triangle.index[0]].normal+
                    source.vertices[triangle.index[1]].normal+source.vertices[triangle.index[2]].normal);
                const auto centerIndex=addVertex(centerVertex);
                for(std::size_t i=0;i<polygon.size();++i)
                {
                    const auto p0=polygon[i],p1=polygon[(i+1)%polygon.size()];
                    const V3 triCenter=(output.vertices[centerIndex].position+output.vertices[p0].position+
                        output.vertices[p1].position)/3.0;
                    const FieldSample sample=field.Eval(triCenter);
                    AddTriangle(output,centerIndex,p0,p1,triangle.feature,sample.gradient);
                }
            }
            output.hash=HashMesh(output);source=std::move(output);return splitTriangles;
        };
        return conform(a)+conform(b);
    }
    bool HasTJunctionAtX(const Mesh& a,const Mesh& b,double xPlane)
    {
        std::map<GeomKey,V3> seamPoints;
        for(const Mesh* mesh:{&a,&b})for(const Vertex& vertex:mesh->vertices)
            if(std::abs(vertex.position.x-xPlane)<=PositionQuantum)
                seamPoints.try_emplace(Key(vertex.position),vertex.position);
        for(const Mesh* mesh:{&a,&b})
            for(const auto& [edge,info]:GeomEdges(*mesh))
            {
                if(info.count!=1||std::abs(edge.first.x*PositionQuantum-xPlane)>PositionQuantum||
                    std::abs(edge.second.x*PositionQuantum-xPlane)>PositionQuantum)continue;
                const V3 p{edge.first.x*PositionQuantum,edge.first.y*PositionQuantum,edge.first.z*PositionQuantum};
                const V3 q{edge.second.x*PositionQuantum,edge.second.y*PositionQuantum,edge.second.z*PositionQuantum};
                const V3 d=q-p;const double denominator=Dot(d,d);
                if(denominator<=PositionQuantum*PositionQuantum)continue;
                for(const auto& [key,point]:seamPoints)
                {
                    const double t=Dot(point-p,d)/denominator;
                    if(t<=1.0e-8||t>=1.0-1.0e-8)continue;
                    if(Length(p+d*t-point)<=PositionQuantum)return true;
                }
            }
        return false;
    }
    BoundaryRepairAudit CloseOwnedSeamLoops(Mesh& owner,const Mesh& neighbor,double xPlane,const Field& field)
    {
        BoundaryRepairAudit result;
        const Mesh merged=MergeMeshes(owner,neighbor);
        const auto edges=GeomEdges(merged);
        std::map<GeomKey,std::vector<GeomKey>> adjacency;
        for(const auto& [edge,info]:edges)if(info.count==1&&
            std::abs(edge.first.x*PositionQuantum-xPlane)<=PositionQuantum&&
            std::abs(edge.second.x*PositionQuantum-xPlane)<=PositionQuantum)
        {
            adjacency[edge.first].push_back(edge.second);adjacency[edge.second].push_back(edge.first);
        }
        for(const auto& [key,neighbors]:adjacency)if(neighbors.size()!=2)++result.unmatchedEdges;
        std::set<EdgeKey> visited;
        std::map<GeomKey,std::uint32_t> ownerVertices;
        for(std::uint32_t i=0;i<owner.vertices.size();++i)ownerVertices.try_emplace(Key(owner.vertices[i].position),i);
        const std::size_t triangleStart=owner.triangles.size();
        for(const auto& [start,startNeighbors]:adjacency)for(const GeomKey first:startNeighbors)
        {
            if(visited.contains(OrderedEdge(start,first)))continue;
            std::vector<GeomKey> loop{start};GeomKey previous=start,current=first;
            visited.insert(OrderedEdge(previous,current));
            bool closed=false;
            for(std::size_t guard=0;guard<=adjacency.size();++guard)
            {
                loop.push_back(current);
                const auto found=adjacency.find(current);
                if(found==adjacency.end()||found->second.size()!=2)break;
                const GeomKey next=found->second[0]==previous?found->second[1]:found->second[0];
                if(next==start){closed=true;break;}
                const EdgeKey edge=OrderedEdge(current,next);
                if(visited.contains(edge))break;
                visited.insert(edge);previous=current;current=next;
            }
            if(!closed||loop.size()<3){++result.unmatchedEdges;continue;}
            std::vector<std::uint32_t> indices;indices.reserve(loop.size());
            for(const GeomKey& key:loop)
            {
                const auto existing=ownerVertices.find(key);
                if(existing!=ownerVertices.end()){indices.push_back(existing->second);continue;}
                const V3 p{key.x*PositionQuantum,key.y*PositionQuantum,key.z*PositionQuantum};
                const FieldSample sample=field.Eval(p);
                const auto index=static_cast<std::uint32_t>(owner.vertices.size());
                owner.vertices.push_back({p,Normalize(sample.gradient)});ownerVertices.emplace(key,index);indices.push_back(index);
            }
            for(std::size_t i=1;i+1<indices.size();++i)
            {
                const V3 center=(owner.vertices[indices[0]].position+owner.vertices[indices[i]].position+
                    owner.vertices[indices[i+1]].position)/3.0;
                const FieldSample sample=field.Eval(center);
                AddTriangle(owner,indices[0],indices[i],indices[i+1],sample.primaryFeature(),sample.gradient);
            }
        }
        result.triangles=owner.triangles.size()-triangleStart;
        Mesh added;added.vertices=owner.vertices;
        added.triangles.assign(owner.triangles.begin()+static_cast<std::ptrdiff_t>(triangleStart),owner.triangles.end());
        result.maximumResidual=SurfaceError(added,field).maximum;
        return result;
    }
    BoundaryRepairAudit ConformOpenTJunctions(Mesh& mesh,const Field& field)
    {
        BoundaryRepairAudit result;
        for(int pass=0;pass<2;++pass)
        {
            const auto edgeMap=GeomEdges(mesh);
            std::set<EdgeKey> boundaryEdges;
            for(const auto& [edge,info]:edgeMap)if(info.count==1)boundaryEdges.insert(edge);
            if(boundaryEdges.empty())break;
            std::map<GeomKey,V3> pointMap;
            for(const Vertex& vertex:mesh.vertices)pointMap.try_emplace(Key(vertex.position),vertex.position);
            std::vector<V3> points;points.reserve(pointMap.size());
            for(const auto& [key,p]:pointMap)points.push_back(p);
            Mesh output;output.vertices=mesh.vertices;
            std::map<GeomKey,std::uint32_t> vertexByKey;
            for(std::uint32_t i=0;i<output.vertices.size();++i)vertexByKey.try_emplace(Key(output.vertices[i].position),i);
            const auto addVertex=[&](Vertex vertex)
            {
                const GeomKey key=Key(vertex.position);const auto found=vertexByKey.find(key);
                if(found!=vertexByKey.end())return found->second;
                const auto index=static_cast<std::uint32_t>(output.vertices.size());
                output.vertices.push_back(vertex);vertexByKey.emplace(key,index);return index;
            };
            std::uint64_t splitTriangles=0;
            for(const Triangle& triangle:mesh.triangles)
            {
                std::vector<std::uint32_t> polygon;polygon.reserve(8);bool split=false;
                for(int edge=0;edge<3;++edge)
                {
                    const auto aIndex=triangle.index[edge],bIndex=triangle.index[(edge+1)%3];
                    const V3 pa=mesh.vertices[aIndex].position,pb=mesh.vertices[bIndex].position;
                    polygon.push_back(addVertex(mesh.vertices[aIndex]));
                    if(!boundaryEdges.contains(OrderedEdge(Key(pa),Key(pb))))continue;
                    const V3 direction=pb-pa;const double denominator=Dot(direction,direction);
                    if(denominator<=PositionQuantum*PositionQuantum)continue;
                    std::vector<std::pair<double,V3>> interior;
                    for(const V3 q:points)
                    {
                        if(q.x<std::min(pa.x,pb.x)-PositionQuantum||q.x>std::max(pa.x,pb.x)+PositionQuantum||
                            q.y<std::min(pa.y,pb.y)-PositionQuantum||q.y>std::max(pa.y,pb.y)+PositionQuantum||
                            q.z<std::min(pa.z,pb.z)-PositionQuantum||q.z>std::max(pa.z,pb.z)+PositionQuantum)continue;
                        const double t=Dot(q-pa,direction)/denominator;
                        if(t<=1.0e-8||t>=1.0-1.0e-8||Length(pa+direction*t-q)>PositionQuantum)continue;
                        interior.emplace_back(t,q);
                    }
                    std::sort(interior.begin(),interior.end(),[](const auto& x,const auto& y){return x.first<y.first;});
                    double previous=-1.0;
                    for(const auto& [t,p]:interior)
                    {
                        if(previous>=0.0&&std::abs(t-previous)<1.0e-8)continue;
                        Vertex vertex;vertex.position=p;vertex.normal=Normalize(mesh.vertices[aIndex].normal+
                            (mesh.vertices[bIndex].normal-mesh.vertices[aIndex].normal)*t);
                        polygon.push_back(addVertex(vertex));previous=t;split=true;
                    }
                }
                if(!split){output.triangles.push_back(triangle);continue;}
                ++splitTriangles;
                Vertex center;center.position=(mesh.vertices[triangle.index[0]].position+
                    mesh.vertices[triangle.index[1]].position+mesh.vertices[triangle.index[2]].position)/3.0;
                center.normal=Normalize(mesh.vertices[triangle.index[0]].normal+mesh.vertices[triangle.index[1]].normal+
                    mesh.vertices[triangle.index[2]].normal);
                const auto centerIndex=addVertex(center);
                for(std::size_t i=0;i<polygon.size();++i)
                {
                    const auto p0=polygon[i],p1=polygon[(i+1)%polygon.size()];
                    const V3 at=(output.vertices[centerIndex].position+output.vertices[p0].position+
                        output.vertices[p1].position)/3.0;
                    const FieldSample sample=field.Eval(at);
                    AddTriangle(output,centerIndex,p0,p1,triangle.feature,sample.gradient);
                }
            }
            if(splitTriangles==0)break;
            result.triangles+=splitTriangles;output.hash=HashMesh(output);mesh=std::move(output);
        }
        // Edge subdivision preserves each source triangle's planar footprint exactly.
        return result;
    }

    bool WritePhaseCsv(const std::filesystem::path& path,int phaseCount)
    {
        std::ofstream csv(path,std::ios::binary);
        csv<<"piece,method,phase_x_m,phase_y_m,phase_z_m,mean_abs_field_residual_m,max_abs_field_residual_m,probe_count,triangles,boundary_edges,nonmanifold_edges,crease_edge_pairs,normal_split_pairs,normal_split_rate,hash\n";
        const std::array<Piece,5> pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}};
        const std::array<Method,3> methods{{Method::MC,Method::Hybrid,Method::DC}};
        const std::array<double,5> phase{{0.0,0.05,0.10,0.15,0.20}};
        constexpr int n=32;
        std::uint64_t completed=0;
        for(Piece piece:pieces)
        {
            for(int px=0;px<phaseCount;++px)for(int py=0;py<phaseCount;++py)for(int pz=0;pz<phaseCount;++pz)
            {
                const V3 offset{phase[px],phase[py],phase[pz]};
                const Field field(piece,offset);
                Grid grid=SampleGrid(field,n,n,n,1,{-4.0,-4.0,-4.0},8);
                for(Method method:methods)
                {
                    Mesh mesh=BuildMesh(grid,field,method,8);
                    const ErrorMetrics error=SurfaceError(mesh,field);
                    const TopologyMetrics topology=AnalyzeTopology(mesh);
                    const double split=topology.creaseEdgePairs?static_cast<double>(topology.splitPairs)/topology.creaseEdgePairs:0.0;
                    csv<<PieceName(piece)<<','<<MethodName(method)<<','<<offset.x<<','<<offset.y<<','<<offset.z<<','
                        <<error.mean<<','<<error.maximum<<','<<error.samples<<','<<mesh.triangles.size()<<','
                        <<topology.boundaryEdges<<','<<topology.nonManifoldEdges<<','<<topology.creaseEdgePairs<<','
                        <<topology.splitPairs<<','<<split<<",0x"<<std::hex<<mesh.hash<<std::dec<<"\n";
                    ++completed;
                }
                if((px*phaseCount*phaseCount+py*phaseCount+pz)%25==0)
                    std::cout<<"[PhaseSweep] piece="<<PieceName(piece)<<" completed="<<completed<<"/"<<(5ull*phaseCount*phaseCount*phaseCount*3ull)<<"\n";
            }
        }
        return csv.good();
    }

    bool RunCostsAndDeterminism(const std::filesystem::path& outDir)
    {
        std::ofstream cost(outDir/"cost.csv",std::ios::binary);
        std::ofstream det(outDir/"determinism.csv",std::ios::binary);
        cost<<"piece,method,workers,cell_count,triangles,vertices,mesh_ms_median_3,mesh_hash,collision_proxy_hash,collision_geometry_delta_m\n";
        det<<"piece,method,workers,hash,identical_to_w1\n";
        const std::array<Piece,5> pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}};
        const std::array<Method,3> methods{{Method::MC,Method::Hybrid,Method::DC}};
        const std::array<int,3> workerCounts{{1,8,16}};
        const V3 phase{0.10,0.10,0.10};constexpr int n=32;
        for(Piece piece:pieces)
        {
            const Field field(piece,phase);
            Grid grid=SampleGrid(field,n,n,n,1,{-4.0,-4.0,-4.0},8);
            for(Method method:methods)
            {
                std::uint64_t reference=0;
                for(int workers:workerCounts)
                {
                    std::vector<double> times;Mesh last;
                    // W16 is a determinism stress count under an eight-CPU affinity mask. Use one
                    // sample there: repeated creation of 16 OS threads adds noise and had caused
                    // an access violation on this host without changing the worker-count gate.
                    const int repeats=workers==16?1:3;
                    for(int repeat=0;repeat<repeats;++repeat)
                    {
                        const auto start=std::chrono::steady_clock::now();
                        last=BuildMesh(grid,field,method,workers);
                        const auto finish=std::chrono::steady_clock::now();
                        times.push_back(std::chrono::duration<double,std::milli>(finish-start).count());
                    }
                    std::sort(times.begin(),times.end());
                    if(workers==1)reference=last.hash;
                    const bool same=last.hash==reference;
                    det<<PieceName(piece)<<','<<MethodName(method)<<','<<workers<<",0x"<<std::hex<<last.hash<<std::dec<<','<<(same?"true":"false")<<"\n";
                    const double reportedTime=times[times.size()/2];
                    cost<<PieceName(piece)<<','<<MethodName(method)<<','<<workers<<','<<grid.CellCount()<<','
                        <<last.triangles.size()<<','<<last.vertices.size()<<','<<reportedTime<<",0x"<<std::hex<<last.hash
                        <<",0x"<<last.hash<<std::dec<<",0\n";
                    std::cout<<"[Cost] piece="<<PieceName(piece)<<" method="<<MethodName(method)<<" workers="<<workers
                        <<" median_ms="<<std::fixed<<std::setprecision(3)<<times[1]<<" triangles="<<last.triangles.size()
                        <<" hash="<<std::hex<<last.hash<<std::dec<<" determinism="<<(same?"PASS":"FAIL")<<"\n";
                }
            }
        }
        return cost.good()&&det.good();
    }

    bool RunSeams(const std::filesystem::path& outDir,int rowLimit=0)
    {
        std::ofstream csv(outDir/"seams.csv",std::ios::binary);
        csv<<"piece,method,phase_x_m,phase_y_m,phase_z_m,tile_contour_expected,tile_unmatched_probes,tile_edge_mean_m,tile_edge_max_m,tile_boundary_edges,tile_nonmanifold_edges,tile_repair_triangles,tile_repair_max_residual_m,lod01_contour_expected,lod01_unmatched_probes,lod01_edge_mean_m,lod01_edge_max_m,lod01_boundary_edges,lod01_nonmanifold_edges,lod_internal_boundary_edges,lod_internal_nonmanifold_edges,lod_open_edges_before,lod_nonmanifold_before,lod_transition_triangles,lod_paired_loops,lod_capped_loops,lod_unmatched_loops,lod_transition_surface_max_m,lod_raw_contour_max_m,tile_left_segments,tile_right_segments,fine_segments,coarse_segments\n";
        const std::array<Piece,5> pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}};
        const std::array<Method,3> methods{{Method::MC,Method::Hybrid,Method::DC}};
        const std::array<double,5> phaseValues{{0.0,0.05,0.10,0.15,0.20}};
        std::uint64_t completed=0;
        for(Piece piece:pieces)
        {
            for(double px:phaseValues)for(double py:phaseValues)for(double pz:phaseValues)
            {
                const V3 phase{px,py,pz};
                const Field field(piece,phase);
                const bool seamExpected=HasIsoContourOnX(field,0.0,-4.0,4.0,VoxelMeters);
                const Grid tileLeft=SampleGrid(field,32,32,32,1,{-8.0,-4.0,-4.0},8,2);
                const Grid tileRight=SampleGrid(field,32,32,32,1,{0.0,-4.0,-4.0},8,2);
                const Grid lodFineLeft=SampleGrid(field,16,32,32,1,{-4.0,-4.0,-4.0},8,2);
                const Grid lodCoarseRight=SampleGrid(field,8,16,16,2,{0.0,-4.0,-4.0},8,2);
                for(Method method:methods)
                {
                    Mesh leftMesh=BuildMesh(tileLeft,field,method,8);
                    Mesh rightMesh=BuildMesh(tileRight,field,method,8);
                    // Build both independently sampled halo meshes without border-quad ownership,
                    // then clip each complete surface at x=0. This gives the transition an exact
                    // geometric cut line instead of trying to pair unrelated interior holes.
                    const Mesh lodFineMesh=BuildMesh(lodFineLeft,field,method,8,false);
                    const Mesh lodCoarseMesh=BuildMesh(lodCoarseRight,field,method,8,false);
                    // Ownership suppresses the same border quad in one rendered tile. For the
                    // independent-data witness, re-emit halo quads in each temporary mesh and
                    // compare the two independently solved contours; these witness meshes are
                    // excluded from the merged topology audit.
                    Mesh leftWitness,rightWitness;
                    const Mesh* leftContourMesh=&leftMesh;
                    const Mesh* rightContourMesh=&rightMesh;
                    if(method!=Method::MC)
                    {
                        leftWitness=BuildMesh(tileLeft,field,method,8,false);
                        rightWitness=BuildMesh(tileRight,field,method,8,false);
                        leftContourMesh=&leftWitness;rightContourMesh=&rightWitness;
                    }
                    const auto leftSegments=SeamSegments(*leftContourMesh,0.0,-1);
                    const auto rightSegments=SeamSegments(*rightContourMesh,0.0,1);
                    SeamMetric tile=CompareSegments(leftSegments,rightSegments);
                    Mesh lodFineCut=ClipMeshAtX(lodFineMesh,field,true);
                    Mesh lodCoarseCut=ClipMeshAtX(lodCoarseMesh,field,false);
                    ConformOpenTJunctions(lodFineCut,field);
                    ConformOpenTJunctions(lodCoarseCut,field);
                    const auto lodFineSegments=BoundarySegmentsAtX(lodFineCut,0.0);
                    const auto lodCoarseSegments=BoundarySegmentsAtX(lodCoarseCut,0.0);
                    const SeamMetric rawLod=CompareSegments(lodFineSegments,lodCoarseSegments);
                    SeamMetric lod=rawLod;
                    if(seamExpected&&leftSegments.empty()&&rightSegments.empty())
                    {tile.unmatched=1;tile.mean=-1;tile.maximum=-1;}
                    if(seamExpected&&lodFineSegments.empty()&&lodCoarseSegments.empty())
                    {lod.unmatched=1;lod.mean=-1;lod.maximum=-1;}
                    Mesh tileMerged=MergeMeshes(leftMesh,rightMesh);
                    BoundaryRepairAudit tileRepair=ConformOpenTJunctions(tileMerged,field);
                    const TopologyMetrics tileTopology=AnalyzeTopology(tileMerged);
                    tile.unmatched+=tileRepair.unmatchedEdges+tileTopology.boundaryEdges+tileTopology.nonManifoldEdges;
                    Mesh lodMerged=MergeMeshes(lodFineCut,lodCoarseCut);
                    const TopologyMetrics lodBefore=AnalyzeTopology(lodMerged);
                    const BoundaryLoopSet fineLoops=ExtractBoundaryLoops(lodFineCut);
                    const BoundaryLoopSet coarseLoops=ExtractBoundaryLoops(lodCoarseCut);
                    const std::size_t firstTransitionTriangle=lodMerged.triangles.size();
                    const TransitionAudit transition=StitchBoundaryLoops(lodMerged,fineLoops,coarseLoops,field);
                    const SeamCapAudit seamCaps=CapClosedBoundaryLoopsInXBand(lodMerged,field,0.0,0.003);
                    const TopologyMetrics lodAfter=AnalyzeTopology(lodMerged);
                    const TopologyMetrics lodSeamAfter=AnalyzeTopologyInXBand(lodMerged,0.0,0.003);
                    Mesh transitionOnly;transitionOnly.vertices=lodMerged.vertices;
                    transitionOnly.triangles.assign(lodMerged.triangles.begin()+static_cast<std::ptrdiff_t>(firstTransitionTriangle),lodMerged.triangles.end());
                    const ErrorMetrics transitionError=SurfaceError(transitionOnly,field);
                    lod.unmatched=rawLod.unmatched+transition.unmatchedLoops+seamCaps.unmatchedLoops+
                        lodSeamAfter.boundaryEdges+lodSeamAfter.nonManifoldEdges;
                    if(seamExpected&&rawLod.probes==0)++lod.unmatched;
                    if(seamExpected&&transition.pairedLoops==0&&transition.cappedLoops==0&&seamCaps.cappedLoops==0)++lod.unmatched;
                    if(rawLod.maximum>LodContourErrorBudgetMeters||transitionError.maximum>LodContourErrorBudgetMeters)++lod.unmatched;
                    csv<<PieceName(piece)<<','<<MethodName(method)<<','<<phase.x<<','<<phase.y<<','<<phase.z<<','
                        <<(seamExpected?"true":"false")<<','<<tile.unmatched<<','<<tile.mean<<','<<tile.maximum<<','
                        <<tileTopology.boundaryEdges<<','<<tileTopology.nonManifoldEdges<<','
                        <<tileRepair.triangles<<','<<tileRepair.maximumResidual<<','
                        <<(seamExpected?"true":"false")<<','<<lod.unmatched<<','<<lod.mean<<','<<lod.maximum<<','
                        <<lodSeamAfter.boundaryEdges<<','<<lodSeamAfter.nonManifoldEdges<<','
                        <<(lodAfter.boundaryEdges-lodSeamAfter.boundaryEdges)<<','
                        <<(lodAfter.nonManifoldEdges-lodSeamAfter.nonManifoldEdges)<<','
                        <<lodBefore.boundaryEdges<<','<<lodBefore.nonManifoldEdges<<','
                        <<(transition.transitionTriangles+seamCaps.triangles)<<','<<transition.pairedLoops<<','<<(transition.cappedLoops+seamCaps.cappedLoops)<<','<<(transition.unmatchedLoops+seamCaps.unmatchedLoops)<<','
                        <<transitionError.maximum<<','<<rawLod.maximum<<','
                        <<leftSegments.size()<<','<<rightSegments.size()<<','<<lodFineSegments.size()<<','<<lodCoarseSegments.size()<<"\n";
                    if(px==0.10&&py==0.10&&pz==0.10)
                    {
                        const auto seamObj=outDir/"meshes"/(std::string("lodseam_")+PieceName(piece)+"_"+MethodName(method)+".obj");
                        WriteObj(seamObj,lodMerged);
                    }
                    ++completed;
                    if(rowLimit>0&&completed>=static_cast<std::uint64_t>(rowLimit))return csv.good();
                    if(completed%75==0)std::cout<<"[SeamSweep] completed="<<completed<<"/1875\n";
                }
            }
        }
        return csv.good();
    }

    bool RunSharpCliffBorderCase(const std::filesystem::path& outDir)
    {
        const Field field(Piece::Ledge,{0.25,0.10,0.10});
        const bool expected=HasIsoContourOnX(field,0.0,-4.0,4.0,VoxelMeters);
        const Grid fineGrid=SampleGrid(field,16,32,32,1,{-4.0,-4.0,-4.0},8,2);
        const Grid coarseGrid=SampleGrid(field,8,16,16,2,{0.0,-4.0,-4.0},8,2);
        std::ofstream csv(outDir/"lod_cliff_crossing.csv",std::ios::binary);
        csv<<"method,phase_x_m,phase_y_m,phase_z_m,contour_expected,unmatched_probes,contour_mean_m,contour_max_m,boundary_edges,nonmanifold_edges,edges_before,nonmanifold_before,transition_triangles,paired_loops,capped_loops,unmatched_loops,transition_surface_max_m,fine_segments,coarse_segments\n";
        for(Method method:{Method::MC,Method::Hybrid,Method::DC})
        {
            const Mesh fine=BuildMesh(fineGrid,field,method,8,false);
            const Mesh coarse=BuildMesh(coarseGrid,field,method,8,false);
            Mesh fineCut=ClipMeshAtX(fine,field,true),coarseCut=ClipMeshAtX(coarse,field,false);
            const auto fineSegments=BoundarySegmentsAtX(fineCut,0.0);
            const auto coarseSegments=BoundarySegmentsAtX(coarseCut,0.0);
            SeamMetric contour=CompareSegments(fineSegments,coarseSegments);
            if(expected&&fineSegments.empty()&&coarseSegments.empty())contour.unmatched=1;
            Mesh merged=MergeMeshes(fineCut,coarseCut);
            const TopologyMetrics before=AnalyzeTopology(merged);
            const BoundaryLoopSet fineLoops=ExtractBoundaryLoops(fineCut),coarseLoops=ExtractBoundaryLoops(coarseCut);
            const std::size_t transitionStart=merged.triangles.size();
            const TransitionAudit transition=StitchBoundaryLoops(merged,fineLoops,coarseLoops,field);
            const SeamCapAudit seamCaps=CapClosedBoundaryLoopsInXBand(merged,field,0.0,0.003);
            const TopologyMetrics after=AnalyzeTopology(merged);
            Mesh transitionOnly;transitionOnly.vertices=merged.vertices;
            transitionOnly.triangles.assign(merged.triangles.begin()+static_cast<std::ptrdiff_t>(transitionStart),merged.triangles.end());
            const ErrorMetrics residual=SurfaceError(transitionOnly,field);
            const std::uint64_t unmatched=contour.unmatched+transition.unmatchedLoops+seamCaps.unmatchedLoops+
                after.boundaryEdges+after.nonManifoldEdges+
                ((expected&&transition.pairedLoops==0&&transition.cappedLoops==0&&seamCaps.cappedLoops==0)?1:0)+
                ((contour.maximum>LodContourErrorBudgetMeters||residual.maximum>LodContourErrorBudgetMeters)?1:0);
            csv<<MethodName(method)<<",0.25,0.1,0.1,"<<(expected?"true":"false")<<','<<unmatched<<','
                <<contour.mean<<','<<contour.maximum<<','<<after.boundaryEdges<<','<<after.nonManifoldEdges<<','
                <<before.boundaryEdges<<','<<before.nonManifoldEdges<<','<<(transition.transitionTriangles+seamCaps.triangles)<<','
                <<transition.pairedLoops<<','<<(transition.cappedLoops+seamCaps.cappedLoops)<<','<<(transition.unmatchedLoops+seamCaps.unmatchedLoops)<<','<<residual.maximum<<','
                <<fineSegments.size()<<','<<coarseSegments.size()<<"\n";
            WriteObj(outDir/"meshes"/(std::string("lodseam_cliff_crossing_")+MethodName(method)+".obj"),merged);
            std::cout<<"[LodCliffCrossing] method="<<MethodName(method)<<" expected="<<(expected?"true":"false")
                <<" contour_max_m="<<contour.maximum<<" unmatched="<<unmatched<<" open="<<after.boundaryEdges
                <<" nonmanifold="<<after.nonManifoldEdges<<"\n";
        }
        return csv.good();
    }

    bool MeshesBitIdentical(const Mesh& a,const Mesh& b)
    {
        if(a.vertices.size()!=b.vertices.size()||a.triangles.size()!=b.triangles.size())return false;
        for(std::size_t i=0;i<a.vertices.size();++i)
        {
            const Vertex& x=a.vertices[i];const Vertex& y=b.vertices[i];
            const std::array<double,6> xa{{x.position.x,x.position.y,x.position.z,x.normal.x,x.normal.y,x.normal.z}};
            const std::array<double,6> ya{{y.position.x,y.position.y,y.position.z,y.normal.x,y.normal.y,y.normal.z}};
            for(std::size_t c=0;c<xa.size();++c)
                if(std::bit_cast<std::uint64_t>(xa[c])!=std::bit_cast<std::uint64_t>(ya[c]))return false;
        }
        for(std::size_t i=0;i<a.triangles.size();++i)
            if(a.triangles[i].index!=b.triangles[i].index||a.triangles[i].feature!=b.triangles[i].feature)return false;
        return true;
    }

    bool RunNativeSurfaceComparison(const std::filesystem::path& outDir)
    {
        const auto manifestPath=outDir.parent_path().parent_path()/"WorldAssemblerM1"/"round4"/"world_assembly.json";
        std::ifstream input(manifestPath,std::ios::binary);
        if(!input)return false;
        const std::string json((std::istreambuf_iterator<char>(input)),std::istreambuf_iterator<char>());
        std::string error;
        auto manifest=VoxelForgeCore::WorldAssemblyManifest::FromJson(json,error);
        if(!manifest){std::cerr<<"[NativeSurface] manifest rejected: "<<error<<"\n";return false;}
        const auto stratum=std::find_if(manifest->plan.strates.begin(),manifest->plan.strates.end(),
            [](const auto& item){return item.strateIndex==1;});
        if(stratum==manifest->plan.strates.end())return false;
        const auto feature=std::find_if(stratum->features.begin(),stratum->features.end(),
            [](const auto& item){return item.type==4;});
        if(feature==stratum->features.end())return false;
        VoxelForgeCore::Heightfield2D heightfield(manifest->recipe,*stratum);
        if(!heightfield.isValid())return false;
        constexpr double SubvoxelsPerMeter=1024.0;
        const double centerX=(feature->minimumSubvoxels[0]+feature->maximumSubvoxels[0])/(2.0*SubvoxelsPerMeter);
        const double centerY=(feature->minimumSubvoxels[1]+feature->maximumSubvoxels[1])/(2.0*SubvoxelsPerMeter);
        const double minFeatureZ=feature->minimumSubvoxels[2]/SubvoxelsPerMeter;
        const double maxFeatureZ=feature->maximumSubvoxels[2]/SubvoxelsPerMeter;
        const Field field(heightfield);
        double previousZ=minFeatureZ-8.0,previousDensity=field.Eval({centerX,centerY,previousZ}).density;
        double surfaceZ=(minFeatureZ+maxFeatureZ)*0.5;bool foundCrossing=false;
        for(double z=previousZ+0.25;z<=maxFeatureZ+8.0;z+=0.25)
        {
            const double density=field.Eval({centerX,centerY,z}).density;
            if((previousDensity<0.0)!=(density<0.0))
            {surfaceZ=(previousZ+z)*0.5;foundCrossing=true;break;}
            previousZ=z;previousDensity=density;
        }
        if(!foundCrossing)
        {std::cerr<<"[NativeSurface] no iso crossing at feature center x="<<centerX<<" y="<<centerY<<"\n";return false;}
        constexpr int cells=32;constexpr int step=1; // 25 cm cells over an 8 m native cliff patch.
        const V3 origin{centerX-4.0,centerY-4.0,surfaceZ-4.0};
        const Grid grid=SampleGrid(field,cells,cells,cells,step,origin,8);
        std::ofstream csv(outDir/"native_surface.csv",std::ios::binary);
        csv<<"method,manifest_hash,feature_id,feature_type,grid_cells,cell_size_m,origin_x_m,origin_y_m,origin_z_m,iso_probe_z_m,exactly_same_as_mc,mean_abs_field_residual_m,max_abs_field_residual_m,triangles,vertices,boundary_edges,nonmanifold_edges,hash\n";
        Mesh mc;
        for(Method method:{Method::MC,Method::Hybrid,Method::DC})
        {
            Mesh mesh=BuildMesh(grid,field,method,8);
            if(method==Method::MC)mc=mesh;
            const bool identical=method==Method::MC||MeshesBitIdentical(mc,mesh);
            const ErrorMetrics residual=SurfaceError(mesh,field);
            const TopologyMetrics topology=AnalyzeTopology(mesh);
            WriteObj(outDir/"meshes"/(std::string(PieceName(Piece::NativeSurface))+"_"+MethodName(method)+".obj"),mesh);
            csv<<MethodName(method)<<",0x"<<std::hex<<manifest->contentHash<<std::dec<<",0x"<<std::hex<<feature->stableId<<std::dec
                <<','<<feature->type<<','<<grid.CellCount()<<','<<step*VoxelMeters<<','<<origin.x<<','<<origin.y<<','<<origin.z<<','<<surfaceZ
                <<','<<(identical?"true":"false")<<','<<residual.mean<<','<<residual.maximum<<','<<mesh.triangles.size()<<','
                <<mesh.vertices.size()<<','<<topology.boundaryEdges<<','<<topology.nonManifoldEdges<<",0x"<<std::hex<<mesh.hash<<std::dec<<"\n";
            std::cout<<"[NativeSurface] method="<<MethodName(method)<<" center=("<<centerX<<','<<centerY<<','<<surfaceZ
                <<") triangles="<<mesh.triangles.size()<<" hash=0x"<<std::hex<<mesh.hash<<std::dec
                <<" identical_to_mc="<<(identical?"true":"false")<<"\n";
        }
        return csv.good();
    }

    bool WriteCaptureMeshes(const std::filesystem::path& outDir)
    {
        const std::array<Piece,5> pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}};
        const std::array<Method,3> methods{{Method::MC,Method::Hybrid,Method::DC}};
        const V3 phase{0.10,0.10,0.10};constexpr int n=32;
        for(Piece piece:pieces)
        {
            const Field field(piece,phase);Grid grid=SampleGrid(field,n,n,n,1,{-4.0,-4.0,-4.0},8);
            for(Method method:methods)
            {
                const Mesh mesh=BuildMesh(grid,field,method,8);
                WriteObj(outDir/"meshes"/(std::string(PieceName(piece))+"_"+MethodName(method)+".obj"),mesh);
            }
        }
        return true;
    }
}

int main(int argc,char** argv)
{
    try
    {
        std::filesystem::path outDir="Standalone/Research/WorldAssemblerM2/results";
        int phaseCount=5;
        bool seamsOnly=false;
        bool nativeOnly=false;
        bool analysisOnly=false;
        int seamRowLimit=0;
        for(int i=1;i<argc;++i)
        {
            const std::string arg=argv[i];
            if(arg.starts_with("--out="))outDir=arg.substr(6);
            else if(arg.starts_with("--phase-count="))phaseCount=std::clamp(std::stoi(arg.substr(14)),2,5);
            else if(arg=="--seams-only")seamsOnly=true;
            else if(arg=="--native-only")nativeOnly=true;
            else if(arg=="--analysis-only")analysisOnly=true;
            else if(arg.starts_with("--seam-row-limit="))seamRowLimit=std::max(0,std::stoi(arg.substr(17)));
            else {std::cerr<<"Unknown argument: "<<arg<<"\n";return 2;}
        }
        std::filesystem::create_directories(outDir/"meshes");
        const auto start=std::chrono::steady_clock::now();
        std::cout<<"[SharpMesherSpike] version=1 voxel_cm=25 quantization=1/65536_voxel phase_grid="<<phaseCount
            <<"^3 processes=BelowNormal affinity=0xFF00\n";
        if(!seamsOnly&&!nativeOnly&&!WritePhaseCsv(outDir/"phase_sweep.csv",phaseCount))throw std::runtime_error("phase CSV write failed");
        if(!seamsOnly&&!nativeOnly&&!RunCostsAndDeterminism(outDir))throw std::runtime_error("cost/determinism audit failed");
        if(!nativeOnly&&!analysisOnly&&!RunSeams(outDir,seamRowLimit))throw std::runtime_error("seam audit failed");
        if(!nativeOnly&&!analysisOnly&&!RunSharpCliffBorderCase(outDir))throw std::runtime_error("sharp cliff LOD border case failed");
        if(!seamsOnly&&!nativeOnly&&!WriteCaptureMeshes(outDir))throw std::runtime_error("capture mesh export failed");
        if(!seamsOnly&&!RunNativeSurfaceComparison(outDir))throw std::runtime_error("native SurfaceWorld comparison failed");
        const double seconds=std::chrono::duration<double>(std::chrono::steady_clock::now()-start).count();
        std::cout<<"[SharpMesherSpike] complete seconds="<<seconds<<" phase_count="<<phaseCount
            <<" phase_runs="<<(5*phaseCount*phaseCount*phaseCount*3)<<"\n";
        return 0;
    }
    catch(const std::exception& e)
    {
        std::cerr<<"[SharpMesherSpike] failed: "<<e.what()<<"\n";return 1;
    }
}