VoxelForge/Standalone/Source/LinuxScreenSpaceGate.h

#pragma once

#include "StreamingWorld.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <limits>
#include <optional>
#include <span>
#include <utility>
#include <vector>

namespace StandaloneScreenSpaceGate
{
    inline constexpr double OnePixelRadians=1.57079632679489661923/2560.0;

    struct FErrorRow
    {
        int Lod=0;
        std::uint64_t Samples=0,Missing=0,Floating=0;
        double P50=0.0,P99=0.0,Maximum=0.0;
        double VerticalP50=0.0,VerticalP99=0.0,VerticalMaximum=0.0;
        double MaximumNormalDistanceMeters=0.0,MaximumX=0.0,MaximumY=0.0;
        double MaximumReferenceZ=0.0,MaximumDisplayedZ=0.0;
        double MaximumSlopeDegrees=0.0;
        std::array<std::uint64_t,4> SlopeSamples{};
        std::array<double,4> SlopeP99{},SlopeMaximum{};
        bool Valid=false;
    };

    struct FPopRow
    {
        int CoarseLod=0;
        std::uint64_t Samples=0,Cracks=0,Floating=0;
        double P50=0.0,P99=0.0,Maximum=0.0;
        double VerticalP50=0.0,VerticalP99=0.0,VerticalMaximum=0.0;
        double MaximumNormalDistanceMeters=0.0,MaximumX=0.0,MaximumY=0.0;
        double MaximumReferenceZ=0.0,MaximumDisplayedZ=0.0;
        double MaximumSlopeDegrees=0.0;
        std::array<std::uint64_t,4> SlopeSamples{};
        std::array<double,4> SlopeP99{},SlopeMaximum{};
        std::vector<double> SamplePixels;
        bool Valid=false;
    };

    struct FSurfaceRoot
    {
        double ZMeters=0.0;
        double SlopeX=0.0;
        double SlopeY=0.0;
    };

    struct FMeshSurfaceHit
    {
        double Z=0.0;
        std::array<double,3> Normal{};
        std::array<std::uint8_t,4> Color{};
        bool Ceiling=false;
    };

    struct FViewBasis
    {
        std::array<double,3> Right{};
        std::array<double,3> Up{};
    };

    struct FErrorAccumulator
    {
        std::vector<double> Pixels,VerticalPixels;
        std::array<std::vector<double>,4> BySlope;
        double Maximum=-1.0,MaximumVertical=0.0,MaximumNormalDistanceMeters=0.0;
        double MaximumX=0.0,MaximumY=0.0,MaximumSlopeDegrees=0.0;
        double MaximumReferenceZ=0.0,MaximumDisplayedZ=0.0;
        std::array<std::uint64_t,4> SlopeSamples{};
    };

    inline double Percentile(std::vector<double>& Values,double P);

    inline FViewBasis MakeViewBasis(double YawDegrees,double PitchDegrees)
    {
        constexpr double Pi=3.14159265358979323846;
        const double Yaw=YawDegrees*Pi/180.0,Pitch=PitchDegrees*Pi/180.0;
        const std::array<double,3> Forward={std::sin(Yaw)*std::cos(Pitch),
            std::cos(Yaw)*std::cos(Pitch),std::sin(Pitch)};
        FViewBasis Basis;
        Basis.Right={std::cos(Yaw),-std::sin(Yaw),0.0};
        Basis.Up={Basis.Right[1]*Forward[2]-Basis.Right[2]*Forward[1],
            Basis.Right[2]*Forward[0]-Basis.Right[0]*Forward[2],
            Basis.Right[0]*Forward[1]-Basis.Right[1]*Forward[0]};
        return Basis;
    }

    inline std::size_t SlopeBin(double Degrees)
    {
        if(Degrees<30.0)return 0;
        if(Degrees<45.0)return 1;
        if(Degrees<60.0)return 2;
        return 3;
    }

    inline void AddError(FErrorAccumulator& Accumulator,double DeltaZMeters,
        const FSurfaceRoot& Root,double PixelBudgetMeters,const FViewBasis& View,
        double X,double Y)
    {
        constexpr double Pi=3.14159265358979323846;
        const double GradientSquared=Root.SlopeX*Root.SlopeX+Root.SlopeY*Root.SlopeY;
        const double NormalLength=std::sqrt(1.0+GradientSquared);
        const double NormalX=-Root.SlopeX/NormalLength;
        const double NormalY=-Root.SlopeY/NormalLength;
        const double NormalZ=1.0/NormalLength;
        const double NormalDistance=std::abs(DeltaZMeters)/NormalLength;
        const double ProjectedNormal=std::hypot(
            NormalX*View.Right[0]+NormalY*View.Right[1]+NormalZ*View.Right[2],
            NormalX*View.Up[0]+NormalY*View.Up[1]+NormalZ*View.Up[2]);
        const double ProjectedVertical=std::hypot(View.Right[2],View.Up[2]);
        const double NormalPixels=NormalDistance*ProjectedNormal/PixelBudgetMeters;
        const double VerticalPixels=std::abs(DeltaZMeters)*ProjectedVertical/PixelBudgetMeters;
        const double SlopeDegrees=std::atan(std::sqrt(GradientSquared))*180.0/Pi;
        Accumulator.Pixels.push_back(NormalPixels);
        Accumulator.VerticalPixels.push_back(VerticalPixels);
        const std::size_t Bin=SlopeBin(SlopeDegrees);
        Accumulator.BySlope[Bin].push_back(NormalPixels);
        ++Accumulator.SlopeSamples[Bin];
        if(NormalPixels>Accumulator.Maximum)
        {
            Accumulator.Maximum=NormalPixels;
            Accumulator.MaximumVertical=VerticalPixels;
            Accumulator.MaximumNormalDistanceMeters=NormalDistance;
            Accumulator.MaximumX=X;Accumulator.MaximumY=Y;
            Accumulator.MaximumSlopeDegrees=SlopeDegrees;
            Accumulator.MaximumReferenceZ=Root.ZMeters;
            Accumulator.MaximumDisplayedZ=Root.ZMeters+DeltaZMeters;
        }
    }

    inline void AddMissing(FErrorAccumulator& Accumulator)
    {
        Accumulator.Pixels.push_back(2.0);
        Accumulator.VerticalPixels.push_back(2.0);
    }

    inline void FinalizeError(FErrorRow& Row,FErrorAccumulator& Accumulator)
    {
        Row.Samples=Accumulator.Pixels.size();
        Row.P50=Percentile(Accumulator.Pixels,0.50);
        Row.P99=Percentile(Accumulator.Pixels,0.99);
        Row.Maximum=Accumulator.Pixels.empty()?0.0
            :*std::max_element(Accumulator.Pixels.begin(),Accumulator.Pixels.end());
        Row.VerticalP50=Percentile(Accumulator.VerticalPixels,0.50);
        Row.VerticalP99=Percentile(Accumulator.VerticalPixels,0.99);
        Row.VerticalMaximum=Accumulator.VerticalPixels.empty()?0.0
            :*std::max_element(Accumulator.VerticalPixels.begin(),Accumulator.VerticalPixels.end());
        Row.MaximumNormalDistanceMeters=Accumulator.MaximumNormalDistanceMeters;
        Row.MaximumX=Accumulator.MaximumX;Row.MaximumY=Accumulator.MaximumY;
        Row.MaximumReferenceZ=Accumulator.MaximumReferenceZ;
        Row.MaximumDisplayedZ=Accumulator.MaximumDisplayedZ;
        Row.MaximumSlopeDegrees=Accumulator.MaximumSlopeDegrees;
        Row.SlopeSamples=Accumulator.SlopeSamples;
        for(std::size_t Bin=0;Bin<Accumulator.BySlope.size();++Bin)
        {
            Row.SlopeP99[Bin]=Percentile(Accumulator.BySlope[Bin],0.99);
            Row.SlopeMaximum[Bin]=Accumulator.BySlope[Bin].empty()?0.0
                :*std::max_element(Accumulator.BySlope[Bin].begin(),Accumulator.BySlope[Bin].end());
        }
    }

    struct FHeightGrid
    {
        int Side=65;
        double MinimumX=0.0,MaximumX=0.0,MinimumY=0.0,MaximumY=0.0;
        std::vector<std::vector<FMeshSurfaceHit>> Hits;

        FHeightGrid()=default;
        FHeightGrid(int InSide,double MinX,double MaxX,double MinY,double MaxY)
            :Side(InSide),MinimumX(MinX),MaximumX(MaxX),MinimumY(MinY),MaximumY(MaxY),
             Hits(static_cast<std::size_t>(InSide)*InSide)
        {}

