#pragma once

#include "OwnerWorldLoader.h"
#include "StandaloneShapeWorld.h"
#include "VoxelTypes.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <limits>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>

struct FFarMeshSegmentIndex
{
    FFarMeshSegmentIndex(const FVoxelMeshData& InMesh,const FIntVector& InOrigin,int32 InExtent,
                         int32 InCells,double InStep)
        :Mesh(InMesh),Origin(InOrigin),Extent(InExtent),Cells(std::max(1,InCells)),Step(InStep),
         ByCell(static_cast<size_t>(Cells)*Cells*Cells),
         TriangleMarks(static_cast<size_t>(InMesh.Triangles.Num()/3),0)
    {
        for(int32 Triangle=0;Triangle+2<Mesh.Triangles.Num();Triangle+=3)
        {
            const int32 A=Mesh.Triangles[Triangle],B=Mesh.Triangles[Triangle+1],C=Mesh.Triangles[Triangle+2];
            if(A<0||B<0||C<0||A>=Mesh.Vertices.Num()||B>=Mesh.Vertices.Num()||C>=Mesh.Vertices.Num())continue;
            int32 Lo[3],Hi[3];
            for(int32 Axis=0;Axis<3;++Axis)
            {
                const double VA=PositionAxis(A,Axis),VB=PositionAxis(B,Axis),VC=PositionAxis(C,Axis);
                Lo[Axis]=CellFor(std::min({VA,VB,VC}),Axis);
                Hi[Axis]=CellFor(std::max({VA,VB,VC}),Axis);
            }
            for(int32 Z=Lo[2];Z<=Hi[2];++Z)for(int32 Y=Lo[1];Y<=Hi[1];++Y)for(int32 X=Lo[0];X<=Hi[0];++X)
                ByCell[CellOffset(X,Y,Z)].push_back(Triangle/3);
        }
    }

    double NearestHitDifference(const FVector& Start,const FVector& Direction,double Length,
                                double TargetDistance)
    {
        double BestDifference=std::numeric_limits<double>::infinity();
        ++QueryStamp;
        if(QueryStamp==0){std::fill(TriangleMarks.begin(),TriangleMarks.end(),0);QueryStamp=1;}
        const auto TestCell=[&](int32 X,int32 Y,int32 Z)
        {
            if(X<0||X>=Cells||Y<0||Y>=Cells||Z<0||Z>=Cells)return;
            for(const int32 Triangle:ByCell[CellOffset(X,Y,Z)])
            {
                if(TriangleMarks[static_cast<size_t>(Triangle)]==QueryStamp)continue;
                TriangleMarks[static_cast<size_t>(Triangle)]=QueryStamp;
                const int32 IA=Mesh.Triangles[Triangle*3],IB=Mesh.Triangles[Triangle*3+1],IC=Mesh.Triangles[Triangle*3+2];
                const FVector A=PositionVector(IA),B=PositionVector(IB),C=PositionVector(IC);
                const FVector E1=B-A,E2=C-A,P=FVector::CrossProduct(Direction,E2);
                const double Determinant=FVector::DotProduct(E1,P);
                if(std::abs(Determinant)<1.0e-8)continue;
                const double Inv=1.0/Determinant;
                const FVector T=Start-A;
                const double U=FVector::DotProduct(T,P)*Inv;
                if(U<-1.0e-5||U>1.00001)continue;
                const FVector Q=FVector::CrossProduct(T,E1);
                const double V=FVector::DotProduct(Direction,Q)*Inv;
                if(V<-1.0e-5||U+V>1.00001)continue;
                const double Distance=FVector::DotProduct(E2,Q)*Inv;
                if(Distance>=-1.0e-4&&Distance<=Length+1.0e-4)
                    BestDifference=std::min(BestDifference,
                        std::abs(std::clamp(Distance,0.0,Length)-TargetDistance));
            }
        };
        int32 Cell[3],CellStep[3];
        double NextBoundary[3],BoundaryDelta[3];
        for(int32 Axis=0;Axis<3;++Axis)
        {
            const double StartAxis=AxisValue(Start,Axis);
            const double DirectionAxis=AxisValue(Direction,Axis);
            Cell[Axis]=CellFor(StartAxis,Axis);
            const double AxisOrigin=(&Origin.X)[Axis];
            if(DirectionAxis>1.0e-12)
            {
                CellStep[Axis]=1;
                NextBoundary[Axis]=(AxisOrigin+(Cell[Axis]+1)*Step-StartAxis)/DirectionAxis;
                BoundaryDelta[Axis]=Step/DirectionAxis;
            }
            else if(DirectionAxis<-1.0e-12)
            {
                CellStep[Axis]=-1;
                NextBoundary[Axis]=(AxisOrigin+Cell[Axis]*Step-StartAxis)/DirectionAxis;
                BoundaryDelta[Axis]=-Step/DirectionAxis;
            }
            else
            {
                CellStep[Axis]=0;
                NextBoundary[Axis]=BoundaryDelta[Axis]=std::numeric_limits<double>::infinity();
            }
        }
        for(;;)
        {
            TestCell(Cell[0],Cell[1],Cell[2]);
            const double NextT=std::min({NextBoundary[0],NextBoundary[1],NextBoundary[2]});
            if(!std::isfinite(NextT)||NextT>Length+1.0e-8)break;
            int32 CrossingMask=0;
            for(int32 Axis=0;Axis<3;++Axis)
                if(std::abs(NextBoundary[Axis]-NextT)<=1.0e-8)CrossingMask|=1<<Axis;
            for(int32 Subset=1;Subset<8;++Subset)
            {
                if(Subset&~CrossingMask)continue;
                int32 X=Cell[0],Y=Cell[1],Z=Cell[2];
                if(Subset&1)X+=CellStep[0];
                if(Subset&2)Y+=CellStep[1];
                if(Subset&4)Z+=CellStep[2];
                TestCell(X,Y,Z);
            }
            for(int32 Axis=0;Axis<3;++Axis)
                if(CrossingMask&(1<<Axis))
                {Cell[Axis]+=CellStep[Axis];NextBoundary[Axis]+=BoundaryDelta[Axis];}
        }
        return BestDifference;
    }

private:
    int32 CellFor(double Value,int32 Axis)const
    {return std::clamp(static_cast<int32>(std::floor((Value-(&Origin.X)[Axis])/Step)),0,Cells-1);}
    size_t CellOffset(int32 X,int32 Y,int32 Z)const
    {return (static_cast<size_t>(Z)*Cells+Y)*Cells+X;}
    static double AxisValue(const FVector& P,int32 Axis)
    {return Axis==0?P.X:(Axis==1?P.Y:P.Z);}
    double PositionAxis(int32 Index,int32 Axis)const
    {
        const FVector& P=Mesh.Vertices[Index];
        return AxisValue(P,Axis)*0.04;
    }
    FVector PositionVector(int32 Index)const
    {
        const FVector& P=Mesh.Vertices[Index];return P*0.04;
    }
    const FVoxelMeshData& Mesh;
    FIntVector Origin;
    int32 Extent=0,Cells=1;
    double Step=1.0;
    std::vector<std::vector<int32>> ByCell;
    std::vector<uint32> TriangleMarks;
    uint32 QueryStamp=0;
};

inline bool ContainsFarTilePoint(const FVoxelTileKey& Key,const FVector& P)
{
    const FIntVector Origin=Key.OriginVoxels();const int32 Extent=Key.ExtentVoxels();
    return P.X>=Origin.X&&P.X<=Origin.X+Extent&&P.Y>=Origin.Y&&P.Y<=Origin.Y+Extent
        &&P.Z>=Origin.Z&&P.Z<=Origin.Z+Extent;
}