        double X(int Ix) const
        {return MinimumX+(MaximumX-MinimumX)*Ix/(Side-1);}
        double Y(int Iy) const
        {return MinimumY+(MaximumY-MinimumY)*Iy/(Side-1);}

        void Add(int Ix,int Iy,const FMeshSurfaceHit& Hit)
        {
            auto& Values=Hits[static_cast<std::size_t>(Iy)*Side+Ix];
            Values.push_back(Hit);
        }

        const std::vector<FMeshSurfaceHit>& At(int Ix,int Iy) const
        {return Hits[static_cast<std::size_t>(Iy)*Side+Ix];}
    };

    inline double Percentile(std::vector<double>& Values,double P)
    {
        if(Values.empty())return 0.0;
        std::sort(Values.begin(),Values.end());
        const std::size_t Index=std::min(Values.size()-1,
            static_cast<std::size_t>(std::ceil(P*Values.size())-1));
        return Values[Index];
    }

    inline double TileAabbDistanceMeters(const FVoxelTileKey& Key,
        const FVector& CameraVoxels,double VoxelMeters)
    {
        const double Extent=static_cast<double>(Key.ExtentVoxels());
        const double Min[3]={static_cast<double>(Key.OriginVoxels().X),
            static_cast<double>(Key.OriginVoxels().Y),static_cast<double>(Key.OriginVoxels().Z)};
        const double Cam[3]={CameraVoxels.X,CameraVoxels.Y,CameraVoxels.Z};
        double Squared=0.0;
        for(int Axis=0;Axis<3;++Axis)
        {
            const double Delta=Cam[Axis]<Min[Axis]?Min[Axis]-Cam[Axis]
                :Cam[Axis]>Min[Axis]+Extent?Cam[Axis]-(Min[Axis]+Extent):0.0;
            Squared+=Delta*Delta;
        }
        return std::sqrt(Squared)*VoxelMeters;
    }

    inline bool TileHasNativeSurface(const FVoxelTileKey& Key,
        const VoxelForgeCore::WorldAssemblyManifest& Assembly,
        double VoxelMeters)
    {
        const double Extent=static_cast<double>(Key.ExtentVoxels());
        const double X=(static_cast<double>(Key.OriginVoxels().X)+Extent*0.5)*256.0;
        const double Y=(static_cast<double>(Key.OriginVoxels().Y)+Extent*0.5)*256.0;
        const double MinZ=static_cast<double>(Key.OriginVoxels().Z);
        const double MaxZ=MinZ+Extent;
        for(const auto& Plan:Assembly.plan.strates)
        {
            VoxelForgeCore::Heightfield2D Field(Assembly.recipe,Plan);
            if(!Field.isValid())continue;
            const auto Column=Field.point(static_cast<std::int64_t>(std::llround(X)),
                static_cast<std::int64_t>(std::llround(Y)),0);
            const double Ground=static_cast<double>(Column.groundHeightSubvoxels)/256.0;
            const double Cap=static_cast<double>(Column.capHeightSubvoxels)/256.0;
            const double Margin=std::max(2.0,0.25/VoxelMeters);
            if(Ground<Cap&&((Ground>=MinZ-Margin&&Ground<=MaxZ+Margin)
                ||(Cap>=MinZ-Margin&&Cap<=MaxZ+Margin)))return true;
        }
        return false;
    }

    inline FHeightGrid Rasterize(std::span<const FStandaloneCpuTile* const> Tiles,
        const FVoxelTileKey& Patch,double VoxelMeters,int Side=65)
    {
        const double Extent=static_cast<double>(Patch.ExtentVoxels())*VoxelMeters;
        const double MinX=static_cast<double>(Patch.OriginVoxels().X)*VoxelMeters;
        const double MinY=static_cast<double>(Patch.OriginVoxels().Y)*VoxelMeters;
        FHeightGrid Grid(Side,MinX,MinX+Extent,MinY,MinY+Extent);
        for(const FStandaloneCpuTile* Tile:Tiles)
        {
            if(!Tile||Tile->Indices.size()<3)continue;
            const std::size_t TriangleCount=Tile->Indices.size()/3;
            const std::size_t CeilingStart=TriangleCount
                -std::min<std::size_t>(TriangleCount,Tile->CeilingTriangleCount);
            for(std::size_t Triangle=0;Triangle<TriangleCount;++Triangle)
            {
                if(Triangle<Tile->TriangleSkirtFaceMasks.size()
                    &&Tile->TriangleSkirtFaceMasks[Triangle]!=0)continue;
                const std::uint32_t I0=Tile->Indices[Triangle*3];
                const std::uint32_t I1=Tile->Indices[Triangle*3+1];
                const std::uint32_t I2=Tile->Indices[Triangle*3+2];
                if(I0>=Tile->Vertices.size()||I1>=Tile->Vertices.size()||I2>=Tile->Vertices.size())continue;
                const auto& A=Tile->Vertices[I0].Position;
                const auto& B=Tile->Vertices[I1].Position;
                const auto& C=Tile->Vertices[I2].Position;
                const double Denominator=(B[1]-C[1])*(A[0]-C[0])+(C[0]-B[0])*(A[1]-C[1]);
                if(std::abs(Denominator)<1e-10)continue;
                const double MinTriX=std::min({A[0],B[0],C[0]}),MaxTriX=std::max({A[0],B[0],C[0]});
                const double MinTriY=std::min({A[1],B[1],C[1]}),MaxTriY=std::max({A[1],B[1],C[1]});
                const double StepX=Extent/(Side-1),StepY=Extent/(Side-1);
                const int Ix0=std::clamp(static_cast<int>(std::ceil((MinTriX-MinX)/StepX)),0,Side-1);
                const int Ix1=std::clamp(static_cast<int>(std::floor((MaxTriX-MinX)/StepX)),0,Side-1);
                const int Iy0=std::clamp(static_cast<int>(std::ceil((MinTriY-MinY)/StepY)),0,Side-1);
                const int Iy1=std::clamp(static_cast<int>(std::floor((MaxTriY-MinY)/StepY)),0,Side-1);
                for(int Iy=Iy0;Iy<=Iy1;++Iy)
                for(int Ix=Ix0;Ix<=Ix1;++Ix)
                {
                    const double X=Grid.X(Ix),Y=Grid.Y(Iy);
                    const double W0=((B[1]-C[1])*(X-C[0])+(C[0]-B[0])*(Y-C[1]))/Denominator;
                    const double W1=((C[1]-A[1])*(X-C[0])+(A[0]-C[0])*(Y-C[1]))/Denominator;
                    const double W2=1.0-W0-W1;
                    if(W0<-1e-6||W1<-1e-6||W2<-1e-6)continue;
                    const auto& VA=Tile->Vertices[I0];
                    const auto& VB=Tile->Vertices[I1];
                    const auto& VC=Tile->Vertices[I2];
                    std::array<double,3> Normal{};
                    for(int Axis=0;Axis<3;++Axis)
                        Normal[Axis]=W0*VA.Normal[Axis]+W1*VB.Normal[Axis]+W2*VC.Normal[Axis];
                    const double NormalLength=std::sqrt(Normal[0]*Normal[0]
                        +Normal[1]*Normal[1]+Normal[2]*Normal[2]);
                    if(NormalLength>1.0e-12)
                        for(double& Component:Normal)Component/=NormalLength;
                    std::array<std::uint8_t,4> Color{};
                    if(Tile->VertexColors.size()==Tile->Vertices.size())
                        for(int Channel=0;Channel<4;++Channel)
                            Color[Channel]=static_cast<std::uint8_t>(std::clamp(std::lround(
                                W0*Tile->VertexColors[I0][Channel]
                                +W1*Tile->VertexColors[I1][Channel]
                                +W2*Tile->VertexColors[I2][Channel]),0l,255l));
                    Grid.Add(Ix,Iy,{W0*A[2]+W1*B[2]+W2*C[2],Normal,Color,Triangle>=CeilingStart});
                }
            }
        }
        for(auto& Values:Grid.Hits)
        {
            std::stable_sort(Values.begin(),Values.end(),[](const FMeshSurfaceHit& A,
                const FMeshSurfaceHit& B){return A.Z<B.Z;});
            Values.erase(std::unique(Values.begin(),Values.end(),[](const FMeshSurfaceHit& A,
                const FMeshSurfaceHit& B){return std::abs(A.Z-B.Z)<0.02;}),Values.end());
        }
        return Grid;
    }

    inline std::vector<FSurfaceRoot> ExactHeightRoots(const FHeightGrid& Grid,int Ix,int Iy,
        const FVoxelTileKey& Patch,double VoxelMeters,
        const VoxelForgeCore::WorldAssemblyManifest& Assembly)
    {
        std::vector<FSurfaceRoot> Roots;
        const double XSubvoxels=Grid.X(Ix)/VoxelMeters*256.0;
        const double YSubvoxels=Grid.Y(Iy)/VoxelMeters*256.0;
        const auto X=static_cast<std::int64_t>(std::llround(XSubvoxels));
        const auto Y=static_cast<std::int64_t>(std::llround(YSubvoxels));
        const double MinZ=static_cast<double>(Patch.OriginVoxels().Z);
        const double MaxZ=MinZ+Patch.ExtentVoxels();
        constexpr std::int64_t NormalProbeSubvoxels=256;
        for(const auto& Plan:Assembly.plan.strates)
        {
            VoxelForgeCore::Heightfield2D Field(Assembly.recipe,Plan);
            if(!Field.isValid())continue;
            const auto Column=Field.point(X,Y,0);
            if(Column.groundHeightSubvoxels>=Column.capHeightSubvoxels)continue;
            const auto XMinus=Field.point(X-NormalProbeSubvoxels,Y,0);
            const auto XPlus=Field.point(X+NormalProbeSubvoxels,Y,0);
            const auto YMinus=Field.point(X,Y-NormalProbeSubvoxels,0);
            const auto YPlus=Field.point(X,Y+NormalProbeSubvoxels,0);
            const auto Slope=[&](std::int64_t Minus,std::int64_t Plus)
            {return static_cast<double>(Plus-Minus)/(2.0*NormalProbeSubvoxels);};
            const double Ground=static_cast<double>(Column.groundHeightSubvoxels)/256.0*VoxelMeters;
            const double Cap=static_cast<double>(Column.capHeightSubvoxels)/256.0*VoxelMeters;
            const double PlanMinZ=static_cast<double>(Plan.bottomChunkZ)*32.0*VoxelMeters;
            const double PlanMaxZ=static_cast<double>(Plan.topChunkZ+1)*32.0*VoxelMeters;
            if(Ground>=MinZ*VoxelMeters-0.25*VoxelMeters
                &&Ground<=MaxZ*VoxelMeters+0.25*VoxelMeters
                &&Ground>=PlanMinZ-0.25*VoxelMeters&&Ground<=PlanMaxZ+0.25*VoxelMeters)
                Roots.push_back({Ground,Slope(XMinus.groundHeightSubvoxels,
                    XPlus.groundHeightSubvoxels),Slope(YMinus.groundHeightSubvoxels,
                    YPlus.groundHeightSubvoxels)});
            if(Cap>=MinZ*VoxelMeters-0.25*VoxelMeters
                &&Cap<=MaxZ*VoxelMeters+0.25*VoxelMeters
                &&Cap>=PlanMinZ-0.25*VoxelMeters&&Cap<=PlanMaxZ+0.25*VoxelMeters)
                Roots.push_back({Cap,Slope(XMinus.capHeightSubvoxels,
                    XPlus.capHeightSubvoxels),Slope(YMinus.capHeightSubvoxels,
                    YPlus.capHeightSubvoxels)});
        }
        std::sort(Roots.begin(),Roots.end(),[](const FSurfaceRoot& A,const FSurfaceRoot& B)
            {return A.ZMeters<B.ZMeters;});
        Roots.erase(std::unique(Roots.begin(),Roots.end(),[](const FSurfaceRoot& A,
            const FSurfaceRoot& B){return std::abs(A.ZMeters-B.ZMeters)<0.01;}),Roots.end());
        return Roots;
    }

    inline bool MatchNearest(const std::vector<FMeshSurfaceHit>& Hits,const std::vector<bool>& Used,
        double Root,std::size_t& OutIndex,double& OutError)
    {
        double Best=std::numeric_limits<double>::infinity();std::size_t BestIndex=Hits.size();
        for(std::size_t I=0;I<Hits.size();++I)
            if(I>=Used.size()||!Used[I])
            {
                const double Error=std::abs(Hits[I].Z-Root);
                if(Error<Best){Best=Error;BestIndex=I;}
            }
        if(BestIndex==Hits.size())return false;
        OutIndex=BestIndex;OutError=Best;return true;
    }

    inline FErrorRow CompareToExact(const FHeightGrid& MeshGrid,const FVoxelTileKey& Patch,
        double VoxelMeters,double PixelBudgetMeters,
        const VoxelForgeCore::WorldAssemblyManifest& Assembly,const FViewBasis& View)
    {
        FErrorRow Row;Row.Lod=Patch.Level;Row.Valid=true;
        FErrorAccumulator Accumulator;
        for(int Iy=0;Iy<MeshGrid.Side;++Iy)
        for(int Ix=0;Ix<MeshGrid.Side;++Ix)
        {
            const auto Roots=ExactHeightRoots(MeshGrid,Ix,Iy,Patch,VoxelMeters,Assembly);
            if(Roots.empty()&&MeshGrid.At(Ix,Iy).empty())continue;
            std::vector<bool> Used(MeshGrid.At(Ix,Iy).size(),false);
            for(const FSurfaceRoot& Root:Roots)
            {
                std::size_t HitIndex=0;double ErrorMeters=0.0;
                if(!MatchNearest(MeshGrid.At(Ix,Iy),Used,Root.ZMeters,HitIndex,ErrorMeters))
                {
                    ++Row.Missing;
                    AddMissing(Accumulator);
                    continue;
                }
                Used[HitIndex]=true;
                AddError(Accumulator,MeshGrid.At(Ix,Iy)[HitIndex].Z-Root.ZMeters,
                    Root,PixelBudgetMeters,View,MeshGrid.X(Ix),MeshGrid.Y(Iy));
                ++Row.Samples;
            }
            // Heightfield2D declares ground and cap roots only. Other mesh intersections may
            // be valid tunnel geometry; the exact-mesh comparison checks extra pieces.
        }
        FinalizeError(Row,Accumulator);
        return Row;
    }

    inline FErrorRow CompareToReference(const FHeightGrid& MeshGrid,
        const FHeightGrid& ReferenceGrid,const FVoxelTileKey& Patch,double VoxelMeters,
        double PixelBudgetMeters,const VoxelForgeCore::WorldAssemblyManifest& Assembly,
        const FViewBasis& View)
    {
        FErrorRow Row;Row.Lod=Patch.Level;Row.Valid=true;
        FErrorAccumulator Accumulator;
        for(int Iy=0;Iy<ReferenceGrid.Side;++Iy)
        for(int Ix=0;Ix<ReferenceGrid.Side;++Ix)
        {
            const auto& Reference=ReferenceGrid.At(Ix,Iy);
            const auto& Mesh=MeshGrid.At(Ix,Iy);
            if(Reference.empty()&&Mesh.empty())continue;
            const auto Roots=ExactHeightRoots(ReferenceGrid,Ix,Iy,Patch,VoxelMeters,Assembly);
            std::vector<bool> UsedReference(Reference.size(),false),UsedMesh(Mesh.size(),false);
            for(std::size_t ReferenceIndex=0;ReferenceIndex<Reference.size();++ReferenceIndex)
            {
                const double ReferenceZ=Reference[ReferenceIndex].Z;
                UsedReference[ReferenceIndex]=true;
                std::size_t MeshIndex=0;double MeshError=0.0;
                if(!MatchNearest(Mesh,UsedMesh,ReferenceZ,MeshIndex,MeshError))
                {
                    ++Row.Missing;AddMissing(Accumulator);continue;
                }
                UsedMesh[MeshIndex]=true;
                FSurfaceRoot Root{ReferenceZ,0.0,0.0};
                if(!Roots.empty())
                {
                    const auto BestRoot=std::min_element(Roots.begin(),Roots.end(),
                        [ReferenceZ](const FSurfaceRoot& A,const FSurfaceRoot& B)
                        {return std::abs(A.ZMeters-ReferenceZ)<std::abs(B.ZMeters-ReferenceZ);});
                    if(BestRoot!=Roots.end()&&std::abs(BestRoot->ZMeters-ReferenceZ)<1.0)
                        Root=*BestRoot;
                }
                AddError(Accumulator,Mesh[MeshIndex].Z-Reference[ReferenceIndex].Z,
                    Root,PixelBudgetMeters,View,ReferenceGrid.X(Ix),ReferenceGrid.Y(Iy));
            }
            for(std::size_t I=0;I<Mesh.size();++I)
                if(!UsedMesh[I]){++Row.Floating;AddMissing(Accumulator);}
        }
        FinalizeError(Row,Accumulator);
        return Row;
    }