inline bool FindFarConnectorSurfaceRadius(const FStandaloneShapeWorld& Shape,const FVector& Center,
    const FVector& Direction,double PreferredRadius,double& OutRadius)
{
    const double SearchLength=PreferredRadius*2.0+32.0;
    const double AcceptBand=std::max(16.0,PreferredRadius*0.75);
    constexpr double Stride=8.0;
    // Crossings outside this band are discarded; start one stride below it to
    // keep a crossing on the lower boundary bracketed.
    const double SearchMin=std::max(0.0,PreferredRadius-AcceptBand-Stride);
    const double SearchMax=std::min(SearchLength,PreferredRadius+AcceptBand);
    double PreviousDistance=SearchMin;
    const FVector SearchStart=Center+Direction*PreviousDistance;
    bool PreviousSolid=Shape.DensityAt(static_cast<float>(SearchStart.X),static_cast<float>(SearchStart.Y),
        static_cast<float>(SearchStart.Z))<0.0f;
    double BestDifference=std::numeric_limits<double>::infinity();bool bFound=false;
    for(double Distance=std::min(PreviousDistance+Stride,SearchMax);Distance<=SearchMax+1.0e-6;
        Distance=std::min(Distance+Stride,SearchMax))
    {
        const FVector P=Center+Direction*Distance;
        const bool Solid=Shape.DensityAt(static_cast<float>(P.X),static_cast<float>(P.Y),
            static_cast<float>(P.Z))<0.0f;
        if(Solid!=PreviousSolid)
        {
            double Lo=PreviousDistance,Hi=Distance;
            for(int32 I=0;I<6;++I)
            {
                const double Mid=(Lo+Hi)*0.5;const FVector M=Center+Direction*Mid;
                const bool MidSolid=Shape.DensityAt(static_cast<float>(M.X),static_cast<float>(M.Y),
                    static_cast<float>(M.Z))<0.0f;
                if(MidSolid==PreviousSolid)Lo=Mid;else Hi=Mid;
            }
            const double Crossing=(Lo+Hi)*0.5,Difference=std::abs(Crossing-PreferredRadius);
            if(Difference<=AcceptBand&&Difference<BestDifference)
            {OutRadius=Crossing;BestDifference=Difference;bFound=true;}
        }
        if(Distance>=SearchMax-1.0e-6)break;
        PreviousDistance=Distance;PreviousSolid=Solid;
    }
    if(!bFound||BestDifference>8.0)
    {
        const double LocalMin=std::max(0.0,PreferredRadius-8.0);
        const double LocalMax=std::min(SearchLength,PreferredRadius+8.0);
        double LocalPrevious=LocalMin;
        bool LocalPreviousSolid=Shape.DensityAt(static_cast<float>(Center.X+Direction.X*LocalPrevious),
            static_cast<float>(Center.Y+Direction.Y*LocalPrevious),
            static_cast<float>(Center.Z+Direction.Z*LocalPrevious))<0.0f;
        for(double Distance=std::min(LocalMin+1.0,LocalMax);Distance<=LocalMax+1.0e-6;
            Distance=std::min(Distance+1.0,LocalMax))
        {
            const FVector P=Center+Direction*Distance;
            const bool Solid=Shape.DensityAt(static_cast<float>(P.X),static_cast<float>(P.Y),
                static_cast<float>(P.Z))<0.0f;
            if(Solid!=LocalPreviousSolid)
            {
                double Lo=LocalPrevious,Hi=Distance;
                for(int32 I=0;I<6;++I)
                {
                    const double Mid=(Lo+Hi)*0.5;const FVector M=Center+Direction*Mid;
                    const bool MidSolid=Shape.DensityAt(static_cast<float>(M.X),static_cast<float>(M.Y),
                        static_cast<float>(M.Z))<0.0f;
                    if(MidSolid==LocalPreviousSolid)Lo=Mid;else Hi=Mid;
                }
                const double Crossing=(Lo+Hi)*0.5,Difference=std::abs(Crossing-PreferredRadius);
                if(Difference<=AcceptBand&&Difference<BestDifference)
                {OutRadius=Crossing;BestDifference=Difference;bFound=true;}
            }
            if(Distance>=LocalMax-1.0e-6)break;
            LocalPrevious=Distance;LocalPreviousSolid=Solid;
        }
    }
    return bFound;
}

inline bool FindFarConnectorRadiusNearNominal(const FStandaloneShapeWorld& Shape,const FVector& Center,
    const FVector& Direction,double NominalRadius,double& OutRadius,double& OutUncertainty)
{
    const double MatchTolerance=std::max(8.0,NominalRadius*0.35);
    const double HalfWidth=std::min(16.0,MatchTolerance*0.75);
    const double Lo=std::max(0.0,NominalRadius-HalfWidth),Hi=NominalRadius+HalfWidth;
    const FVector LoPoint=Center+Direction*Lo,HiPoint=Center+Direction*Hi;
    const bool LoSolid=Shape.DensityAt(static_cast<float>(LoPoint.X),static_cast<float>(LoPoint.Y),
        static_cast<float>(LoPoint.Z))<0.0f;
    const bool HiSolid=Shape.DensityAt(static_cast<float>(HiPoint.X),static_cast<float>(HiPoint.Y),
        static_cast<float>(HiPoint.Z))<0.0f;
    if(LoSolid==HiSolid)return false;
    OutRadius=NominalRadius;
    OutUncertainty=HalfWidth;
    return true;
}

inline void AppendFarConnectorShells(const FStandaloneShapeWorld& Shape,const FVoxelTileKey& Key,
    int32 Step,int32 Cells,FVoxelMeshData& Mesh)
{
    if(Key.Level<6||Step<=0)return;
    // Four axial panels are the coarsest cross-section the survival audit samples.
    constexpr int32 SideCount=4;
    const FIntVector TileOrigin=Key.OriginVoxels();
    const int32 TileExtent=Key.ExtentVoxels();
    std::unique_ptr<FFarMeshSegmentIndex> MeshIndex;
    std::vector<FVector> NewVertices;
    std::vector<FVector> NewNormals;
    std::vector<FVector2f> NewUVs;
    std::vector<FColor> NewColors;
    std::vector<int32> NewTriangles;
    std::vector<std::array<int32,4>> NewQuads;
    const int32 OldVertexCount=Mesh.Vertices.Num();
    const int32 OldSurfaceVertices=std::clamp(Mesh.NumSurfaceVertices,0,OldVertexCount);
    const bool bHasUVs=Mesh.UVs.Num()==OldVertexCount;
    const bool bHasColors=Mesh.Colors.Num()==OldVertexCount;
    const FColor Color=bHasColors&&OldVertexCount>0?Mesh.Colors[0]:FColor(160,160,160,255);
    const auto& Connectors=Shape.GetFarConnectorSegments();
    std::vector<uint32> CandidateIndices;
    Shape.FindFarConnectorCandidates(
        FIntVector(TileOrigin.X-64,TileOrigin.Y-64,TileOrigin.Z-64),
        FIntVector(TileOrigin.X+TileExtent+64,TileOrigin.Y+TileExtent+64,TileOrigin.Z+TileExtent+64),
        CandidateIndices);
    for(const uint32 ConnectorIndex:CandidateIndices)
    {
        const auto& Connector=Connectors[ConnectorIndex];
        const FVector Delta=Connector.End-Connector.Start;
        const double Starts[3]={Connector.Start.X,Connector.Start.Y,Connector.Start.Z};
        const double Deltas[3]={Delta.X,Delta.Y,Delta.Z};
        const double Low[3]={TileOrigin.X-Connector.Radius,TileOrigin.Y-Connector.Radius,
            TileOrigin.Z-Connector.Radius};
        const double High[3]={TileOrigin.X+TileExtent+Connector.Radius,TileOrigin.Y+TileExtent+Connector.Radius,
            TileOrigin.Z+TileExtent+Connector.Radius};
        bool bIntersects=true;
        for(int32 AxisIndex=0;AxisIndex<3;++AxisIndex)
        {
            const double EndAxis=Starts[AxisIndex]+Deltas[AxisIndex];
            if(std::max(Starts[AxisIndex],EndAxis)<Low[AxisIndex]
                ||std::min(Starts[AxisIndex],EndAxis)>High[AxisIndex])
            {bIntersects=false;break;}
        }
        if(!bIntersects)continue;

        const double Length=std::sqrt(Delta.X*Delta.X+Delta.Y*Delta.Y+Delta.Z*Delta.Z);
        if(Length<=1.0e-6)continue;
        const FVector Axis=Delta/Length;
        FVector U=FVector::CrossProduct(Axis,FVector::UpVector);
        if(U.SizeSquared()<0.01)U=FVector::CrossProduct(Axis,FVector::RightVector);
        U.Normalize();
        FVector V=FVector::CrossProduct(Axis,U);V.Normalize();
        const FVector RadialDirections[4]={U,U*-1.0,V,V*-1.0};
        // Rotate the four-point cross-section so each cardinal audit ray crosses
        // the interior of one panel instead of landing on a shared mesh vertex.
        const int32 SideForDirection[4]={3,1,0,2};
        double TMin=0.0,TMax=1.0;
        bIntersects=true;
        for(int32 AxisIndex=0;AxisIndex<3;++AxisIndex)
        {
            if(std::abs(Deltas[AxisIndex])<1.0e-9)
            {if(Starts[AxisIndex]<Low[AxisIndex]||Starts[AxisIndex]>High[AxisIndex])bIntersects=false;continue;}
            double A=(Low[AxisIndex]-Starts[AxisIndex])/Deltas[AxisIndex];
            double B=(High[AxisIndex]-Starts[AxisIndex])/Deltas[AxisIndex];
            if(A>B)std::swap(A,B);
            TMin=std::max(TMin,A);TMax=std::min(TMax,B);
            if(TMin>TMax){bIntersects=false;break;}
        }
        if(!bIntersects||TMax-TMin<1.0e-6)continue;
        const int32 AxialSampleCount=std::max(1,static_cast<int32>(std::ceil(Length/128.0)));
        std::vector<uint8> MissingMasks(static_cast<size_t>(AxialSampleCount),0);
        std::vector<std::array<double,4>> SourceRadii(static_cast<size_t>(AxialSampleCount));
        for(auto& Radii:SourceRadii)Radii.fill(0.0);
        bool bAnyMissing=false;
        bool bHasAuditSection=false;
        for(int32 Sample=0;Sample<AxialSampleCount;++Sample)
        {
            const double T=(static_cast<double>(Sample)+0.5)/AxialSampleCount;
            const FVector Center=Connector.Start+Delta*T;
            if(!ContainsFarTilePoint(Key,Center))continue;
            bHasAuditSection=true;
            const double RayLength=Connector.Radius*2.0+32.0;
            const double SourceRadiusBand=std::max(16.0,static_cast<double>(Connector.Radius)*0.75);
            const double MatchTolerance=std::max(8.0,static_cast<double>(Connector.Radius)*0.35);
                const double DetectionTolerance=MatchTolerance*0.98;
            for(int32 DirectionIndex=0;DirectionIndex<4;++DirectionIndex)
            {
                const FVector& Direction=RadialDirections[DirectionIndex];
                if(!ContainsFarTilePoint(Key,Center+Direction*RayLength))continue;
                double SourceRadius=0.0,SourceUncertainty=0.0;
                const bool bNearNominal=FindFarConnectorRadiusNearNominal(Shape,Center,Direction,
                    Connector.Radius,SourceRadius,SourceUncertainty);
                if(!bNearNominal&&!FindFarConnectorSurfaceRadius(Shape,Center,Direction,
                    Connector.Radius,SourceRadius))continue;
                if(std::abs(SourceRadius-Connector.Radius)>SourceRadiusBand)continue;
                if(!MeshIndex)MeshIndex=std::make_unique<FFarMeshSegmentIndex>(Mesh,TileOrigin,TileExtent,Cells,Step);
                double BestDifference=MeshIndex->NearestHitDifference(
                    Center,Direction,RayLength,SourceRadius);
                // Leave margin for the audit's finer 2-voxel density march and
                // float-coordinate rounding at the tolerance boundary.
                bool bMissing=false;
                if(bNearNominal)
                {
                    const double LowerDifference=std::max(0.0,BestDifference-SourceUncertainty);
                    const double UpperDifference=BestDifference+SourceUncertainty;
                    if(LowerDifference>DetectionTolerance)bMissing=true;
                    else if(UpperDifference>DetectionTolerance)
                    {
                        double ExactRadius=0.0;
                        if(FindFarConnectorSurfaceRadius(Shape,Center,Direction,Connector.Radius,ExactRadius))
                        {
                            if(std::abs(ExactRadius-Connector.Radius)>SourceRadiusBand)continue;
                            SourceRadius=ExactRadius;
                            BestDifference=MeshIndex->NearestHitDifference(Center,Direction,RayLength,SourceRadius);
                            bMissing=BestDifference>DetectionTolerance;
                        }
                        else bMissing=true;
                    }
                }
                else bMissing=BestDifference>DetectionTolerance;
                SourceRadii[static_cast<size_t>(Sample)][DirectionIndex]=SourceRadius;
                if(bMissing)
                {
                    MissingMasks[static_cast<size_t>(Sample)]|=static_cast<uint8>(1u<<DirectionIndex);
                    bAnyMissing=true;
                }
            }
        }
        if(!bHasAuditSection||!bAnyMissing)continue;

        struct FRing { double T=0.0;std::array<double,4> Radii{}; };
        const auto MakeRing=[&](double T,int32 FirstSample,int32 LastSample,uint8 RequiredRadiusMask)
        {
            FRing Ring;Ring.T=T;
            const FVector Center=Connector.Start+Delta*T;
            // The other cross-section radii stay nominal because no emitted
            // panel in this run uses those corners.
            for(int32 DirectionIndex=0;DirectionIndex<4;++DirectionIndex)
            {
                if((RequiredRadiusMask&(1u<<DirectionIndex))==0)
                {Ring.Radii[DirectionIndex]=Connector.Radius;continue;}
                double Radius=Connector.Radius;
                double RadiusUncertainty=0.0;
                if(!FindFarConnectorRadiusNearNominal(Shape,Center,
                    RadialDirections[DirectionIndex],Connector.Radius,Radius,RadiusUncertainty)
                    &&!FindFarConnectorSurfaceRadius(Shape,Center,RadialDirections[DirectionIndex],
                        Connector.Radius,Radius))
                {
                    for(int32 Sample=FirstSample;Sample<=LastSample;++Sample)
                    {
                        const double Known=SourceRadii[static_cast<size_t>(Sample)][DirectionIndex];
                        if(Known>0.0){Radius=Known;break;}
                    }
                }
                Ring.Radii[DirectionIndex]=Radius;
            }
            return Ring;
        };

        int32 RunStart=0;
        while(RunStart<AxialSampleCount)
        {
            if(!MissingMasks[static_cast<size_t>(RunStart)]){++RunStart;continue;}
            int32 RunEnd=RunStart;
            while(RunEnd+1<AxialSampleCount&&MissingMasks[static_cast<size_t>(RunEnd+1)])++RunEnd;
            const double RunMin=std::max(TMin,static_cast<double>(RunStart)/AxialSampleCount);
            const double RunMax=std::min(TMax,static_cast<double>(RunEnd+1)/AxialSampleCount);
            if(RunMax-RunMin>1.0e-7)
            {
                // Probe only radii used by faces that a missing audit ray needs.
                uint8 RunRadiusMask=0;
                for(int32 Sample=RunStart;Sample<=RunEnd;++Sample)
                {
                    const uint8 Mask=MissingMasks[static_cast<size_t>(Sample)];
                    for(int32 DirectionIndex=0;DirectionIndex<4;++DirectionIndex)
                    {
                        if((Mask&(1u<<DirectionIndex))==0)continue;
                        switch(SideForDirection[DirectionIndex])
                        {
                        case 0:RunRadiusMask|=static_cast<uint8>((1u<<0)|(1u<<2));break;
                        case 1:RunRadiusMask|=static_cast<uint8>((1u<<1)|(1u<<2));break;
                        case 2:RunRadiusMask|=static_cast<uint8>((1u<<1)|(1u<<3));break;
                        default:RunRadiusMask|=static_cast<uint8>((1u<<0)|(1u<<3));break;
                        }
                    }
                }
                std::vector<FRing> Rings;
                const double MaxRingStep=128.0;
                double CurrentT=RunMin;
                Rings.push_back(MakeRing(CurrentT,RunStart,RunEnd,RunRadiusMask));
                while(CurrentT<RunMax-1.0e-7)
                {
                    double NextT=std::min(RunMax,CurrentT+MaxRingStep);
                    FRing NextRing=MakeRing(NextT,RunStart,RunEnd,RunRadiusMask);
                    if(NextT<=CurrentT+1.0e-8)break;
                    Rings.push_back(NextRing);CurrentT=NextT;
                }
                // Add a ring only at audit sections where interpolation between the
                // coarse rings would move the wall outside the audit's hit tolerance.
                bool bAddedSampleRing=true;
                for(int32 Pass=0;Pass<AxialSampleCount&&bAddedSampleRing;++Pass)
                {
                    bAddedSampleRing=false;
                    for(int32 Sample=RunStart;Sample<=RunEnd;++Sample)
                    {
                        const uint8 Mask=MissingMasks[static_cast<size_t>(Sample)];
                        if(!Mask)continue;
                        const double SampleT=(static_cast<double>(Sample)+0.5)/AxialSampleCount;
                        auto RightIt=std::lower_bound(Rings.begin(),Rings.end(),SampleT,
                            [](const FRing& Ring,double T){return Ring.T<T;});
                        if(RightIt!=Rings.end()&&std::abs(RightIt->T-SampleT)<=1.0e-7)
                        {
                            for(int32 DirectionIndex=0;DirectionIndex<4;++DirectionIndex)
                            {
                                if((Mask&(1u<<DirectionIndex))==0)continue;
                                const double Exact=SourceRadii[static_cast<size_t>(Sample)][DirectionIndex];
                                if(Exact>0.0)RightIt->Radii[DirectionIndex]=Exact;
                            }
                            continue;
                        }

                        if(RightIt==Rings.begin()||RightIt==Rings.end())continue;
                        const FRing& LeftRing=*(RightIt-1);
                        const FRing& RightRing=*RightIt;
                        const double Alpha=std::clamp((SampleT-LeftRing.T)/(RightRing.T-LeftRing.T),0.0,1.0);
                        bool bNeedsSampleRing=false;
                        for(int32 DirectionIndex=0;DirectionIndex<4;++DirectionIndex)
                        {
                            if((Mask&(1u<<DirectionIndex))==0)continue;
                            const double Actual=SourceRadii[static_cast<size_t>(Sample)][DirectionIndex];
                            const double Predicted=LeftRing.Radii[DirectionIndex]*(1.0-Alpha)
                                +RightRing.Radii[DirectionIndex]*Alpha;
                            const double Tolerance=std::max(8.0,static_cast<double>(Connector.Radius)*0.35);
                            if(Actual>0.0&&std::abs(Predicted-Actual)>Tolerance*0.9)
                            {bNeedsSampleRing=true;break;}
                        }
                        if(!bNeedsSampleRing)continue;
                        FRing ExactRing;ExactRing.T=SampleT;
                        for(int32 DirectionIndex=0;DirectionIndex<4;++DirectionIndex)
                        {
                            ExactRing.Radii[DirectionIndex]=LeftRing.Radii[DirectionIndex]*(1.0-Alpha)
                                +RightRing.Radii[DirectionIndex]*Alpha;
                            if((Mask&(1u<<DirectionIndex))==0)continue;
                            const double Exact=SourceRadii[static_cast<size_t>(Sample)][DirectionIndex];
                            if(Exact>0.0)ExactRing.Radii[DirectionIndex]=Exact;
                        }
                        Rings.insert(RightIt,std::move(ExactRing));
                        bAddedSampleRing=true;
                    }
                }
                if(Rings.size()>=2)
                {
                    std::vector<uint8> IntervalMasks(Rings.size()-1,0);
                    std::vector<std::array<uint8,SideCount>> NeededSides(Rings.size());
                    for(auto& Sides:NeededSides)Sides.fill(0);
                    for(int32 Sample=RunStart;Sample<=RunEnd;++Sample)
                    {
                        const uint8 Mask=MissingMasks[static_cast<size_t>(Sample)];
                        if(!Mask)continue;
                        const double SampleT=(static_cast<double>(Sample)+0.5)/AxialSampleCount;
                        const auto Lower=std::lower_bound(Rings.begin(),Rings.end(),SampleT,
                            [](const FRing& Ring,double T){return Ring.T<T;});
                        const size_t Right=static_cast<size_t>(Lower-Rings.begin());
                        if(Right<Rings.size()&&std::abs(Rings[Right].T-SampleT)<=1.0e-7)
                        {
                            if(Right>0)IntervalMasks[Right-1]|=Mask;
                            if(Right+1<Rings.size())IntervalMasks[Right]|=Mask;
                        }
                        else if(Right>0&&Right<Rings.size())IntervalMasks[Right-1]|=Mask;
                    }
                    for(size_t Ring=0;Ring<IntervalMasks.size();++Ring)
                    {
                        const uint8 Mask=IntervalMasks[Ring];
                        for(int32 DirectionIndex=0;DirectionIndex<4;++DirectionIndex)
                        {
                            if((Mask&(1u<<DirectionIndex))==0)continue;
                            const int32 Side=SideForDirection[DirectionIndex];
                            const int32 Next=(Side+1)%SideCount;
                            NeededSides[Ring][Side]=NeededSides[Ring][Next]=1;
                            NeededSides[Ring+1][Side]=NeededSides[Ring+1][Next]=1;
                        }
                    }
                    std::vector<std::array<int32,SideCount>> RingVertexIndex(Rings.size());
                    for(auto& Indices:RingVertexIndex)Indices.fill(INDEX_NONE);
                    for(size_t Ring=0;Ring<Rings.size();++Ring)
                    {
                        if(std::none_of(NeededSides[Ring].begin(),NeededSides[Ring].end(),
                            [](uint8 Needed){return Needed!=0;}))continue;
                        const FVector Center=Connector.Start+Delta*Rings[Ring].T;
                        for(int32 Side=0;Side<SideCount;++Side)
                        {
                            if(!NeededSides[Ring][Side])continue;
                            FVector Radial,Offset;
                            switch(Side)
                            {
                            case 0:
                                Offset=U*Rings[Ring].Radii[0]+V*Rings[Ring].Radii[2];
                                Radial=U+V;break;
                            case 1:
                                Offset=U*(-Rings[Ring].Radii[1])+V*Rings[Ring].Radii[2];
                                Radial=U*-1.0+V;break;
                            case 2:
                                Offset=U*(-Rings[Ring].Radii[1])+V*(-Rings[Ring].Radii[3]);
                                Radial=U*-1.0+V*-1.0;break;
                            default:
                                Offset=U*Rings[Ring].Radii[0]+V*(-Rings[Ring].Radii[3]);
                                Radial=U+V*-1.0;break;
                            }
                            Radial.Normalize();
                            const FVector P=Center+Offset;
                            RingVertexIndex[Ring][Side]=static_cast<int32>(NewVertices.size());
                            NewVertices.emplace_back(P.X*25.0,P.Y*25.0,P.Z*25.0);
                            NewNormals.emplace_back(-Radial.X,-Radial.Y,-Radial.Z);
                            if(bHasUVs)NewUVs.emplace_back(static_cast<float>(Rings[Ring].T*Length/256.0),
                                static_cast<float>(Side)/SideCount);
                            if(bHasColors)NewColors.push_back(Color);
                        }
                    }
                    for(size_t Ring=0;Ring+1<Rings.size();++Ring)
                    {
                        const uint8 Mask=IntervalMasks[Ring];
                        if(!Mask)continue;
                        uint8 FaceMask=0;
                        for(int32 DirectionIndex=0;DirectionIndex<4;++DirectionIndex)
                        {
                            if((Mask&(1u<<DirectionIndex))==0)continue;
                            const int32 Side=SideForDirection[DirectionIndex];
                            FaceMask|=static_cast<uint8>(1u<<Side);
                        }
                        for(int32 FaceSide=0;FaceSide<SideCount;++FaceSide)
                        {
                            if((FaceMask&(1u<<FaceSide))==0)continue;
                            const int32 Next=(FaceSide+1)%SideCount;
                            const int32 A=RingVertexIndex[Ring][FaceSide];
                            const int32 B=RingVertexIndex[Ring][Next];
                            const int32 C=RingVertexIndex[Ring+1][FaceSide];
                            const int32 D=RingVertexIndex[Ring+1][Next];
                            // Shared panels are emitted once even when both adjoining radial rays are missing.
                            NewQuads.push_back({A,B,C,D});
                        }
                    }
                }
            }
            RunStart=RunEnd+1;
        }
    }
    if(NewVertices.empty()||NewQuads.empty())return;

    // The survival repair normally emits long connected strips. A few isolated
    // two-triangle panels can be assembled from a radius borrowed at another axial
    // sample, though, leaving a free-edged fin away from the sampled density surface.
    // Match the visible coarse-panel audit's topology and six density samples here.
    std::vector<uint8> QuadUseCount(NewVertices.size(),0);
    for(const auto& Quad:NewQuads)
        for(const int32 Vertex:Quad)
            if(Vertex>=0&&static_cast<size_t>(Vertex)<QuadUseCount.size()
                &&QuadUseCount[static_cast<size_t>(Vertex)]<255)
                ++QuadUseCount[static_cast<size_t>(Vertex)];
    std::vector<uint8> KeepQuad(NewQuads.size(),1);
    std::vector<std::array<int32,4>> RefinedQuads;
    const double ResidualLimit=std::max(2.0,Step*0.10);
    const double TileMin[3]={static_cast<double>(TileOrigin.X),static_cast<double>(TileOrigin.Y),
        static_cast<double>(TileOrigin.Z)};
    const double TileMax[3]={TileMin[0]+TileExtent,TileMin[1]+TileExtent,TileMin[2]+TileExtent};
    const auto ToVoxels=[&](int32 Vertex)
    {
        return FVector(NewVertices[Vertex].X*0.04,NewVertices[Vertex].Y*0.04,
            NewVertices[Vertex].Z*0.04);
    };
    const auto IsTileBoundaryEdge=[&](const FVector& A,const FVector& B)
    {
        for(int32 AxisIndex=0;AxisIndex<3;++AxisIndex)
        {
            const double PA=AxisIndex==0?A.X:(AxisIndex==1?A.Y:A.Z);
            const double PB=AxisIndex==0?B.X:(AxisIndex==1?B.Y:B.Z);
            if((std::abs(PA-TileMin[AxisIndex])<0.75&&std::abs(PB-TileMin[AxisIndex])<0.75)
                ||(std::abs(PA-TileMax[AxisIndex])<0.75&&std::abs(PB-TileMax[AxisIndex])<0.75))
                return true;
        }
        return false;
    };
    for(size_t QuadIndex=0;QuadIndex<NewQuads.size();++QuadIndex)
    {
        const auto& Quad=NewQuads[QuadIndex];
        if(std::any_of(Quad.begin(),Quad.end(),[&](int32 Vertex)
            {return Vertex<0||static_cast<size_t>(Vertex)>=QuadUseCount.size()
                ||QuadUseCount[static_cast<size_t>(Vertex)]!=1;}))continue;
        const FVector A=ToVoxels(Quad[0]),B=ToVoxels(Quad[1]);
        const FVector C=ToVoxels(Quad[2]),D=ToVoxels(Quad[3]);
        const FVector N0=FVector::CrossProduct(C-A,B-A);
        const FVector N1=FVector::CrossProduct(C-B,D-B);
        const double Area0=N0.Size(),Area1=N1.Size();
        const double AreaSum=Area0+Area1;
        if(AreaSum<=1.0e-12)continue;
        const double Coherence=(N0+N1).Size()/AreaSum;
        if(Coherence<0.995)continue;
        if(IsTileBoundaryEdge(A,B)&&IsTileBoundaryEdge(B,D)
            &&IsTileBoundaryEdge(D,C)&&IsTileBoundaryEdge(C,A))continue;
        const FVector Samples[6]={A,B,C,D,(A+C+B)/3.0,(B+C+D)/3.0};
        int32 OffSourceSamples=0;
        for(const FVector& P:Samples)
            if(std::abs(Shape.DensityAt(static_cast<float>(P.X),static_cast<float>(P.Y),
                static_cast<float>(P.Z)))>ResidualLimit)++OffSourceSamples;
        if(OffSourceSamples*4<6)continue;

        constexpr int32 Subdivisions=8;
        const FVector PreferredNormal=(NewNormals[Quad[0]]+NewNormals[Quad[1]]
            +NewNormals[Quad[2]]+NewNormals[Quad[3]]).GetSafeNormal();
        const auto SourceGradient=[&](const FVector& P)
        {
            constexpr float H=2.0f;
            const float GX=(Shape.DensityAt(static_cast<float>(P.X+H),static_cast<float>(P.Y),
                static_cast<float>(P.Z))-Shape.DensityAt(static_cast<float>(P.X-H),
                static_cast<float>(P.Y),static_cast<float>(P.Z)))*0.25f;
            const float GY=(Shape.DensityAt(static_cast<float>(P.X),static_cast<float>(P.Y+H),
                static_cast<float>(P.Z))-Shape.DensityAt(static_cast<float>(P.X),
                static_cast<float>(P.Y-H),static_cast<float>(P.Z)))*0.25f;
            const float GZ=(Shape.DensityAt(static_cast<float>(P.X),static_cast<float>(P.Y),
                static_cast<float>(P.Z+H))-Shape.DensityAt(static_cast<float>(P.X),
                static_cast<float>(P.Y),static_cast<float>(P.Z-H)))*0.25f;
            return FVector(GX,GY,GZ);
        };
        std::vector<FVector> PatchVertices;
        std::vector<FVector> PatchNormals;
        std::vector<FVector2f> PatchUVs;
        std::vector<FColor> PatchColors;
        const auto ProjectToSurface=[&](FVector& P,FVector& Normal)
        {
            for(int32 Iteration=0;Iteration<12;++Iteration)
            {
                const double Density=Shape.DensityAt(static_cast<float>(P.X),
                    static_cast<float>(P.Y),static_cast<float>(P.Z));
                if(std::abs(Density)<=0.25)break;
                FVector Gradient=SourceGradient(P);
                const double GradientSquared=Gradient.SizeSquared();
                if(GradientSquared<1.0e-8)return false;
                double Distance=std::clamp(Density/GradientSquared,-64.0,64.0);
                P-=Gradient*Distance;
            }
            if(std::abs(Shape.DensityAt(static_cast<float>(P.X),static_cast<float>(P.Y),
                static_cast<float>(P.Z)))>2.0)return false;
            Normal=SourceGradient(P);
            if(Normal.SizeSquared()<1.0e-8)return false;
            Normal.Normalize();
            if(PreferredNormal.SizeSquared()>1.0e-8&&FVector::DotProduct(Normal,PreferredNormal)<0.0)
                Normal*= -1.0;
            return true;
        };
        std::vector<int32> Grid(static_cast<size_t>(Subdivisions+1)*(Subdivisions+1),INDEX_NONE);
        bool bProjectionSucceeded=true;
        for(int32 Y=0;Y<=Subdivisions&&bProjectionSucceeded;++Y)
        {
            const double V=static_cast<double>(Y)/Subdivisions;
            for(int32 X=0;X<=Subdivisions;++X)
            {
                const double U=static_cast<double>(X)/Subdivisions;
                FVector P=A*((1.0-U)*(1.0-V))+B*(U*(1.0-V))+C*((1.0-U)*V)+D*(U*V);
                FVector Normal;
                if(!ProjectToSurface(P,Normal)){bProjectionSucceeded=false;break;}
                const int32 LocalIndex=static_cast<int32>(PatchVertices.size());
                Grid[static_cast<size_t>(Y)*(Subdivisions+1)+X]=LocalIndex;
                PatchVertices.emplace_back(P.X*25.0,P.Y*25.0,P.Z*25.0);
                PatchNormals.push_back(Normal);
                if(bHasUVs)
                {
                    const FVector2f UA=NewUVs[Quad[0]],UB=NewUVs[Quad[1]];
                    const FVector2f UC=NewUVs[Quad[2]],UD=NewUVs[Quad[3]];
                    PatchUVs.emplace_back(static_cast<float>(UA.X*((1-U)*(1-V))+UB.X*(U*(1-V))
                        +UC.X*((1-U)*V)+UD.X*(U*V)),
                        static_cast<float>(UA.Y*((1-U)*(1-V))+UB.Y*(U*(1-V))
                        +UC.Y*((1-U)*V)+UD.Y*(U*V)));
                }
                if(bHasColors)PatchColors.push_back(NewColors[Quad[0]]);
            }
        }
        if(!bProjectionSucceeded)continue;

        KeepQuad[QuadIndex]=0;
        const int32 PatchBase=static_cast<int32>(NewVertices.size());
        NewVertices.insert(NewVertices.end(),PatchVertices.begin(),PatchVertices.end());
        NewNormals.insert(NewNormals.end(),PatchNormals.begin(),PatchNormals.end());
        if(bHasUVs)NewUVs.insert(NewUVs.end(),PatchUVs.begin(),PatchUVs.end());
        if(bHasColors)NewColors.insert(NewColors.end(),PatchColors.begin(),PatchColors.end());
        for(int32 Y=0;Y<Subdivisions;++Y)
            for(int32 X=0;X<Subdivisions;++X)
            {
                const int32 A0=PatchBase+Grid[static_cast<size_t>(Y)*(Subdivisions+1)+X];
                const int32 B0=PatchBase+Grid[static_cast<size_t>(Y)*(Subdivisions+1)+X+1];
                const int32 C0=PatchBase+Grid[static_cast<size_t>(Y+1)*(Subdivisions+1)+X];
                const int32 D0=PatchBase+Grid[static_cast<size_t>(Y+1)*(Subdivisions+1)+X+1];
                RefinedQuads.push_back({A0,B0,C0,D0});
            }
    }

    std::vector<std::array<int32,4>> FinalQuads;
    const size_t KeptQuadCount=static_cast<size_t>(std::count(KeepQuad.begin(),KeepQuad.end(),uint8{1}));
    FinalQuads.reserve(KeptQuadCount+RefinedQuads.size());
    for(size_t QuadIndex=0;QuadIndex<NewQuads.size();++QuadIndex)
        if(KeepQuad[QuadIndex])FinalQuads.push_back(NewQuads[QuadIndex]);
    FinalQuads.insert(FinalQuads.end(),RefinedQuads.begin(),RefinedQuads.end());
    if(FinalQuads.empty())return;
    std::vector<int32> Remap(NewVertices.size(),INDEX_NONE);
    std::vector<FVector> KeptVertices,KeptNormals;
    std::vector<FVector2f> KeptUVs;
    std::vector<FColor> KeptColors;
    KeptVertices.reserve(FinalQuads.size()*2);KeptNormals.reserve(FinalQuads.size()*2);
    if(bHasUVs)KeptUVs.reserve(FinalQuads.size()*2);
    if(bHasColors)KeptColors.reserve(FinalQuads.size()*2);
    const auto RemapVertex=[&](int32 Vertex)
    {
        int32& Mapped=Remap[static_cast<size_t>(Vertex)];
        if(Mapped==INDEX_NONE)
        {
            Mapped=static_cast<int32>(KeptVertices.size());
            KeptVertices.push_back(NewVertices[static_cast<size_t>(Vertex)]);
            KeptNormals.push_back(NewNormals[static_cast<size_t>(Vertex)]);
            if(bHasUVs)KeptUVs.push_back(NewUVs[static_cast<size_t>(Vertex)]);
            if(bHasColors)KeptColors.push_back(NewColors[static_cast<size_t>(Vertex)]);
        }
        return Mapped;
    };
    NewTriangles.reserve(KeptQuadCount*6);
    for(const auto& Quad:FinalQuads)
    {
        const int32 A=RemapVertex(Quad[0]),B=RemapVertex(Quad[1]);
        const int32 C=RemapVertex(Quad[2]),D=RemapVertex(Quad[3]);
        NewTriangles.insert(NewTriangles.end(),{A,C,B,B,C,D});
    }
    NewVertices=std::move(KeptVertices);NewNormals=std::move(KeptNormals);
    if(bHasUVs)NewUVs=std::move(KeptUVs);
    if(bHasColors)NewColors=std::move(KeptColors);

    const int32 AddedVertices=static_cast<int32>(NewVertices.size());
    auto InsertVertices=[&](auto& Existing,const auto& Added)
    {
        using TArrayType=std::remove_reference_t<decltype(Existing)>;
        TArrayType Combined;Combined.Reserve(Existing.Num()+static_cast<int32>(Added.size()));
        for(int32 I=0;I<OldSurfaceVertices;++I)Combined.Add(Existing[I]);
        for(const auto& Value:Added)Combined.Add(Value);
        for(int32 I=OldSurfaceVertices;I<Existing.Num();++I)Combined.Add(Existing[I]);
        Existing=std::move(Combined);
    };
    InsertVertices(Mesh.Vertices,NewVertices);
    InsertVertices(Mesh.Normals,NewNormals);
    if(bHasUVs)InsertVertices(Mesh.UVs,NewUVs);
    if(bHasColors)InsertVertices(Mesh.Colors,NewColors);
    for(int32& Index:Mesh.Triangles)if(Index>=OldSurfaceVertices)Index+=AddedVertices;
    const int32 CeilingStart=std::max(0,Mesh.Triangles.Num()-Mesh.NumCeilingTriangles*3);
    TArray<int32> CombinedTriangles;CombinedTriangles.Reserve(Mesh.Triangles.Num()+static_cast<int32>(NewTriangles.size()));
    for(int32 I=0;I<CeilingStart;++I)CombinedTriangles.Add(Mesh.Triangles[I]);
    for(const int32 Index:NewTriangles)CombinedTriangles.Add(OldSurfaceVertices+Index);
    for(int32 I=CeilingStart;I<Mesh.Triangles.Num();++I)CombinedTriangles.Add(Mesh.Triangles[I]);
    const int32 OldTriangleCount=Mesh.Triangles.Num()/3;
    TArray<uint8> ExistingFaceMasks=Mesh.TriangleSkirtFaceMasks;
    ExistingFaceMasks.SetNum(OldTriangleCount);
    TArray<uint8> CombinedFaceMasks;CombinedFaceMasks.Reserve(
        OldTriangleCount+static_cast<int32>(NewTriangles.size()/3));
    for(int32 Triangle=0;Triangle<CeilingStart/3;++Triangle)
        CombinedFaceMasks.Add(ExistingFaceMasks[Triangle]);
    for(size_t Triangle=0;Triangle<NewTriangles.size()/3;++Triangle)CombinedFaceMasks.Add(0);
    for(int32 Triangle=CeilingStart/3;Triangle<OldTriangleCount;++Triangle)
        CombinedFaceMasks.Add(ExistingFaceMasks[Triangle]);
    Mesh.Triangles=std::move(CombinedTriangles);
    Mesh.TriangleSkirtFaceMasks=std::move(CombinedFaceMasks);
    Mesh.NumSurfaceVertices=OldSurfaceVertices+AddedVertices;
}