    struct FMeshComparisonRow
    {
        FErrorRow Geometry;
        std::uint64_t Holes=0,Cracks=0,Floating=0,MaterialMismatches=0;
        double ShadeP50=0.0,ShadeP99=0.0,ShadeMaximum=0.0;
        double NormalDegreesP50=0.0,NormalDegreesP99=0.0,NormalDegreesMaximum=0.0;
        std::vector<double> GeometryPixels,ShadeValues,NormalAngles;
        bool Valid=false;
    };

    inline std::array<double,3> UnitNormal(const FMeshSurfaceHit& Hit)
    {
        std::array<double,3> N=Hit.Normal;
        const double Length=std::sqrt(N[0]*N[0]+N[1]*N[1]+N[2]*N[2]);
        if(Length<=1.0e-12)return {0.0,0.0,Hit.Ceiling?-1.0:1.0};
        for(double& Component:N)Component/=Length;
        return N;
    }

    inline FMeshComparisonRow CompareMeshes(const FHeightGrid& Reference,
        const FHeightGrid& Candidate,const FVoxelTileKey& Patch,double VoxelMeters,
        double PixelBudgetMeters,const FViewBasis& View)
    {
        FMeshComparisonRow Row;Row.Valid=Reference.Side==Candidate.Side;
        if(!Row.Valid)return Row;
        FErrorAccumulator Geometry;
        std::vector<double> ShadeValues,NormalAngles;
        const std::size_t CellCount=static_cast<std::size_t>(Reference.Side)*Reference.Side;
        ShadeValues.reserve(CellCount);NormalAngles.reserve(CellCount);
        const double LightLength=std::sqrt(0.3*0.3+(-0.4)*(-0.4)+0.85*0.85);
        const std::array<double,3> Light={0.3/LightLength,-0.4/LightLength,0.85/LightLength};
        for(int Iy=0;Iy<Reference.Side;++Iy)
        for(int Ix=0;Ix<Reference.Side;++Ix)
        {
            const auto& Before=Reference.At(Ix,Iy);
            const auto& After=Candidate.At(Ix,Iy);
            std::vector<bool> Used(After.size(),false);
            for(std::size_t I=0;I<Before.size();++I)
            {
                std::size_t Matched=0;double Difference=0.0;
                if(!MatchNearest(After,Used,Before[I].Z,Matched,Difference))
                {
                    ++Row.Holes;AddMissing(Geometry);ShadeValues.push_back(1.0);
                    NormalAngles.push_back(180.0);continue;
                }
                Used[Matched]=true;
                const FMeshSurfaceHit& A=Before[I];
                const FMeshSurfaceHit& B=After[Matched];
                const auto NA=UnitNormal(A),NB=UnitNormal(B);
                const double Dot=std::clamp(NA[0]*NB[0]+NA[1]*NB[1]+NA[2]*NB[2],-1.0,1.0);
                constexpr double Pi=3.14159265358979323846;
                const double Angle=std::acos(Dot)*180.0/Pi;
                NormalAngles.push_back(Angle);
                const double DiffuseA=std::max(0.0,NA[0]*Light[0]+NA[1]*Light[1]+NA[2]*Light[2]);
                const double DiffuseB=std::max(0.0,NB[0]*Light[0]+NB[1]*Light[1]+NB[2]*Light[2]);
                const bool bMaterialChanged=A.Color[0]!=B.Color[0]||A.Color[3]!=B.Color[3]
                    ||A.Ceiling!=B.Ceiling;
                if(bMaterialChanged)++Row.MaterialMismatches;
                const double MaterialDelta=bMaterialChanged?1.0:std::max(
                    std::abs(static_cast<double>(A.Color[1])-B.Color[1]),
                    std::abs(static_cast<double>(A.Color[2])-B.Color[2]))/255.0;
                ShadeValues.push_back(std::max(std::abs(DiffuseA-DiffuseB),MaterialDelta));
                const double NX=NA[2]<0.0?NA[0]/std::max(1.0e-12,-NA[2])
                    :-NA[0]/std::max(1.0e-12,NA[2]);
                const double NY=NA[2]<0.0?NA[1]/std::max(1.0e-12,-NA[2])
                    :-NA[1]/std::max(1.0e-12,NA[2]);
                AddError(Geometry,B.Z-A.Z,{A.Z,NX,NY},PixelBudgetMeters,View,
                    Reference.X(Ix),Reference.Y(Iy));
            }
            for(std::size_t I=0;I<After.size();++I)
                if(!Used[I])
                {++Row.Floating;AddMissing(Geometry);ShadeValues.push_back(1.0);NormalAngles.push_back(180.0);}
        }
        for(int Iy=0;Iy<Reference.Side;++Iy)
        for(int Ix=0;Ix<Reference.Side;++Ix)
        {
            const std::size_t Index=static_cast<std::size_t>(Iy)*Reference.Side+Ix;
            const auto CheckEdge=[&](int NX,int NY)
            {
                if(NX>=Reference.Side||NY>=Reference.Side)return;
                const auto& A=Reference.Hits[Index];
                const auto& B=Reference.At(NX,NY);
                const auto& CA=Candidate.Hits[Index];
                const auto& CB=Candidate.At(NX,NY);
                if(!A.empty()&&!B.empty()&&(CA.empty()!=CB.empty()))++Row.Cracks;
            };
            CheckEdge(Ix+1,Iy);CheckEdge(Ix,Iy+1);
        }
        Row.Geometry.Lod=Patch.Level;Row.Geometry.Valid=Row.Valid;
        Row.Geometry.Missing=Row.Holes;Row.Geometry.Floating=Row.Floating;
        FinalizeError(Row.Geometry,Geometry);
        Row.ShadeP50=Percentile(ShadeValues,0.50);
        Row.ShadeP99=Percentile(ShadeValues,0.99);
        Row.ShadeMaximum=ShadeValues.empty()?0.0:*std::max_element(ShadeValues.begin(),ShadeValues.end());
        Row.NormalDegreesP50=Percentile(NormalAngles,0.50);
        Row.NormalDegreesP99=Percentile(NormalAngles,0.99);
        Row.NormalDegreesMaximum=NormalAngles.empty()?0.0
            :*std::max_element(NormalAngles.begin(),NormalAngles.end());
        Row.GeometryPixels=std::move(Geometry.Pixels);
        Row.ShadeValues=std::move(ShadeValues);
        Row.NormalAngles=std::move(NormalAngles);
        return Row;
    }

    inline FPopRow ComparePop(const FHeightGrid& Coarse,const FHeightGrid& Fine,
        const FVoxelTileKey& Patch,double VoxelMeters,double PixelBudgetMeters,
        const VoxelForgeCore::WorldAssemblyManifest& Assembly,const FViewBasis& View)
    {
        FPopRow Row;Row.CoarseLod=Patch.Level;Row.Valid=true;
        FErrorAccumulator Accumulator;
        for(int Iy=0;Iy<Coarse.Side;++Iy)
        for(int Ix=0;Ix<Coarse.Side;++Ix)
        {
            const auto& CoarseHits=Coarse.At(Ix,Iy);
            const auto& FineHits=Fine.At(Ix,Iy);
            const auto Roots=ExactHeightRoots(Coarse,Ix,Iy,Patch,VoxelMeters,Assembly);
            if(Roots.empty()&&CoarseHits.empty()&&FineHits.empty())continue;
            std::vector<bool> CoarseUsed(CoarseHits.size(),false);
            std::vector<bool> FineUsed(FineHits.size(),false);
            for(const FSurfaceRoot& Root:Roots)
            {
                std::size_t CoarseIndex=0,FineIndex=0;
                double CoarseError=0.0,FineError=0.0;
                if(!MatchNearest(CoarseHits,CoarseUsed,Root.ZMeters,CoarseIndex,CoarseError)
                    ||!MatchNearest(FineHits,FineUsed,Root.ZMeters,FineIndex,FineError))
                {
                    ++Row.Cracks;AddMissing(Accumulator);continue;
                }
                CoarseUsed[CoarseIndex]=true;
                FineUsed[FineIndex]=true;
                AddError(Accumulator,CoarseHits[CoarseIndex].Z-FineHits[FineIndex].Z,
                    Root,PixelBudgetMeters,View,Coarse.X(Ix),Coarse.Y(Iy));
                ++Row.Samples;
            }
        }
        Row.SamplePixels=Accumulator.Pixels;
        FErrorRow Stats;Stats.Lod=Patch.Level;Stats.Valid=true;
        FinalizeError(Stats,Accumulator);
        Row.Samples=Stats.Samples;Row.P50=Stats.P50;Row.P99=Stats.P99;Row.Maximum=Stats.Maximum;
        Row.VerticalP50=Stats.VerticalP50;Row.VerticalP99=Stats.VerticalP99;
        Row.VerticalMaximum=Stats.VerticalMaximum;
        Row.MaximumNormalDistanceMeters=Stats.MaximumNormalDistanceMeters;
        Row.MaximumX=Stats.MaximumX;Row.MaximumY=Stats.MaximumY;
        Row.MaximumReferenceZ=Stats.MaximumReferenceZ;
        Row.MaximumDisplayedZ=Stats.MaximumDisplayedZ;
        Row.MaximumSlopeDegrees=Stats.MaximumSlopeDegrees;
        Row.SlopeSamples=Stats.SlopeSamples;Row.SlopeP99=Stats.SlopeP99;
        Row.SlopeMaximum=Stats.SlopeMaximum;
        return Row;
    }

    inline bool PopNotWorseAtAnySample(const FPopRow& Before,const FPopRow& After)
    {
        if(!Before.Valid||!After.Valid||Before.Samples<16
            ||Before.SamplePixels.size()!=After.SamplePixels.size()
            ||After.Cracks>Before.Cracks||After.Floating>Before.Floating)return false;
        for(std::size_t I=0;I<Before.SamplePixels.size();++I)
            if(After.SamplePixels[I]>Before.SamplePixels[I]+1.0e-4)return false;
        return true;
    }

    inline std::array<double,9> ClosestDisplayedDistancesMeters(
        const FStandaloneOwnerWorld& World,const FVector& CameraVoxels,
        std::uint32_t VoxelSizeMicrons)
    {
        const double VoxelMeters=static_cast<double>(VoxelSizeMicrons)/1000000.0;
        std::array<double,9> Minimums{};
        Minimums.fill(std::numeric_limits<double>::infinity());
        for(int Lod=0;Lod<=8;++Lod)
        {
            const int PhaseCount=1<<Lod;
            for(int Phase=0;Phase<PhaseCount;++Phase)
            {
                FVector Probe=CameraVoxels;
                Probe.X=static_cast<float>(Phase*32.0+32.0-1.0/256.0);
                const FStandaloneDesiredTileSet Desired=StandaloneWorldMesher::BuildDesiredSet(World,Probe);
                for(const FVoxelTileKey& Tile:Desired.Tiles)
                {
                    if(Tile.Level!=Lod)continue;
                    Minimums[Lod]=std::min(Minimums[Lod],
                        TileAabbDistanceMeters(Tile,Probe,VoxelMeters));
                }
            }
        }
        for(double& Distance:Minimums)
            if(!std::isfinite(Distance))Distance=0.0;
        return Minimums;
    }

    inline std::array<double,9> ClipmapClosestDisplayedDistancesMeters(
        const FStandaloneOwnerWorld& World,std::uint32_t VoxelSizeMicrons)
    {
        const double VoxelMeters=static_cast<double>(VoxelSizeMicrons)/1000000.0;
        std::array<double,9> Distances{};
        for(int Lod=1;Lod<=8;++Lod)
            Distances[Lod]=static_cast<double>(World.ClipRadius)
                *static_cast<double>(CHUNK_SIZE)*(1ull<<(Lod-1))*VoxelMeters;
        return Distances;
    }

    inline void PrintSlopeTail(const char* Kind,int Lod,const std::array<std::uint64_t,4>& Counts,
        const std::array<double,4>& P99,const std::array<double,4>& Maximum)
    {
        constexpr const char* Labels[4]={"lt30","30to45","45to60","ge60"};
        for(std::size_t Bin=0;Bin<Counts.size();++Bin)
            std::printf("[ScreenSpaceSlopeTail] kind=%s lod=%d slope_deg=%s samples=%llu p99_px=%.4f max_px=%.4f\n",
                Kind,Lod,Labels[Bin],static_cast<unsigned long long>(Counts[Bin]),P99[Bin],Maximum[Bin]);
    }

    inline void PrintKnownFailureControl(double CellMeters,double DisplayDistanceMeters,
        const FViewBasis& View)
    {
        constexpr double Slope=100.0;
        const double LateralMeters=CellMeters/200.0;
        const double VerticalMeters=Slope*LateralMeters;
        const double NormalLength=std::sqrt(1.0+Slope*Slope);
        const double NormalDistance=VerticalMeters/NormalLength;
        const std::array<double,3> Normal={-Slope/NormalLength,0.0,1.0/NormalLength};
        const double ProjectedNormal=std::hypot(
            Normal[0]*View.Right[0]+Normal[1]*View.Right[1]+Normal[2]*View.Right[2],
            Normal[0]*View.Up[0]+Normal[1]*View.Up[1]+Normal[2]*View.Up[2]);
        const double PixelBudgetMeters=DisplayDistanceMeters*OnePixelRadians;
        const double VerticalPixels=VerticalMeters/PixelBudgetMeters;
        const double NormalPixels=NormalDistance*ProjectedNormal/PixelBudgetMeters;
        const bool Passed=VerticalPixels>4.0&&NormalPixels<1.0;
        std::printf("[ScreenSpaceKnownFailureControl] case=steep-plane-lateral-offset slope_deg=%.4f cell_m=%.6f lateral_m=%.6f cell_fraction=%.6f vertical_error_m=%.6f vertical_height_px=%.4f closest_normal_m=%.6f projected_normal_px=%.4f status=%s\n",
            std::atan(Slope)*180.0/3.14159265358979323846,CellMeters,LateralMeters,
            LateralMeters/CellMeters,VerticalMeters,VerticalPixels,NormalDistance,NormalPixels,
            Passed?"EXPECTED_VERTICAL_FALSE_POSITIVE":"CONTROL_FAILED");
    }