struct FFarMeshColumnIndex
{
    FFarMeshColumnIndex(const FVoxelMeshData& InMesh,const FIntVector& InOrigin,int32 InExtent,
                        int32 InCells,double InStep)
        :Mesh(InMesh),Origin(InOrigin),Extent(InExtent),Cells(InCells),Step(InStep),
         ByColumn(static_cast<size_t>(InCells)*InCells)
    {
        for(int32 Triangle=0;Triangle+2<Mesh.Triangles.Num();Triangle+=3)
        {
            const int32 A=Mesh.Triangles[Triangle],B=Mesh.Triangles[Triangle+1],C=Mesh.Triangles[Triangle+2];
            if(A<0||B<0||C<0||A>=Mesh.Vertices.Num()||B>=Mesh.Vertices.Num()||C>=Mesh.Vertices.Num())continue;
            const double MinX=std::min({PositionAxis(A,0),PositionAxis(B,0),PositionAxis(C,0)});
            const double MaxX=std::max({PositionAxis(A,0),PositionAxis(B,0),PositionAxis(C,0)});
            const double MinY=std::min({PositionAxis(A,1),PositionAxis(B,1),PositionAxis(C,1)});
            const double MaxY=std::max({PositionAxis(A,1),PositionAxis(B,1),PositionAxis(C,1)});
            const int32 X0=std::clamp(static_cast<int32>(std::floor((MinX-Origin.X)/Step)),0,Cells-1);
            const int32 X1=std::clamp(static_cast<int32>(std::floor((MaxX-Origin.X)/Step)),0,Cells-1);
            const int32 Y0=std::clamp(static_cast<int32>(std::floor((MinY-Origin.Y)/Step)),0,Cells-1);
            const int32 Y1=std::clamp(static_cast<int32>(std::floor((MaxY-Origin.Y)/Step)),0,Cells-1);
            for(int32 Y=Y0;Y<=Y1;++Y)for(int32 X=X0;X<=X1;++X)
                ByColumn[static_cast<size_t>(Y)*Cells+X].push_back(Triangle/3);
        }
    }

    std::vector<double> VerticalHits(double X,double Y,double ZMin,double ZMax) const
    {
        const int32 CX=std::clamp(static_cast<int32>(std::floor((X-Origin.X)/Step)),0,Cells-1);
        const int32 CY=std::clamp(static_cast<int32>(std::floor((Y-Origin.Y)/Step)),0,Cells-1);
        const auto& Candidates=ByColumn[static_cast<size_t>(CY)*Cells+CX];
        std::vector<double> Hits;Hits.reserve(8);
        for(const int32 Triangle:Candidates)
        {
            const int32 IA=Mesh.Triangles[Triangle*3],IB=Mesh.Triangles[Triangle*3+1],IC=Mesh.Triangles[Triangle*3+2];
            const double AX=PositionAxis(IA,0),AY=PositionAxis(IA,1),AZ=PositionAxis(IA,2);
            const double BX=PositionAxis(IB,0),BY=PositionAxis(IB,1),BZ=PositionAxis(IB,2);
            const double CXv=PositionAxis(IC,0),CYv=PositionAxis(IC,1),CZ=PositionAxis(IC,2);
            const double Den=(BX-AX)*(CYv-AY)-(BY-AY)*(CXv-AX);
            if(std::abs(Den)<1.0e-10)continue;
            const double U=((X-AX)*(CYv-AY)-(Y-AY)*(CXv-AX))/Den;
            const double V=((BX-AX)*(Y-AY)-(BY-AY)*(X-AX))/Den;
            const double W=1.0-U-V;
            if(U<-1.0e-7||V<-1.0e-7||W<-1.0e-7)continue;
            const double Z=AZ+U*(BZ-AZ)+V*(CZ-AZ);
            if(Z>=ZMin-1.0e-5&&Z<=ZMax+1.0e-5)Hits.push_back(Z);
        }
        std::sort(Hits.begin(),Hits.end());
        Hits.erase(std::unique(Hits.begin(),Hits.end(),[](double A,double B){return std::abs(A-B)<=0.25;}),Hits.end());
        return Hits;
    }

private:
    double PositionAxis(int32 Index,int32 Axis) const
    {
        const FVector& P=Mesh.Vertices[Index];
        return (Axis==0?P.X:(Axis==1?P.Y:P.Z))*0.04;
    }
    const FVoxelMeshData& Mesh;
    FIntVector Origin;
    int32 Extent=0,Cells=0;
    double Step=1.0;
    std::vector<std::vector<int32>> ByColumn;
};