    inline bool Run(const std::array<std::optional<FStandaloneCpuTile>,9>& SampleTiles,
        FStandaloneStreamingGenerator& Streamer,
        const std::vector<FVoxelTileKey>& DesiredKeys,
        const VoxelForgeCore::WorldAssemblyManifest& Assembly,
        const std::array<double,9>& Distances,std::uint32_t VoxelSizeMicrons,int MaxLod,
        const std::array<std::int64_t,9>& PixelBudgetsSubvoxels,
        double CameraYawDegrees,double CameraPitchDegrees)
    {
        const double VoxelMeters=static_cast<double>(VoxelSizeMicrons)/1000000.0;
        const FViewBasis View=MakeViewBasis(CameraYawDegrees,CameraPitchDegrees);
        std::array<double,9> Budgets{};
        VoxelForgeCore::TileGenerationOptions BaselineOptions;
        BaselineOptions.overrideFarDisplayMesh=true;
        VoxelForgeCore::TileGenerationOptions CandidateOptions;
        CandidateOptions.overrideFarDisplayMesh=true;
        CandidateOptions.farDisplayGeometryBudgetMillipixels=Assembly.recipe.farDisplayGeometryBudgetMillipixels;
        CandidateOptions.farRefinementGridScale=Assembly.recipe.farRefinementGridScale;
        CandidateOptions.farRefinementLodMask=Assembly.recipe.farRefinementLodMask;
        CandidateOptions.farConnectorStationScale=Assembly.recipe.farConnectorStationScale;
        CandidateOptions.farConnectorRayMask=Assembly.recipe.farConnectorRayMask;
        CandidateOptions.farConnectorReuseStationRadii=Assembly.recipe.farConnectorReuseStationRadii;
        bool GatePass=true;
        constexpr double ShadeLimit=1.0/255.0;
        std::printf("[ScreenSpaceGate] angular_pixel=%.9f rad (horizontal_fov=90deg,width=2560) lod>=3_geometry_limit=1.0px same_lod_reference=today mesh sampling=65x65 skirts=excluded shading=max(lambert-delta,material-selector-or-blend-delta) shade_limit=1/255 normal_probe_vox=1 view_yaw_deg=%.3f view_pitch_deg=%.3f\n",
            OnePixelRadians,CameraYawDegrees,CameraPitchDegrees);
        for(int Lod=0;Lod<=MaxLod&&Lod<9;++Lod)
        {
            Budgets[Lod]=Distances[Lod]*OnePixelRadians;
            std::printf("[ScreenSpaceDistance] lod=%d closest_m=%.3f pixel_budget_m=%.6f pixel_budget_cm=%.3f pixel_budget_subvoxels=%lld rule=%s\n",
                Lod,Distances[Lod],Budgets[Lod],Budgets[Lod]*100.0,
                static_cast<long long>(PixelBudgetsSubvoxels[Lod]),Lod<3?"exact-2cm":"same_lod_1px");
        }
        const double LOD3CellMeters=static_cast<double>(CHUNK_SIZE*(1<<3))*VoxelMeters/64.0;
        PrintKnownFailureControl(LOD3CellMeters,Distances[3],View);
        VoxelForgeCore::Heightfield2D::resetThreadColumnCacheStats(false,true,0);
        for(int Lod=3;Lod<=MaxLod&&Lod<9;++Lod)
        {
            if(!SampleTiles[Lod]||SampleTiles[Lod]->Vertices.empty())
            {
                std::printf("[ScreenSpaceGateLOD] lod=%d status=NO_SURFACE_TILE\n",Lod);
                GatePass=false;continue;
            }
            const FStandaloneCpuTile& CoarseTile=*SampleTiles[Lod];
            const FStandaloneCpuTile* CoarsePtr=&CoarseTile;
            const std::span<const FStandaloneCpuTile* const> CoarseSpan(&CoarsePtr,1);
            const FHeightGrid CoarseGrid=Rasterize(CoarseSpan,CoarseTile.Key,VoxelMeters);
            FStandaloneCpuTile BaselineTile=Streamer.GenerateDiagnostic(
                CoarseTile.Key,0,1,true,&BaselineOptions,true);
            const FStandaloneCpuTile* BaselinePointer=&BaselineTile;
            const std::span<const FStandaloneCpuTile* const> BaselineSpan(&BaselinePointer,1);
            const FHeightGrid BaselineGrid=Rasterize(BaselineSpan,BaselineTile.Key,VoxelMeters);
            const FMeshComparisonRow MeshDelta=CompareMeshes(BaselineGrid,CoarseGrid,
                CoarseTile.Key,VoxelMeters,Budgets[Lod],View);
            const FErrorRow& Error=MeshDelta.Geometry;
            const bool LodPass=Error.Valid&&Error.Missing==0&&Error.Floating==0
                &&MeshDelta.Cracks==0&&Error.Samples>=16&&Error.Maximum<=1.0
                &&MeshDelta.ShadeMaximum<=ShadeLimit;
            GatePass=GatePass&&LodPass;
            std::printf("[ScreenSpaceGateLOD] lod=%d key=(%d,%d,%d) samples=%llu geometry_px_p50=%.4f p99=%.4f max=%.4f shading_delta_p50=%.6f p99=%.6f max=%.6f normal_deg_p50=%.4f p99=%.4f max=%.4f material_mismatch=%llu holes=%llu cracks=%llu floating=%llu status=%s\n",
                Lod,CoarseTile.Key.Coord.X,CoarseTile.Key.Coord.Y,CoarseTile.Key.Coord.Z,
                static_cast<unsigned long long>(Error.Samples),Error.P50,Error.P99,Error.Maximum,
                MeshDelta.ShadeP50,MeshDelta.ShadeP99,MeshDelta.ShadeMaximum,
                MeshDelta.NormalDegreesP50,MeshDelta.NormalDegreesP99,MeshDelta.NormalDegreesMaximum,
                static_cast<unsigned long long>(MeshDelta.MaterialMismatches),
                static_cast<unsigned long long>(MeshDelta.Holes),
                static_cast<unsigned long long>(MeshDelta.Cracks),
                static_cast<unsigned long long>(MeshDelta.Floating),
                LodPass?"PASS":"REJECT");

            const FErrorRow FieldError=CompareToExact(CoarseGrid,CoarseTile.Key,VoxelMeters,
                Budgets[Lod],Assembly,View);
            const bool FieldPass=FieldError.Valid&&FieldError.Missing==0&&FieldError.Floating==0
                &&FieldError.Samples>=16&&FieldError.Maximum<=1.0;
            std::printf("[ScreenSpaceGateFieldLOD] lod=%d samples=%llu field_geometry_px_p50=%.4f p99=%.4f max=%.4f missing=%llu floating=%llu status=INFO_ONLY\n",
                Lod,static_cast<unsigned long long>(FieldError.Samples),FieldError.P50,
                FieldError.P99,FieldError.Maximum,
                static_cast<unsigned long long>(FieldError.Missing),
                static_cast<unsigned long long>(FieldError.Floating));

            bool FaultControlsPass=true;
            const FHeightGrid MissingGrid(BaselineGrid.Side,BaselineGrid.MinimumX,BaselineGrid.MaximumX,
                BaselineGrid.MinimumY,BaselineGrid.MaximumY);
            const FMeshComparisonRow MissingError=CompareMeshes(BaselineGrid,MissingGrid,
                CoarseTile.Key,VoxelMeters,Budgets[Lod],View);
            const bool MissingRejected=MissingError.Holes>0&&MissingError.Geometry.Maximum>1.0;
            FaultControlsPass=FaultControlsPass&&MissingRejected;
            std::printf("[ScreenSpaceFaultInjection] lod=%d case=missing-displayed-surface holes=%llu status=%s\n",
                Lod,static_cast<unsigned long long>(MissingError.Holes),
                MissingRejected?"REJECTED_AS_REQUIRED":"GATE_FALSE_NEGATIVE");

            FStandaloneCpuTile FloatingTile=BaselineTile;
            const std::uint32_t FloatingVertexBase=static_cast<std::uint32_t>(FloatingTile.Vertices.size());
            for(FStandaloneCpuVertex Vertex:BaselineTile.Vertices)
            {
                Vertex.Position[2]+=2.0f;
                FloatingTile.Vertices.push_back(Vertex);
            }
            FloatingTile.VertexColors.insert(FloatingTile.VertexColors.end(),
                BaselineTile.VertexColors.begin(),BaselineTile.VertexColors.end());
            for(const std::uint32_t Index:BaselineTile.Indices)
                FloatingTile.Indices.push_back(FloatingVertexBase+Index);
            FloatingTile.TriangleSkirtFaceMasks.insert(FloatingTile.TriangleSkirtFaceMasks.end(),
                BaselineTile.TriangleSkirtFaceMasks.begin(),BaselineTile.TriangleSkirtFaceMasks.end());
            const FStandaloneCpuTile* FloatingPointer=&FloatingTile;
            const std::span<const FStandaloneCpuTile* const> FloatingSpan(&FloatingPointer,1);
            const FHeightGrid FloatingGrid=Rasterize(FloatingSpan,BaselineTile.Key,VoxelMeters);
            const FMeshComparisonRow FloatingError=CompareMeshes(BaselineGrid,FloatingGrid,
                CoarseTile.Key,VoxelMeters,Budgets[Lod],View);
            const bool FloatingRejected=FloatingError.Floating>0||FloatingError.Geometry.Maximum>1.0;
            FaultControlsPass=FaultControlsPass&&FloatingRejected;
            std::printf("[ScreenSpaceFaultInjection] lod=%d case=2m-floating-surface floating=%llu max_px=%.4f status=%s\n",
                Lod,static_cast<unsigned long long>(FloatingError.Floating),FloatingError.Geometry.Maximum,
                FloatingRejected?"REJECTED_AS_REQUIRED":"GATE_FALSE_NEGATIVE");

            std::vector<FStandaloneCpuTile> FineBeforeTiles,FineAfterTiles;
            FineBeforeTiles.reserve(8);FineAfterTiles.reserve(8);
            for(int Dz=0;Dz<=1;++Dz)
            for(int Dy=0;Dy<=1;++Dy)
            for(int Dx=0;Dx<=1;++Dx)
            {
                const FVoxelTileKey FineKey(FIntVector(2*CoarseTile.Key.Coord.X+Dx,
                    2*CoarseTile.Key.Coord.Y+Dy,2*CoarseTile.Key.Coord.Z+Dz),Lod-1);
                FStandaloneCpuTile FineBefore=Streamer.GenerateDiagnostic(
                    FineKey,0,1,true,&BaselineOptions,true);
                FStandaloneCpuTile FineAfter=Streamer.GenerateDiagnostic(
                    FineKey,0,1,true,nullptr,true);
                if(!FineBefore.Vertices.empty())FineBeforeTiles.push_back(std::move(FineBefore));
                if(!FineAfter.Vertices.empty())FineAfterTiles.push_back(std::move(FineAfter));
            }
            std::vector<const FStandaloneCpuTile*> FineBeforePointers,FineAfterPointers;
            FineBeforePointers.reserve(FineBeforeTiles.size());FineAfterPointers.reserve(FineAfterTiles.size());
            for(const auto& Tile:FineBeforeTiles)FineBeforePointers.push_back(&Tile);
            for(const auto& Tile:FineAfterTiles)FineAfterPointers.push_back(&Tile);
            const auto FineBeforeSpan=std::span<const FStandaloneCpuTile* const>(
                FineBeforePointers.data(),FineBeforePointers.size());
            const auto FineAfterSpan=std::span<const FStandaloneCpuTile* const>(
                FineAfterPointers.data(),FineAfterPointers.size());
            const FHeightGrid FineBeforeGrid=Rasterize(FineBeforeSpan,CoarseTile.Key,VoxelMeters);
            const FHeightGrid FineAfterGrid=Rasterize(FineAfterSpan,CoarseTile.Key,VoxelMeters);
            const FPopRow PopBefore=ComparePop(BaselineGrid,FineBeforeGrid,CoarseTile.Key,VoxelMeters,
                Budgets[Lod],Assembly,View);
            const FPopRow PopAfter=ComparePop(CoarseGrid,FineAfterGrid,CoarseTile.Key,VoxelMeters,
                Budgets[Lod],Assembly,View);
            const bool PopPass=PopNotWorseAtAnySample(PopBefore,PopAfter);
            GatePass=GatePass&&PopPass;
            std::printf("[ScreenSpaceSwap] swap=%d_to_%d distance_m=%.3f samples_before=%llu pop_before_px_p50=%.4f p99=%.4f max=%.4f pop_after_px_p50=%.4f p99=%.4f max=%.4f missing_root_samples_before=%llu after=%llu status=%s\n",
                Lod,Lod-1,Distances[Lod],static_cast<unsigned long long>(PopBefore.Samples),
                PopBefore.P50,PopBefore.P99,PopBefore.Maximum,PopAfter.P50,PopAfter.P99,PopAfter.Maximum,
                static_cast<unsigned long long>(PopBefore.Cracks),static_cast<unsigned long long>(PopAfter.Cracks),
                PopPass?"PASS":"REJECT");

            const FHeightGrid MissingFineGrid(BaselineGrid.Side,BaselineGrid.MinimumX,
                BaselineGrid.MaximumX,BaselineGrid.MinimumY,BaselineGrid.MaximumY);
            const FPopRow MissingFineError=ComparePop(BaselineGrid,MissingFineGrid,
                CoarseTile.Key,VoxelMeters,Budgets[Lod],Assembly,View);
            const bool CrackRejected=MissingFineError.Cracks>0&&MissingFineError.Maximum>1.0;
            FaultControlsPass=FaultControlsPass&&CrackRejected;
            std::printf("[ScreenSpaceFaultInjection] lod=%d case=missing-fine-ring-surface cracks=%llu max_px=%.4f status=%s\n",
                Lod,static_cast<unsigned long long>(MissingFineError.Cracks),
                MissingFineError.Maximum,CrackRejected?"REJECTED_AS_REQUIRED":"GATE_FALSE_NEGATIVE");

            bool NegativeControlsPass=true;
            const auto CheckNegative=[&](const char* Name,std::uint32_t Drop,std::uint32_t Stride)
            {
                FStandaloneCpuTile Bad=Streamer.GenerateDiagnostic(
                    CoarseTile.Key,Drop,Stride,true,nullptr,true);
                const FStandaloneCpuTile* BadPointer=&Bad;
                const std::span<const FStandaloneCpuTile* const> BadSpan(&BadPointer,1);
                const FHeightGrid BadGrid=Rasterize(BadSpan,Bad.Key,VoxelMeters);
                const FMeshComparisonRow BadError=CompareMeshes(BaselineGrid,BadGrid,
                    CoarseTile.Key,VoxelMeters,Budgets[Lod],View);
                const bool Rejected=BadError.Holes>0||BadError.Floating>0||BadError.Cracks>0
                    ||BadError.Geometry.Samples<16||BadError.Geometry.Maximum>1.0
                    ||BadError.ShadeMaximum>ShadeLimit;
                NegativeControlsPass=NegativeControlsPass&&Rejected;
                std::printf("[ScreenSpaceNegative] lod=%d case=%s p50_px=%.4f p99_px=%.4f max_px=%.4f shade_max=%.6f holes=%llu cracks=%llu floating=%llu status=%s\n",
                    Lod,Name,BadError.Geometry.P50,BadError.Geometry.P99,BadError.Geometry.Maximum,
                    BadError.ShadeMaximum,static_cast<unsigned long long>(BadError.Holes),
                    static_cast<unsigned long long>(BadError.Cracks),
                    static_cast<unsigned long long>(BadError.Floating),
                    Rejected?"REJECTED_AS_REQUIRED":"GATE_FALSE_NEGATIVE");
            };
            CheckNegative("s13-global-octave-drop",1,1);
            CheckNegative("stride-2-column-interpolation",0,2);
            GatePass=GatePass&&FaultControlsPass&&NegativeControlsPass;
        }

        struct FCoverageRow
        {
            std::uint64_t Tiles=0,Holes=0,Cracks=0,Floating=0,MaterialMismatches=0;
            std::vector<double> GeometryPixels,ShadeValues,NormalAngles;
        };
        std::array<FCoverageRow,9> Coverage;
        for(const FVoxelTileKey& Key:DesiredKeys)
        {
            const int Lod=Key.Level;
            if(Lod<3||Lod>MaxLod||Lod>=9)continue;
            FStandaloneCpuTile BaselineTile=Streamer.GenerateDiagnostic(
                Key,0,1,true,&BaselineOptions,true);
            FStandaloneCpuTile CandidateTile=Streamer.GenerateDiagnostic(
                Key,0,1,true,&CandidateOptions,true);
            const FStandaloneCpuTile* BaselinePointer=&BaselineTile;
            const FStandaloneCpuTile* CandidatePointer=&CandidateTile;
            const std::span<const FStandaloneCpuTile* const> BaselineSpan(&BaselinePointer,1);
            const std::span<const FStandaloneCpuTile* const> CandidateSpan(&CandidatePointer,1);
            const FHeightGrid BaselineGrid=Rasterize(BaselineSpan,Key,VoxelMeters);
            const FHeightGrid CandidateGrid=Rasterize(CandidateSpan,Key,VoxelMeters);
            const FMeshComparisonRow Delta=CompareMeshes(BaselineGrid,CandidateGrid,
                Key,VoxelMeters,Budgets[Lod],View);
            auto& Row=Coverage[Lod];
            ++Row.Tiles;
            Row.Holes+=Delta.Holes;Row.Cracks+=Delta.Cracks;Row.Floating+=Delta.Floating;
            Row.MaterialMismatches+=Delta.MaterialMismatches;
            Row.GeometryPixels.insert(Row.GeometryPixels.end(),
                Delta.GeometryPixels.begin(),Delta.GeometryPixels.end());
            Row.ShadeValues.insert(Row.ShadeValues.end(),Delta.ShadeValues.begin(),Delta.ShadeValues.end());
            Row.NormalAngles.insert(Row.NormalAngles.end(),Delta.NormalAngles.begin(),Delta.NormalAngles.end());
        }
        for(int Lod=3;Lod<=MaxLod&&Lod<9;++Lod)
        {
            auto& Row=Coverage[Lod];
            const double GeometryP50=Percentile(Row.GeometryPixels,0.50);
            const double GeometryP99=Percentile(Row.GeometryPixels,0.99);
            const double GeometryMaximum=Row.GeometryPixels.empty()?0.0
                :*std::max_element(Row.GeometryPixels.begin(),Row.GeometryPixels.end());
            const double ShadeP50=Percentile(Row.ShadeValues,0.50);
            const double ShadeP99=Percentile(Row.ShadeValues,0.99);
            const double ShadeMaximum=Row.ShadeValues.empty()?0.0
                :*std::max_element(Row.ShadeValues.begin(),Row.ShadeValues.end());
            const double NormalP50=Percentile(Row.NormalAngles,0.50);
            const double NormalP99=Percentile(Row.NormalAngles,0.99);
            const double NormalMaximum=Row.NormalAngles.empty()?0.0
                :*std::max_element(Row.NormalAngles.begin(),Row.NormalAngles.end());
            const bool LodPass=Row.Tiles>0&&Row.GeometryPixels.size()>=16
                &&GeometryMaximum<=1.0&&ShadeMaximum<=ShadeLimit
                &&Row.Holes==0&&Row.Cracks==0&&Row.Floating==0;
            GatePass=GatePass&&LodPass;
            std::printf("[ScreenSpaceGateCoverageLOD] lod=%d tiles=%llu geometry_samples=%llu geometry_px_p50=%.4f p99=%.4f max=%.4f shading_delta_p50=%.6f p99=%.6f max=%.6f normal_deg_p50=%.4f p99=%.4f max=%.4f material_mismatch=%llu holes=%llu cracks=%llu floating=%llu status=%s\n",
                Lod,static_cast<unsigned long long>(Row.Tiles),
                static_cast<unsigned long long>(Row.GeometryPixels.size()),
                GeometryP50,GeometryP99,GeometryMaximum,ShadeP50,ShadeP99,ShadeMaximum,
                NormalP50,NormalP99,NormalMaximum,
                static_cast<unsigned long long>(Row.MaterialMismatches),
                static_cast<unsigned long long>(Row.Holes),
                static_cast<unsigned long long>(Row.Cracks),
                static_cast<unsigned long long>(Row.Floating),LodPass?"PASS":"REJECT");
        }

        std::uint64_t RepairOpportunities=0,BaselineWallHoles=0,BaselineWallFloating=0;
        std::uint64_t CandidateWallRays=0,CandidateWallHoles=0,CandidateWallFloating=0;
        std::vector<std::uint64_t> BaselineHoleIds,CandidateHoleIds;
        std::vector<std::uint64_t> BaselineFloatingIds,CandidateFloatingIds;
        std::vector<std::uint64_t> CandidatePanelRemovalHoleIds;
        bool bRanPanelRemovalControl=false;
        for(const FVoxelTileKey& Key:DesiredKeys)
        {
            if(Key.Level<6||Key.Level>8)continue;
            VoxelForgeCore::TileGenerationOptions NoPanelBaseline=BaselineOptions;
            NoPanelBaseline.diagnosticRemoveFarConnectorPanels=true;
            NoPanelBaseline.diagnosticAuditFarConnectorCoverage=true;
            const FStandaloneCpuTile BaseNoPanels=Streamer.GenerateDiagnostic(Key,0,1,true,&NoPanelBaseline);
            RepairOpportunities+=BaseNoPanels.CoreDiagnostics.farConnectorAuditHoles;

            VoxelForgeCore::TileGenerationOptions AuditedBaseline=BaselineOptions;
            AuditedBaseline.diagnosticAuditFarConnectorCoverage=true;
            const FStandaloneCpuTile BaseRepaired=Streamer.GenerateDiagnostic(Key,0,1,true,&AuditedBaseline);
            BaselineWallHoles+=BaseRepaired.CoreDiagnostics.farConnectorAuditHoles;
            BaselineWallFloating+=BaseRepaired.CoreDiagnostics.farConnectorAuditFloating;
            BaselineHoleIds.insert(BaselineHoleIds.end(),
                BaseRepaired.CoreDiagnostics.farConnectorAuditHoleIds.begin(),
                BaseRepaired.CoreDiagnostics.farConnectorAuditHoleIds.end());
            BaselineFloatingIds.insert(BaselineFloatingIds.end(),
                BaseRepaired.CoreDiagnostics.farConnectorAuditFloatingIds.begin(),
                BaseRepaired.CoreDiagnostics.farConnectorAuditFloatingIds.end());

            VoxelForgeCore::TileGenerationOptions AuditedCandidate=CandidateOptions;
            AuditedCandidate.diagnosticAuditFarConnectorCoverage=true;
            const FStandaloneCpuTile CandidateRepaired=Streamer.GenerateDiagnostic(Key,0,1,true,&AuditedCandidate);
            CandidateWallRays+=CandidateRepaired.CoreDiagnostics.farConnectorAuditRays;
            CandidateWallHoles+=CandidateRepaired.CoreDiagnostics.farConnectorAuditHoles;
            CandidateWallFloating+=CandidateRepaired.CoreDiagnostics.farConnectorAuditFloating;
            CandidateHoleIds.insert(CandidateHoleIds.end(),
                CandidateRepaired.CoreDiagnostics.farConnectorAuditHoleIds.begin(),
                CandidateRepaired.CoreDiagnostics.farConnectorAuditHoleIds.end());
            CandidateFloatingIds.insert(CandidateFloatingIds.end(),
                CandidateRepaired.CoreDiagnostics.farConnectorAuditFloatingIds.begin(),
                CandidateRepaired.CoreDiagnostics.farConnectorAuditFloatingIds.end());

            if(!bRanPanelRemovalControl&&BaseNoPanels.CoreDiagnostics.farConnectorAuditHoles>0)
            {
                AuditedCandidate.diagnosticRemoveFarConnectorPanels=true;
                const FStandaloneCpuTile CandidateNoPanels=Streamer.GenerateDiagnostic(
                    Key,0,1,true,&AuditedCandidate);
                CandidatePanelRemovalHoleIds=CandidateNoPanels.CoreDiagnostics.farConnectorAuditHoleIds;
                bRanPanelRemovalControl=true;
            }
        }
        std::sort(BaselineHoleIds.begin(),BaselineHoleIds.end());
        std::sort(CandidateHoleIds.begin(),CandidateHoleIds.end());
        std::sort(BaselineFloatingIds.begin(),BaselineFloatingIds.end());
        std::sort(CandidateFloatingIds.begin(),CandidateFloatingIds.end());
        const auto IsSubset=[](const std::vector<std::uint64_t>& A,const std::vector<std::uint64_t>& B)
        {return std::includes(B.begin(),B.end(),A.begin(),A.end());};
        std::sort(CandidatePanelRemovalHoleIds.begin(),CandidatePanelRemovalHoleIds.end());
        const bool bNoNewHoles=IsSubset(CandidateHoleIds,BaselineHoleIds);
        const bool bNoNewFloating=IsSubset(CandidateFloatingIds,BaselineFloatingIds);
        const bool bPanelRemovalDetected=bRanPanelRemovalControl
            &&!IsSubset(CandidatePanelRemovalHoleIds,CandidateHoleIds);
        std::vector<std::uint64_t> NewHoleIds,NewFloatingIds,PanelRemovalNewHoleIds;
        std::set_difference(CandidateHoleIds.begin(),CandidateHoleIds.end(),
            BaselineHoleIds.begin(),BaselineHoleIds.end(),std::back_inserter(NewHoleIds));
        std::set_difference(CandidateFloatingIds.begin(),CandidateFloatingIds.end(),
            BaselineFloatingIds.begin(),BaselineFloatingIds.end(),std::back_inserter(NewFloatingIds));
        std::set_difference(CandidatePanelRemovalHoleIds.begin(),CandidatePanelRemovalHoleIds.end(),
            CandidateHoleIds.begin(),CandidateHoleIds.end(),std::back_inserter(PanelRemovalNewHoleIds));
        const bool bRepairRowsPass=RepairOpportunities>0&&bNoNewHoles&&bNoNewFloating
            &&bPanelRemovalDetected;
        GatePass=GatePass&&bRepairRowsPass;
        std::printf("[ScreenSpaceConnectorGate] lods=6..8 baseline_missing_ray_stations=%llu candidate_audited_rays=%llu baseline_holes=%llu candidate_holes=%llu new_holes=%llu baseline_floating=%llu candidate_floating=%llu new_floating=%llu panel_removal_new_holes=%llu status=%s\n",
            static_cast<unsigned long long>(RepairOpportunities),
            static_cast<unsigned long long>(CandidateWallRays),
            static_cast<unsigned long long>(BaselineWallHoles),
            static_cast<unsigned long long>(CandidateWallHoles),
            static_cast<unsigned long long>(NewHoleIds.size()),
            static_cast<unsigned long long>(BaselineWallFloating),
            static_cast<unsigned long long>(CandidateWallFloating),
            static_cast<unsigned long long>(NewFloatingIds.size()),
            static_cast<unsigned long long>(PanelRemovalNewHoleIds.size()),
            bRepairRowsPass?"PASS":"REJECT");
        std::printf("[ScreenSpaceGateResult] status=%s\n",GatePass?"PASS":"REJECT");
        return GatePass;
    }
}