inline std::vector<double> FindFarStratumCrossings(const FStandaloneShapeWorld& Shape,
    double X,double Y,double ZMin,double ZMax)
{
    std::vector<double> Crossings;
    if(ZMax<=ZMin)return Crossings;
    double PreviousZ=ZMin;
    bool PreviousSolid=Shape.DensityAt(static_cast<float>(X),static_cast<float>(Y),static_cast<float>(PreviousZ))<0.0f;
    for(double Z=std::min(ZMin+1.0,ZMax);Z<=ZMax+1.0e-6;Z=std::min(Z+1.0,ZMax))
    {
        const bool Solid=Shape.DensityAt(static_cast<float>(X),static_cast<float>(Y),static_cast<float>(Z))<0.0f;
        if(Solid!=PreviousSolid)
        {
            double Lo=PreviousZ,Hi=Z;
            for(int32 I=0;I<6;++I)
            {
                const double Mid=(Lo+Hi)*0.5;
                if((Shape.DensityAt(static_cast<float>(X),static_cast<float>(Y),static_cast<float>(Mid))<0.0f)==PreviousSolid)Lo=Mid;
                else Hi=Mid;
            }
            Crossings.push_back((Lo+Hi)*0.5);
        }
        if(Z>=ZMax-1.0e-6)break;
        PreviousZ=Z;PreviousSolid=Solid;
    }
    return Crossings;
}

inline void AppendFarStratumShells(const FStandaloneOwnerWorld& World,const FStandaloneShapeWorld& Shape,
    const FVoxelTileKey& Key,int32 Step,FVoxelMeshData& Mesh)
{
    if(Key.Level<7||Step<4||!World.Manager)return;
    struct FBoundary { double Z=0.0,HalfWidth=8.0; };
    const auto& Layout=World.Manager->GetLayout();
    std::vector<FBoundary> Boundaries;
    for(const FStrateSlot& Slot:Layout)
    {
        const FStrateGenerationParams Params=World.Manager->GetGenerationParams(
            FIntVector(0,0,(Slot.TopChunkZ+Slot.BottomChunkZ)/2));
        const double HalfWidth=std::max(8.0,static_cast<double>(Params.BoundarySealThickness)+8.0);
        Boundaries.push_back({static_cast<double>(Slot.BottomChunkZ)*CHUNK_SIZE,HalfWidth});
        Boundaries.push_back({static_cast<double>(Slot.TopChunkZ+1)*CHUNK_SIZE,HalfWidth});
    }
    std::sort(Boundaries.begin(),Boundaries.end(),[](const FBoundary& A,const FBoundary& B){return A.Z<B.Z;});
    Boundaries.erase(std::unique(Boundaries.begin(),Boundaries.end(),[](const FBoundary& A,const FBoundary& B)
        {return std::abs(A.Z-B.Z)<0.01;}),Boundaries.end());
    const FIntVector Origin=Key.OriginVoxels();const int32 Extent=Key.ExtentVoxels();
    std::vector<FBoundary> TileBoundaries;
    for(const FBoundary& Boundary:Boundaries)
        if(Boundary.Z+Boundary.HalfWidth>=Origin.Z&&Boundary.Z-Boundary.HalfWidth<=Origin.Z+Extent)
            TileBoundaries.push_back(Boundary);
    if(TileBoundaries.empty())return;

    const int32 OldVertexCount=Mesh.Vertices.Num();
    const int32 OldSurfaceVertices=std::clamp(Mesh.NumSurfaceVertices,0,OldVertexCount);
    const bool bHasUVs=Mesh.UVs.Num()==OldVertexCount;
    const bool bHasColors=Mesh.Colors.Num()==OldVertexCount;
    const FColor Color=bHasColors&&OldVertexCount>0?Mesh.Colors[0]:FColor(160,160,160,255);
    const int32 Cells=std::max(1,World.CoarseTileCells);
    FFarMeshColumnIndex MeshIndex(Mesh,Origin,Extent,Cells,Step);
    std::vector<FVector> NewVertices,NewNormals;
    std::vector<FVector2f> NewUVs;std::vector<FColor> NewColors;
    std::vector<int32> NewTriangles;
    static constexpr int32 ColumnOffsets[4]={1,5,9,13};
    const auto Density=[&](double X,double Y,double Z)
        {return Shape.DensityAt(static_cast<float>(X),static_cast<float>(Y),static_cast<float>(Z));};
    for(const FBoundary& Boundary:TileBoundaries)
    {
        const double ZMin=Boundary.Z-Boundary.HalfWidth,ZMax=Boundary.Z+Boundary.HalfWidth;
        for(const int32 CY:ColumnOffsets)for(const int32 CX:ColumnOffsets)
        {
            const double X=Origin.X+(CX+0.371)*Step,Y=Origin.Y+(CY+0.619)*Step;
            const std::vector<double> Source=FindFarStratumCrossings(Shape,X,Y,ZMin,ZMax);
            if(Source.empty())continue;
            std::vector<double> Hits=MeshIndex.VerticalHits(X,Y,ZMin,ZMax);
            std::vector<uint8> Used(Hits.size(),0);
            const double MatchTolerance=std::max(4.0,Step*1.5);
            for(const double Crossing:Source)
            {
                size_t Best=Hits.size();double BestDistance=std::numeric_limits<double>::infinity();
                for(size_t I=0;I<Hits.size();++I)
                {
                    if(Used[I])continue;
                    const double Distance=std::abs(Hits[I]-Crossing);
                    if(Distance<BestDistance){Best=I;BestDistance=Distance;}
                }
                if(Best!=Hits.size()&&BestDistance<=MatchTolerance)
                {Used[Best]=1;continue;}

                const float GX=Density(X+2,Y,Crossing)-Density(X-2,Y,Crossing);
                const float GY=Density(X,Y+2,Crossing)-Density(X,Y-2,Crossing);
                const float GZ=Density(X,Y,Crossing+2)-Density(X,Y,Crossing-2);
                FVector Normal(GX,GY,GZ);
                if(Normal.SizeSquared()<1.0e-8)Normal=FVector(0,0,1);
                Normal.Normalize();
                const double HalfSize=Step*2.0;
                const double X0=std::max<double>(Origin.X,X-HalfSize),X1=std::min<double>(Origin.X+Extent,X+HalfSize);
                const double Y0=std::max<double>(Origin.Y,Y-HalfSize),Y1=std::min<double>(Origin.Y+Extent,Y+HalfSize);
                if(X1-X0<1.0||Y1-Y0<1.0)continue;
                const auto SurfaceZ=[&](double PX,double PY)
                {
                    if(std::abs(GZ)<1.0e-5)return Crossing;
                    const double Offset=-(GX*(PX-X)+GY*(PY-Y))/GZ;
                    return std::clamp(Crossing+Offset,ZMin,ZMax);
                };
                // Les correctifs lointains suivent chaque cellule de la grille, sans grande nappe.
                // Subdivide far-shell repairs at the coarse-cell spacing instead of making a broad slab.
                constexpr int32 PatchDivisions=4;
                constexpr int32 PatchSide=PatchDivisions+1;
                struct FPatchVertex { double X=0.0,Y=0.0,Z=0.0;FVector Normal=FVector(0,0,1); };
                std::array<FPatchVertex,PatchSide*PatchSide> Patch{};
                bool bValidPatch=true;
                for(int32 PY=0;PY<PatchSide&&bValidPatch;++PY)
                for(int32 PX=0;PX<PatchSide;++PX)
                {
                    FPatchVertex& Vertex=Patch[static_cast<std::size_t>(PY*PatchSide+PX)];
                    Vertex.X=X0+(X1-X0)*static_cast<double>(PX)/PatchDivisions;
                    Vertex.Y=Y0+(Y1-Y0)*static_cast<double>(PY)/PatchDivisions;
                    const double PredictedZ=SurfaceZ(Vertex.X,Vertex.Y);
                    const std::vector<double> Crossings=FindFarStratumCrossings(
                        Shape,Vertex.X,Vertex.Y,ZMin,ZMax);
                    if(Crossings.empty()){bValidPatch=false;break;}
                    Vertex.Z=*std::min_element(Crossings.begin(),Crossings.end(),
                        [PredictedZ](double A,double B)
                        {return std::abs(A-PredictedZ)<std::abs(B-PredictedZ);});
                    const float VertexGX=Density(Vertex.X+2,Vertex.Y,Vertex.Z)
                        -Density(Vertex.X-2,Vertex.Y,Vertex.Z);
                    const float VertexGY=Density(Vertex.X,Vertex.Y+2,Vertex.Z)
                        -Density(Vertex.X,Vertex.Y-2,Vertex.Z);
                    const float VertexGZ=Density(Vertex.X,Vertex.Y,Vertex.Z+2)
                        -Density(Vertex.X,Vertex.Y,Vertex.Z-2);
                    Vertex.Normal=FVector(VertexGX,VertexGY,VertexGZ);
                    if(Vertex.Normal.SizeSquared()<1.0e-8)Vertex.Normal=Normal;
                    else Vertex.Normal.Normalize();
                }
                // Une seule traversee centrale ne justifie pas une grande nappe de terrain.
                // A center crossing alone cannot justify a broad patch that bridges horizon relief.
                // Require a source-field crossing at every coarse-cell sample before emitting it.
                if(!bValidPatch)continue;
                const int32 Base=static_cast<int32>(NewVertices.size());
                for(const FPatchVertex& Vertex:Patch)
                {
                    NewVertices.emplace_back(Vertex.X*25.0,Vertex.Y*25.0,Vertex.Z*25.0);
                    NewNormals.push_back(Vertex.Normal);
                    if(bHasUVs)NewUVs.emplace_back(static_cast<float>((Vertex.X-Origin.X)/256.0),
                        static_cast<float>((Vertex.Y-Origin.Y)/256.0));
                    if(bHasColors)NewColors.push_back(Color);
                }
                const bool bBelowSolid=Density(X,Y,Crossing-1.0)<0.0f;
                for(int32 PY=0;PY<PatchDivisions;++PY)
                for(int32 PX=0;PX<PatchDivisions;++PX)
                {
                    const int32 I0=Base+PY*PatchSide+PX;
                    const int32 I1=I0+1;
                    const int32 I2=I1+PatchSide;
                    const int32 I3=I0+PatchSide;
                    if(bBelowSolid)NewTriangles.insert(NewTriangles.end(),{I0,I1,I3,I1,I2,I3});
                    else NewTriangles.insert(NewTriangles.end(),{I0,I3,I1,I1,I3,I2});
                }
            }
        }
    }
    if(NewVertices.empty()||NewTriangles.empty())return;
    const int32 AddedVertices=static_cast<int32>(NewVertices.size());
    auto InsertVertices=[&](auto& Existing,const auto& Added)
    {
        using TArrayType=std::remove_reference_t<decltype(Existing)>;
        TArrayType Combined;Combined.Reserve(Existing.Num()+static_cast<int32>(Added.size()));
        for(int32 I=0;I<OldSurfaceVertices;++I)Combined.Add(Existing[I]);
        for(const auto& Value:Added)Combined.Add(Value);
        for(int32 I=OldSurfaceVertices;I<Existing.Num();++I)Combined.Add(Existing[I]);
        Existing=std::move(Combined);
    };
    InsertVertices(Mesh.Vertices,NewVertices);InsertVertices(Mesh.Normals,NewNormals);
    if(bHasUVs)InsertVertices(Mesh.UVs,NewUVs);
    if(bHasColors)InsertVertices(Mesh.Colors,NewColors);
    for(int32& Index:Mesh.Triangles)if(Index>=OldSurfaceVertices)Index+=AddedVertices;
    const int32 CeilingStart=std::max(0,Mesh.Triangles.Num()-Mesh.NumCeilingTriangles*3);
    TArray<int32> CombinedTriangles;CombinedTriangles.Reserve(Mesh.Triangles.Num()+static_cast<int32>(NewTriangles.size()));
    for(int32 I=0;I<CeilingStart;++I)CombinedTriangles.Add(Mesh.Triangles[I]);
    for(const int32 Index:NewTriangles)CombinedTriangles.Add(OldSurfaceVertices+Index);
    for(int32 I=CeilingStart;I<Mesh.Triangles.Num();++I)CombinedTriangles.Add(Mesh.Triangles[I]);
    const int32 OldTriangleCount=Mesh.Triangles.Num()/3;
    TArray<uint8> ExistingFaceMasks=Mesh.TriangleSkirtFaceMasks;
    ExistingFaceMasks.SetNum(OldTriangleCount);
    TArray<uint8> CombinedFaceMasks;CombinedFaceMasks.Reserve(
        OldTriangleCount+static_cast<int32>(NewTriangles.size()/3));
    for(int32 Triangle=0;Triangle<CeilingStart/3;++Triangle)
        CombinedFaceMasks.Add(ExistingFaceMasks[Triangle]);
    for(size_t Triangle=0;Triangle<NewTriangles.size()/3;++Triangle)CombinedFaceMasks.Add(0);
    for(int32 Triangle=CeilingStart/3;Triangle<OldTriangleCount;++Triangle)
        CombinedFaceMasks.Add(ExistingFaceMasks[Triangle]);
    Mesh.Triangles=std::move(CombinedTriangles);
    Mesh.TriangleSkirtFaceMasks=std::move(CombinedFaceMasks);
    Mesh.NumSurfaceVertices=OldSurfaceVertices+AddedVertices;
}