VoxelForge/Standalone/Source/StandaloneShapeWorld.h

#pragma once

#include "OwnerWorldLoader.h"
#include "StandaloneResearchTunnels.h"
#include "VoxelCaveMorphology.h"
#include "VoxelNoise.h"
#include "VoxelStrateManager.h"

#include <algorithm>
#include <array>
#include <atomic>
#include <bit>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <limits>
#include <memory>
#include <mutex>
#include <immintrin.h>
#include <thread>
#include <type_traits>
#include <unordered_map>
#include <vector>

// Standalone-only world compiler. It compiles the runtime room/tunnel/passages into
// axis-line intervals, composes rock/cave signs into packed Z-column words, then displaces
// only the conservative surface band. The copied Unreal generator remains the A/B reference.
class FStandaloneShapeWorld
{
public:
    enum class EDetailVariant : uint8 { Coherent, Smooth, PerVoxel };
    struct FEditBrush
    {
        enum class EShape : uint8 { Sphere, Capsule, RoundedBox };
        enum class EOperation : uint8 { Carve, AddRock };
        EShape Shape=EShape::Sphere;
        EOperation Operation=EOperation::Carve;
        std::array<double,3> A{};
        std::array<double,3> B{};
        std::array<double,3> HalfExtent{};
        double Radius=0.0;
        double Rounding=0.0;
    };

    struct FParallelTiming
    {
        double WallSeconds = 0.0;
        double WorkerSeconds = 0.0;
    };

    struct FDetailWorkerProfile
    {
        uint64 Lines=0,Candidates=0,DetailEvaluations=0,Flips=0;
        uint64 LocalBandSkips=0;
        uint64 DistanceSamples=0,DistanceSampleNanoseconds=0;
        uint64 DetailSamples=0,DetailSampleNanoseconds=0;
        uint64 ReliefSampleNanoseconds=0,LayerLookupSampleNanoseconds=0;
        uint64 OwnerRoughnessSampleNanoseconds=0,OwnerRoughnessFieldCalls=0;
        uint64 OwnerRoughnessActiveCalls=0,OwnerNoiseCacheHits=0,OwnerNoiseCacheMisses=0;
    };

    struct FBuildStats
    {
        double ShapeSeconds = 0.0;
        double IntervalSeconds = 0.0;
        double SurfaceSeconds = 0.0;
        double FarShapeSeconds = 0.0;
        uint64 PrimitiveCount = 0;
        uint64 RoomCount = 0;
        uint64 RoomFloorJoinCount = 0;
        uint64 TunnelCount = 0;
        uint64 PassageCount = 0;
        uint64 SignNodes = 0;
        uint64 SurfaceCandidates = 0;
        uint64 SurfaceProxyNodes = 0;
        uint64 DetailFlips = 0;
        uint64 WalkSamples = 0;
        uint64 WalkFailures = 0;
        double MaxCapsuleJoinRadiusDelta = 0.0;
        double MaxCapsuleJoinRadiusSlopeDelta = 0.0;
        double MaxCapsuleJoinFloorDelta = 0.0;
        double MaxCapsuleJoinFloorSlopeDelta = 0.0;
        double MaxRawUnmarkedFloorSlopeDelta = 0.0;
        double MaxUnmarkedFloorSlopeDelta = 0.0;
        double MaxUnmarkedFloorDirectionSlopeDelta = 0.0;
        double MaxMarkedFloorSlopeDelta = 0.0;
        double MaxFloorCurveDeviation = 0.0;
        double LocalPrimitiveSeconds = 0.0;
        double LocalLineBinSeconds = 0.0;
        double CandidateMergeSeconds = 0.0;
        double FarPrimitiveSeconds = 0.0;
        double FarBinSeconds = 0.0;
        FParallelTiming RockIntervalLines;
        FParallelTiming CandidateAxisLines;
        FParallelTiming CandidateMergeLines;
        FParallelTiming DetailBandLines;
        uint64 CurvedFloorSegments = 0;
        uint64 MarkedFloorSlopeJoins = 0;
        double LegacySegmentRadiusStep = 0.0;
        double SmoothedProfileRadiusStep = 0.0;
        double MaxTunnelBendDegrees = 0.0;
        FVector MaxTunnelBendVoxels = FVector::ZeroVector;
        FVector MaxTunnelBendDirection = FVector::ZeroVector;
        uint64 CapsuleJoinCount = 0;
    };

    struct FPrimitive
    {
        int32 CenterSub[3]{};
        int32 CapsuleStartSub[3]{};
        int32 CapsuleEndSub[3]{};
        double CapsuleStartWorld[3]{};
        double CapsuleAxisDelta[3]{};
        double CapsuleAxisLengthSquared = 0.0;
        double CapsuleAxisLength = 0.0;
        int64 Q[3][3]{};
        int32 Min[3]{};
        int32 Max[3]{};
        float MinRadius = 1.0f;
        float CapsuleStartRadius = 1.0f;
        float CapsuleEndRadius = 1.0f;
        float CapsuleProfileStartRadius = 1.0f;
        float CapsuleProfileEndRadius = 1.0f;
        double CapsuleProfileStartT = 0.0;
        double CapsuleProfileEndT = 1.0;
        double CapsuleRadiusSlopeStart = 0.0;
        double CapsuleRadiusSlopeEnd = 0.0;
        double BlendRadius = 0.0;
        double FloorBlendRadius = 0.0;
        bool bCarve = true;
        bool bCritical = false;
        bool bCapsule = false;
        bool bSmoothCapsuleRadius = false;
        bool bRoundedBox = false;
        bool bHasFloorClip = false;
        bool bCurvedFloorClip = false;
        double RoundedHalfExtent[3]{};
        double Rounding = 0.0;
        double FloorCenterZ = 0.0;
        double FloorSlopeX = 0.0;
        double FloorSlopeY = 0.0;
        double FloorCurveOriginX = 0.0;
        double FloorCurveOriginY = 0.0;
        double FloorCurveDeltaX = 0.0;
        double FloorCurveDeltaY = 0.0;
        double FloorCurveLength = 0.0;
        double FloorCurveStartZ = 0.0;
        double FloorCurveEndZ = 0.0;
        double FloorCurveTangentStart = 0.0;
        double FloorCurveTangentEnd = 0.0;
        double FloorCurveStartBlend = 0.0;
        double FloorCurveEndBlend = 0.0;
        double FloorCurveMaxDeviation = 0.0;
    };

    struct FInterval { int32 Lo = 0, Hi = -1; };
    struct FBinKey { int32 X=0,Y=0,Z=0; bool operator==(const FBinKey&)const=default; };
    struct FBinHash
    {
        size_t operator()(const FBinKey& K) const noexcept
        {
            uint64 H=static_cast<uint32>(K.X)*0x9E3779B185EBCA87ull;
            H^=static_cast<uint64>(static_cast<uint32>(K.Y))*0xC2B2AE3D27D4EB4Full;
            H^=static_cast<uint64>(static_cast<uint32>(K.Z))*0x165667B19E3779F9ull;
            H^=H>>29;H*=0xBF58476D1CE4E5B9ull;H^=H>>31;
            return static_cast<size_t>(H);
        }
    };
    struct FBinItems
    {
        static constexpr size_t InlineCapacity=4;
        struct FHeapStorage {uint32* Data=nullptr;uint32 Capacity=0;uint32 Reserved=0;};
        union FStorage
        {
            std::array<uint32,InlineCapacity> Inline;
            FHeapStorage Heap;
            FStorage():Inline{}{}
            ~FStorage(){}
        } Storage;
        uint32 Count=0;
        uint32 bHeap=0;

        FBinItems()=default;
        FBinItems(const FBinItems&)=delete;
        FBinItems& operator=(const FBinItems&)=delete;
        FBinItems(FBinItems&& Other) noexcept {MoveFrom(Other);}
        FBinItems& operator=(FBinItems&& Other) noexcept
        {
            if(this!=&Other){Release();MoveFrom(Other);}
            return *this;
        }
        ~FBinItems(){Release();}

        void push_back(uint32 Id)
        {
            if(!bHeap&&Count<InlineCapacity)
            {
                Storage.Inline[Count++]=Id;
                return;
            }
            if(!bHeap)
            {
                const std::array<uint32,InlineCapacity> Previous=Storage.Inline;
                Storage.Heap={new uint32[8],8,0};
                std::copy(Previous.begin(),Previous.end(),Storage.Heap.Data);
                bHeap=1;
            }
            else if(Count==Storage.Heap.Capacity)
            {
                const uint32 NewCapacity=Storage.Heap.Capacity*2;
                uint32* NewData=new uint32[NewCapacity];
                std::copy_n(Storage.Heap.Data,Count,NewData);
                delete[] Storage.Heap.Data;
                Storage.Heap={NewData,NewCapacity,0};
            }
            Storage.Heap.Data[Count++]=Id;
        }
        size_t size() const {return Count;}
        size_t capacity() const {return bHeap?Storage.Heap.Capacity:InlineCapacity;}
        const uint32* begin() const {return bHeap?Storage.Heap.Data:Storage.Inline.data();}
        const uint32* end() const {return begin()+Count;}

    private:
        void Release()
        {
            if(bHeap)delete[] Storage.Heap.Data;
            bHeap=0;Count=0;
        }
        void MoveFrom(FBinItems& Other)
        {
            Count=Other.Count;bHeap=Other.bHeap;
            if(bHeap)
            {
                Storage.Heap=Other.Storage.Heap;
                Other.bHeap=0;Other.Count=0;
            }
            else
            {
                Storage.Inline=Other.Storage.Inline;
                Other.Count=0;
            }
        }
    };

    static constexpr int32 QScale = 1 << 28;
    static constexpr int32 SubPerVoxel = 16;
    static constexpr int32 DetailBand = 6;
    static constexpr float DetailSupremum = 4.0f;
    static constexpr double JoinBlendRadius = 1.5;
    static constexpr double RoomUnionBlendRadius = 12.0;
    static constexpr double FloorClipBlendRadius = 16.0;
    static constexpr float StrataCaveBlendRadius = 12.0f;
    static constexpr int32 LineBinSize = 32;
    static constexpr int32 SpatialBinSize = 256;

    static void SetDetailProfileEnabled(bool bEnabled)
    {ProfileDetailBand.store(bEnabled,std::memory_order_release);}

    bool Build(FStandaloneOwnerWorld& World, const FIntVector& MinVoxels,
               const FIntVector& MaxVoxels, const FIntVector& ClipMinVoxels,
               const FIntVector& ClipMaxVoxels, int32 Workers, EDetailVariant Variant,
               bool bIncludeFarWorld,
               std::string& Error)
    {
        using Clock=std::chrono::steady_clock;
        const auto TotalStart=Clock::now();
        Owner=&World;WorkerCount=std::clamp(Workers,1,16);DetailVariant=Variant;
        for(int32 A=0;A<3;++A)
        {
            Min[A]=(&MinVoxels.X)[A]-2;
            Max[A]=(&MaxVoxels.X)[A]+2;
            ClipMin[A]=(&ClipMinVoxels.X)[A];
            ClipMax[A]=(&ClipMaxVoxels.X)[A];
            Size[A]=Max[A]-Min[A]+1;
            if(Size[A]<=0||Size[A]>4096||ClipMax[A]<ClipMin[A])
            {Error="shape interval region is invalid or outside the bounded 4096-node compiler envelope";return false;}
        }
        const uint64 NodeCount=static_cast<uint64>(Size[0])*Size[1]*Size[2];
        if(NodeCount>160000000ull)
        {Error="shape interval region exceeds the 160-million-node safety bound";return false;}
        Stats={};Stats.SignNodes=NodeCount;

        auto Stage=Clock::now();
        BuildLayerBounds();
        auto Substage=Clock::now();
        BuildOwnerPrimitives();
        Stats.LocalPrimitiveSeconds=std::chrono::duration<double>(Clock::now()-Substage).count();
        Substage=Clock::now();
        if(FStandaloneResearchTunnels::Enabled.load(std::memory_order_acquire))
            BuildRoundTunnelPrimitives();
        BuildSpatialBins(LineBinSize);
        Stats.LocalLineBinSeconds=std::chrono::duration<double>(Clock::now()-Substage).count();
        Stats.ShapeSeconds=std::chrono::duration<double>(Clock::now()-Stage).count();
        Stats.PrimitiveCount=Primitives.size();
        if(Primitives.empty())
        {Error="owner shape extraction produced no rooms, tunnels, or passages";return false;}
        Stage=Clock::now();
        const size_t NZ=static_cast<size_t>(Size[2]);
        WordsPerZ=(NZ+63u)/64u;
        RockWords.assign(static_cast<size_t>(Size[0])*Size[1]*WordsPerZ,0);
        SurfaceDistanceQ=std::make_unique_for_overwrite<int8[]>(static_cast<size_t>(NodeCount));
        BuildPackedRockIntervals(Stats.RockIntervalLines);
        BuildCandidateBands(Stats.CandidateAxisLines,Stats.CandidateMergeLines,
            Stats.CandidateMergeSeconds);
        Stats.IntervalSeconds=std::chrono::duration<double>(Clock::now()-Stage).count();

        Stage=Clock::now();
        ApplyDetailBand(Stats.DetailBandLines);
        Stats.SurfaceSeconds=std::chrono::duration<double>(Clock::now()-Stage).count();
        ValidatePassageRoutes();
        Built=true;
        if(bIncludeFarWorld)BuildFarWorld();
        const uint64 WordBytes=static_cast<uint64>(RockWords.size()+CandidateWords.size())*sizeof(uint64);
        uint64 CandidateAxisBytes=0;
        for(const auto& AxisWords:CandidateAxisWords)
            CandidateAxisBytes+=static_cast<uint64>(AxisWords.capacity())*sizeof(uint64);
        uint64 BinReferences=0,BinCapacity=0,InlineBins=0,OverflowBins=0;
        for(const auto& Bin:SpatialBins)
        {
            BinReferences+=Bin.second.size();BinCapacity+=Bin.second.capacity();
            if(Bin.second.size()<=FBinItems::InlineCapacity)++InlineBins;
            else ++OverflowBins;
        }
        std::printf("[ShapeProfile] local_primitives=%.3fms line_bins=%.3fms rock_lines_wall=%.3fms rock_lines_worker=%.3fms candidate_axes_wall=%.3fms candidate_axes_worker=%.3fms candidate_merge_wall=%.3fms candidate_merge_worker=%.3fms candidate_merge_serial=%.3fms detail_lines_wall=%.3fms detail_lines_worker=%.3fms\n",
            Stats.LocalPrimitiveSeconds*1000.0,Stats.LocalLineBinSeconds*1000.0,
            Stats.RockIntervalLines.WallSeconds*1000.0,Stats.RockIntervalLines.WorkerSeconds*1000.0,
            Stats.CandidateAxisLines.WallSeconds*1000.0,Stats.CandidateAxisLines.WorkerSeconds*1000.0,
            Stats.CandidateMergeLines.WallSeconds*1000.0,Stats.CandidateMergeLines.WorkerSeconds*1000.0,
            Stats.CandidateMergeSeconds*1000.0,Stats.DetailBandLines.WallSeconds*1000.0,
            Stats.DetailBandLines.WorkerSeconds*1000.0);
        std::printf("[ShapeAllocation] words=%.2fMiB candidate_axis_words=%.2fMiB proxy=%.2fMiB primitives=%.2fMiB spatial_bins=%llu spatial_bin_refs=%llu spatial_bin_capacity=%llu inline_bins=%llu overflow_bins=%llu far_primitive_build=%.3fms far_bin_rebuild=%.3fms\n",
            static_cast<double>(WordBytes)/1048576.0,static_cast<double>(CandidateAxisBytes)/1048576.0,
            static_cast<double>(Stats.SignNodes*sizeof(int8))/1048576.0,
            static_cast<double>(Primitives.capacity()*sizeof(FPrimitive))/1048576.0,
            static_cast<unsigned long long>(SpatialBins.size()),
            static_cast<unsigned long long>(BinReferences),
            static_cast<unsigned long long>(BinCapacity),
            static_cast<unsigned long long>(InlineBins),
            static_cast<unsigned long long>(OverflowBins),Stats.FarPrimitiveSeconds*1000.0,
            Stats.FarBinSeconds*1000.0);
        std::printf("[ShapeCompile] primitives=%llu rooms=%llu tunnels=%llu passages=%llu nodes=%llu packed=%.2f MiB surface_proxy_nodes=%llu proxy_address_capacity=%.2f MiB workers=%d\n",
            static_cast<unsigned long long>(Stats.PrimitiveCount),
            static_cast<unsigned long long>(Stats.RoomCount),
            static_cast<unsigned long long>(Stats.TunnelCount),
            static_cast<unsigned long long>(Stats.PassageCount),
            static_cast<unsigned long long>(Stats.SignNodes),
            static_cast<double>(WordBytes)/1048576.0,
            static_cast<unsigned long long>(Stats.SurfaceProxyNodes),
            static_cast<double>(Stats.SignNodes*sizeof(int8))/1048576.0,WorkerCount);
        std::printf("[ShapeStage] local-shapes=%.3f ms intervals+signs=%.3f ms surface-band=%.3f ms far-parents=%.3f ms total=%.3f ms candidates=%llu flips=%llu walk=%llu/%llu\n",
            Stats.ShapeSeconds*1000.0,Stats.IntervalSeconds*1000.0,Stats.SurfaceSeconds*1000.0,Stats.FarShapeSeconds*1000.0,
            std::chrono::duration<double>(Clock::now()-TotalStart).count()*1000.0,
            static_cast<unsigned long long>(Stats.SurfaceCandidates),
            static_cast<unsigned long long>(Stats.DetailFlips),
            static_cast<unsigned long long>(Stats.WalkSamples-Stats.WalkFailures),
            static_cast<unsigned long long>(Stats.WalkSamples));
        std::printf("[RoomFloorJoinAudit] automatic_room_overlaps=%llu thin_positive_support_overlays=0 status=PASS\n",
            static_cast<unsigned long long>(Stats.RoomFloorJoinCount));
        const double SmoothUnionGrowthBound=JoinBlendRadius/6.0;
        const bool bFloorSlopePass=Stats.MaxUnmarkedFloorSlopeDelta<=0.75;
        const bool bFloorBlendPass=Stats.MaxFloorCurveDeviation<=2.0;
        std::printf("[FloorSlopeJoinAudit] raw_unmarked_max_delta=%.5f vox generated_tangent_delta=%.5f vox bend_direction_slope_change=%.5f vox marked_max_delta=%.5f vox marked_joins=%llu threshold=0.750 status=%s\n",
            Stats.MaxRawUnmarkedFloorSlopeDelta,Stats.MaxUnmarkedFloorSlopeDelta,
            Stats.MaxUnmarkedFloorDirectionSlopeDelta,Stats.MaxMarkedFloorSlopeDelta,
            static_cast<unsigned long long>(Stats.MarkedFloorSlopeJoins),bFloorSlopePass?"PASS":"FAIL");
        std::printf("[FloorProfileBlend] c1_segments=%llu max_curve_deviation=%.3f vox threshold=2.000 status=%s\n",
            static_cast<unsigned long long>(Stats.CurvedFloorSegments),Stats.MaxFloorCurveDeviation,
            bFloorBlendPass?"PASS":"FAIL");
        std::printf("[TunnelRadiusAudit] shared_joins=%llu legacy_fixed_segment_max_step=%.4f vox smoothed_profile_max_step=%.4f vox tapered_endpoint_max_delta=%.6f vox radius_slope_join_delta=%.8f floor_endpoint_max_delta=%.6f vox floor_slope_join_delta=%.5f unmarked_floor_slope_delta=%.5f marked_floor_slope_delta=%.5f marked_slope_joins=%llu smooth_union_growth_bound=%.4f vox thresholds=(0.005,0.000001,0.005,0.25) status=%s\n",
            static_cast<unsigned long long>(Stats.CapsuleJoinCount),Stats.LegacySegmentRadiusStep,
            Stats.SmoothedProfileRadiusStep,Stats.MaxCapsuleJoinRadiusDelta,
            Stats.MaxCapsuleJoinRadiusSlopeDelta,
            Stats.MaxCapsuleJoinFloorDelta,Stats.MaxCapsuleJoinFloorSlopeDelta,
            Stats.MaxUnmarkedFloorSlopeDelta,Stats.MaxMarkedFloorSlopeDelta,
            static_cast<unsigned long long>(Stats.MarkedFloorSlopeJoins),SmoothUnionGrowthBound,
            Stats.MaxCapsuleJoinRadiusDelta<=0.005&&Stats.MaxCapsuleJoinRadiusSlopeDelta<=1.0e-6
                &&Stats.MaxCapsuleJoinFloorDelta<=0.005
                &&SmoothUnionGrowthBound<=0.25?"PASS":"FAIL");
        const double BendHorizontal=std::sqrt(Stats.MaxTunnelBendDirection.X*Stats.MaxTunnelBendDirection.X
            +Stats.MaxTunnelBendDirection.Y*Stats.MaxTunnelBendDirection.Y);
        const double BendYaw=std::atan2(Stats.MaxTunnelBendDirection.X,Stats.MaxTunnelBendDirection.Y)*57.29577951308232;
        const double BendPitch=std::atan2(Stats.MaxTunnelBendDirection.Z,BendHorizontal)*57.29577951308232;
        std::printf("[TunnelBendAudit] sharpest_bend=%.2f deg centerline=(%.2f,%.2f,%.2f) vox tangent=(%.4f,%.4f,%.4f) camera_yaw=%.2f camera_pitch=%.2f\n",
            Stats.MaxTunnelBendDegrees,Stats.MaxTunnelBendVoxels.X,
            Stats.MaxTunnelBendVoxels.Y,Stats.MaxTunnelBendVoxels.Z,
            Stats.MaxTunnelBendDirection.X,Stats.MaxTunnelBendDirection.Y,
            Stats.MaxTunnelBendDirection.Z,BendYaw,BendPitch);
        std::fflush(stdout);
        return true;
    }

    const FBuildStats& GetStats() const { return Stats; }
    bool IsBuilt() const { return Built; }
    bool HasFarWorld() const { return FarBuilt; }
    bool DebugCriticalAirAt(double X,double Y,double Z) const { return IsInsideCriticalAir(X,Y,Z); }
    bool DebugSolidOverlayAt(double X,double Y,double Z) const { return IsInsideSolidOverlay(X,Y,Z); }

    void SetEditTransaction(std::vector<FEditBrush> Brushes)
    {
        EditBrushes=std::move(Brushes);
        EditBoundsMin={MAX_int32,MAX_int32,MAX_int32};
        EditBoundsMax={MIN_int32,MIN_int32,MIN_int32};
        for(const FEditBrush& Brush:EditBrushes)
        {
            double Extent[3]{};
            if(Brush.Shape==FEditBrush::EShape::Sphere)
                for(int32 Axis=0;Axis<3;++Axis)Extent[Axis]=Brush.Radius;
            else if(Brush.Shape==FEditBrush::EShape::Capsule)
                for(int32 Axis=0;Axis<3;++Axis)Extent[Axis]=Brush.Radius;
            else
                for(int32 Axis=0;Axis<3;++Axis)Extent[Axis]=Brush.HalfExtent[Axis]+Brush.Rounding;
            for(int32 Axis=0;Axis<3;++Axis)
            {
                const double Lo=Brush.Shape==FEditBrush::EShape::Capsule
                    ?std::min(Brush.A[Axis],Brush.B[Axis])-Extent[Axis]
                    :Brush.A[Axis]-Extent[Axis];
                const double Hi=Brush.Shape==FEditBrush::EShape::Capsule
                    ?std::max(Brush.A[Axis],Brush.B[Axis])+Extent[Axis]
                    :Brush.A[Axis]+Extent[Axis];
                (&EditBoundsMin.X)[Axis]=std::min((&EditBoundsMin.X)[Axis],static_cast<int32>(std::floor(Lo)));
                (&EditBoundsMax.X)[Axis]=std::max((&EditBoundsMax.X)[Axis],static_cast<int32>(std::ceil(Hi)));
            }
        }
    }

    bool HasEditTransaction() const { return !EditBrushes.empty(); }
    FIntVector GetEditBoundsMin() const { return EditBoundsMin; }
    FIntVector GetEditBoundsMax() const { return EditBoundsMax; }

    void BuildFarWorld()
    {
        if(FarBuilt)return;
        auto Stage=std::chrono::steady_clock::now();
        BuildFarParentShapes();
        Stats.FarPrimitiveSeconds=std::chrono::duration<double>(std::chrono::steady_clock::now()-Stage).count();
        Stage=std::chrono::steady_clock::now();
        BuildSpatialBins(SpatialBinSize,false);
        Stats.FarBinSeconds=std::chrono::duration<double>(std::chrono::steady_clock::now()-Stage).count();
        Stats.FarShapeSeconds=Stats.FarPrimitiveSeconds+Stats.FarBinSeconds;
        Stats.PrimitiveCount=Primitives.size();FarBuilt=true;
        std::printf("[FarShapeStage] parent_shapes=%.3f ms spatial_bins=%.3f ms total=%.3f ms primitives=%llu\n",
            Stats.FarPrimitiveSeconds*1000.0,Stats.FarBinSeconds*1000.0,Stats.FarShapeSeconds*1000.0,
            static_cast<unsigned long long>(Stats.PrimitiveCount));
        std::fflush(stdout);
    }

    // Lattice signs come from packed intervals; a bounded analytic distance proxy supplies
    // edge interpolation and normals to the existing mesher.
    float DensityAt(float X,float Y,float Z) const
    {
        int32 IX=0,IY=0,IZ=0;
        const bool bInteger=NearInteger(X,IX)&&NearInteger(Y,IY)&&NearInteger(Z,IZ);
        const bool bPacked=bInteger&&Contains(IX,IY,IZ);
        if(bPacked)
            return PackedDensityAt(IX,IY,IZ,IsRock(IX,IY,IZ),X,Y,Z);
        return AnalyticDensityAt(X,Y,Z);
    }

    float DensityAtLattice(int32 X,int32 Y,int32 Z) const
    {
        const bool bPacked=Contains(X,Y,Z);
        if(bPacked)
            return PackedDensityAt(X,Y,Z,IsRock(X,Y,Z),static_cast<float>(X),
                static_cast<float>(Y),static_cast<float>(Z));
        return AnalyticDensityAt(static_cast<float>(X),static_cast<float>(Y),
            static_cast<float>(Z));
    }

private:
    float PackedDensityAt(int32 X,int32 Y,int32 Z,bool bRock,float FX,float FY,float FZ) const
    {
        const size_t Column=static_cast<size_t>(Y-Min[1])*Size[0]+static_cast<size_t>(X-Min[0]);
        const size_t Zi=static_cast<size_t>(Z-Min[2]);
        const size_t Word=Column*WordsPerZ+(Zi>>6);
        const uint64 Mask=uint64{1}<<(Zi&63u);
        if(CandidateWords[Word]&Mask)
        {
            const size_t Node=Zi*static_cast<size_t>(Size[1])*Size[0]+Column;
            const int8 Proxy=SurfaceDistanceQ[Node];
            float Density=static_cast<float>(Proxy)*(1.0f/16.0f);
            if((Density<0.0f)!=bRock)
                Density=(bRock?-1.0f:1.0f)*std::max(std::abs(Density),1.0f/16.0f);
            return ApplyEditBrushes(Density,FX,FY,FZ);
        }
        // Match the fixed fallback to the edge of the stored proxy band. The
        // displacement is tapered to zero there, so this avoids a magnitude
        // jump along the candidate bit-mask boundary.
        constexpr float ProxyBandEdge=DetailSupremum+0.125f;
        return ApplyEditBrushes(bRock?-ProxyBandEdge:ProxyBandEdge,FX,FY,FZ);
    }

    float AnalyticDensityAt(float X,float Y,float Z) const
    {
        const bool bInsideFine=X>=Min[0]&&X<=Max[0]&&Y>=Min[1]&&Y<=Max[1]
            &&Z>=Min[2]&&Z<=Max[2];
        if(!bInsideFine)
        {
            bool bCriticalAir=false;
            float RockDistance=CoarseRockDistance(X,Y,Z,&bCriticalAir);
            if(bCriticalAir)RockDistance=std::min(RockDistance,-1.0f/16.0f);
            const bool bRock=RockDistance>=0.0f;
            float Density=-RockDistance;
            if((Density<0.0f)!=bRock||std::abs(Density)<0.05f)
                Density=(bRock?-1.0f:1.0f)*std::max(std::abs(Density),0.05f);
            return ApplyEditBrushes(Density,X,Y,Z);
        }
        ESurfaceFamily SurfaceFamily=ESurfaceFamily::Strata;float SurfaceSDF=100000.0f;
        bool bFloorSurface=false,bCriticalAir=false;float CaveInfluence=0.0f,FloorInfluence=0.0f;
        float RockDistance=EffectiveRockDistance(X,Y,Z,&SurfaceFamily,&SurfaceSDF,
            &bFloorSurface,&CaveInfluence,&FloorInfluence,&bCriticalAir);
        if(bInsideFine&&DetailVariant!=EDetailVariant::Smooth
            &&std::abs(RockDistance)<=DetailSupremum+0.125f)
        {
            const float Fade=DetailBandFadeAt(X,Y,Z,RockDistance,SurfaceSDF,
                CaveInfluence,FloorInfluence);
            if(Fade>0.0f)RockDistance+=DetailNoise(X,Y,Z,SurfaceSDF,CaveInfluence,FloorInfluence)*Fade;
        }
        if(bCriticalAir)RockDistance=std::min(RockDistance,-1.0f/16.0f);
        const bool bRock=RockDistance>=0.0f;
        float Density=-RockDistance;
        if((Density<0.0f)!=bRock||std::abs(Density)<0.05f)
            Density=(bRock?-1.0f:1.0f)*std::max(std::abs(Density),0.05f);
        return ApplyEditBrushes(Density,X,Y,Z);
    }

public:
    EVoxelTileClass ClassifyBox(const FIntVector& Origin,int32 Step,int32 Cells) const
    {
        const int32 Pad=std::max(Step,1);
        int32 BMin[3]={Origin.X-Pad,Origin.Y-Pad,Origin.Z-Pad};
        int32 BMax[3]={Origin.X+(Cells+1)*Step,Origin.Y+(Cells+1)*Step,Origin.Z+(Cells+1)*Step};
        if(!EditBrushes.empty()&&EditBoundsMin.X<=BMax[0]&&EditBoundsMax.X>=BMin[0]
            &&EditBoundsMin.Y<=BMax[1]&&EditBoundsMax.Y>=BMin[1]
            &&EditBoundsMin.Z<=BMax[2]&&EditBoundsMax.Z>=BMin[2])
            return EVoxelTileClass::Mixed;
        for(int32 BZ=FloorDiv(BMin[2],ActiveSpatialBinSize);BZ<=FloorDiv(BMax[2],ActiveSpatialBinSize);++BZ)
        for(int32 BY=FloorDiv(BMin[1],ActiveSpatialBinSize);BY<=FloorDiv(BMax[1],ActiveSpatialBinSize);++BY)
        for(int32 BX=FloorDiv(BMin[0],ActiveSpatialBinSize);BX<=FloorDiv(BMax[0],ActiveSpatialBinSize);++BX)
        {
            const auto It=SpatialBins.find({BX,BY,BZ});if(It==SpatialBins.end())continue;
            for(uint32 Id:It->second)
            {
                const FPrimitive& P=Primitives[Id];bool Overlap=true;
                for(int32 A=0;A<3;++A)
                    Overlap &= P.Max[A]+DetailBand>=BMin[A]&&P.Min[A]-DetailBand<=BMax[A];
                if(Overlap)return EVoxelTileClass::Mixed;
            }
        }
        for(const FInterval& R:LayerRanges)
            if(R.Lo<=BMax[2]&&R.Hi>=BMin[2]
                &&(R.Lo>=BMin[2]||R.Hi<=BMax[2]))
                return EVoxelTileClass::Mixed;
        const bool bRock=EvaluateFast((BMin[0]+BMax[0])*0.5,
            (BMin[1]+BMax[1])*0.5,(BMin[2]+BMax[2])*0.5);
        return bRock?EVoxelTileClass::AllSolid:EVoxelTileClass::AllAir;
    }

    bool IsRock(int32 X,int32 Y,int32 Z) const
    {
        if(!Contains(X,Y,Z))return EvaluateFast(X,Y,Z);
        const size_t Line=static_cast<size_t>(Y-Min[1])*Size[0]
            +static_cast<size_t>(X-Min[0]);
        const size_t Zi=static_cast<size_t>(Z-Min[2]);
        return ((RockWords[Line*WordsPerZ+(Zi>>6)]>>(Zi&63u))&1u)!=0;
    }

    bool IsDetailAuditPassing() const
    {
        const double Correlation=DetailAutocorrelation(DetailVariant);
        const bool bFlips=Stats.DetailFlips>=1000;
        const bool bCoherent=Correlation>0.45;
        const bool bFamilies=AuditSurfaceBandFamilies();
        std::printf("[DetailAudit] mode=%s supremum=%.2f vox band=%d changed_signs=%llu adjacent_correlation=%.4f surface_detail_check=%s spatial_scale_check=%s\n",
            DetailVariant==EDetailVariant::Coherent?"coherent":DetailVariant==EDetailVariant::Smooth?"smooth":"per-voxel",
            DetailSupremum,DetailBand,static_cast<unsigned long long>(Stats.DetailFlips),Correlation,
            bFlips&&bFamilies?"PASS":"FAIL",bCoherent?"PASS":"FAIL");
        return bFlips&&bCoherent&&bFamilies;
    }

    bool AuditSurfaceBandFamilies() const
    {
        struct FPatchKey
        {
            int32 X=0,Y=0,Z=0;uint8 Family=0;
            bool operator==(const FPatchKey&)const=default;
        };
        struct FPatchHash
        {
            size_t operator()(const FPatchKey& K)const noexcept
            {
                uint64 H=static_cast<uint32>(K.X)*0x9E3779B185EBCA87ull;
                H^=static_cast<uint64>(static_cast<uint32>(K.Y))*0xC2B2AE3D27D4EB4Full;
                H^=static_cast<uint64>(static_cast<uint32>(K.Z))*0x165667B19E3779F9ull;
                H^=static_cast<uint64>(K.Family)*0x94D049BB133111EBull;
                H^=H>>29;H*=0xBF58476D1CE4E5B9ull;H^=H>>31;
                return static_cast<size_t>(H);
            }
        };
        struct FPatchStats {uint64 Count=0;double Sum=0.0,SumSq=0.0;};
        struct FFamilyStats {uint64 Count=0,NearZero=0;double Sum=0.0,SumSq=0.0;};
        std::unordered_map<FPatchKey,FPatchStats,FPatchHash> Patches;
        std::array<FFamilyStats,3> Families{};
        for(size_t Line=0;Line<static_cast<size_t>(Size[0])*Size[1];++Line)
        {
            const int32 X=Min[0]+static_cast<int32>(Line%Size[0]);
            const int32 Y=Min[1]+static_cast<int32>(Line/Size[0]);
            const size_t WordBase=Line*WordsPerZ;
            for(size_t W=0;W<WordsPerZ;++W)
            {
                uint64 Bits=CandidateWords[WordBase+W];
                while(Bits)
                {
                    const uint32 Bit=std::countr_zero(Bits);Bits&=Bits-1;
                    const size_t Zi=W*64+Bit;if(Zi>=static_cast<size_t>(Size[2]))continue;
                    const int32 Z=Min[2]+static_cast<int32>(Zi);
                    ESurfaceFamily Family=ESurfaceFamily::Strata;float SurfaceSDF=100000.0f;
                    bool bFloorSurface=false;float CaveInfluence=0.0f,FloorInfluence=0.0f;
                    const float Raw=EffectiveRockDistance(X,Y,Z,&Family,&SurfaceSDF,
                        &bFloorSurface,&CaveInfluence,&FloorInfluence);
                    if(std::abs(Raw)>=0.5f)continue;
                    const float Fade=DetailBandFadeAt(X,Y,Z,Raw,SurfaceSDF,
                        CaveInfluence,FloorInfluence);
                    const float Detail=Fade>0.0f
                        ?DetailNoise(X,Y,Z,SurfaceSDF,CaveInfluence,FloorInfluence)*Fade:0.0f;
                    const size_t FamilyIndex=static_cast<size_t>(Family);
                    FFamilyStats& F=Families[FamilyIndex];++F.Count;F.Sum+=Detail;F.SumSq+=Detail*Detail;
                    if(std::abs(Detail)<0.05f)++F.NearZero;
                    const FPatchKey Key{FloorDiv(X,16),FloorDiv(Y,16),FloorDiv(Z,16),
                        static_cast<uint8>(FamilyIndex)};
                    FPatchStats& P=Patches[Key];++P.Count;P.Sum+=Detail;P.SumSq+=Detail*Detail;
                }
            }
        }
        std::array<uint64,3> FlatCounts{};uint64 LargestFlatPatch=0;
        for(const auto& Entry:Patches)
        {
            const FPatchStats& P=Entry.second;if(P.Count<16)continue;
            const double Mean=P.Sum/P.Count;
            const double StdDev=std::sqrt(std::max(0.0,P.SumSq/P.Count-Mean*Mean));
            if(StdDev<0.10)
            {
                FlatCounts[Entry.first.Family]+=P.Count;
                LargestFlatPatch=std::max(LargestFlatPatch,P.Count);
            }
        }
        bool bPass=true;uint64 Tested=0,FlatTotal=0,Total=0;
        for(size_t I=0;I<Families.size();++I)
        {
            const FFamilyStats& F=Families[I];Total+=F.Count;FlatTotal+=FlatCounts[I];
            if(F.Count<16)continue;
            ++Tested;
            const double Mean=F.Sum/F.Count;
            const double StdDev=std::sqrt(std::max(0.0,F.SumSq/F.Count-Mean*Mean));
            const double ZeroShare=static_cast<double>(F.NearZero)/F.Count;
            const bool FamilyPass=StdDev>=0.15&&ZeroShare<0.60&&FlatCounts[I]<F.Count/2;
            bPass&=FamilyPass;
            const char* Name=I==0?"strata":I==1?"cave":"solid_overlay";
            std::printf("[SurfaceBandFamilyAudit] family=%s samples=%llu detail_stddev=%.4f near_zero_share=%.1f%% flat_patch_samples=%llu status=%s\n",
                Name,static_cast<unsigned long long>(F.Count),StdDev,ZeroShare*100.0,
                static_cast<unsigned long long>(FlatCounts[I]),FamilyPass?"PASS":"FAIL");
        }
        const double FlatShare=Total?100.0*static_cast<double>(FlatTotal)/Total:100.0;
        std::printf("[FlatPatchAudit] patch_width=16 vox min_patch_samples=16 flat_stddev_lt=0.10 flat_share=%.1f%% largest_flat_patch=%llu status=%s\n",
            FlatShare,static_cast<unsigned long long>(LargestFlatPatch),
            Tested>0&&FlatShare<20.0?"PASS":"FAIL");
        return Tested>0&&FlatShare<20.0&&bPass;
    }

    bool AuditViewSurface(double CameraX,double CameraY,double CameraZ,
                          float Yaw,float Pitch) const
    {
        struct FFamilyStats
        {
            uint64 Hits=0,NearZero=0,AxisLocked=0;
            double Sum=0.0,SumSq=0.0;
        };
        std::array<FFamilyStats,3> FamilyStats{};
        constexpr std::array<float,13> YawOffsets={-36.0f,-30.0f,-24.0f,-18.0f,-12.0f,-6.0f,
            0.0f,6.0f,12.0f,18.0f,24.0f,30.0f,36.0f};
        constexpr std::array<float,11> PitchOffsets={-80.0f,-70.0f,-60.0f,-50.0f,-40.0f,
            -30.0f,-20.0f,-10.0f,0.0f,10.0f,20.0f};
        const double OX=CameraX*4.0,OY=CameraY*4.0,OZ=CameraZ*4.0;
        auto Distance=[&](double X,double Y,double Z,ESurfaceFamily* OutFamily=nullptr,
                          float* OutRaw=nullptr,float* OutDetail=nullptr,
                          bool* OutFloorSurface=nullptr)
        {
            ESurfaceFamily Family=ESurfaceFamily::Strata;float SurfaceSDF=100000.0f;
            bool bFloorSurface=false;float CaveInfluence=0.0f,FloorInfluence=0.0f;
            float D=EffectiveRockDistance(X,Y,Z,&Family,&SurfaceSDF,
                &bFloorSurface,&CaveInfluence,&FloorInfluence);
            float Detail=0.0f;
            const bool bInside=X>=Min[0]&&X<=Max[0]&&Y>=Min[1]&&Y<=Max[1]
                &&Z>=Min[2]&&Z<=Max[2];
            const bool bCriticalAir=IsInsideCriticalAir(X,Y,Z);
            if(bInside&&DetailVariant!=EDetailVariant::Smooth&&std::abs(D)<=DetailSupremum+0.125f)
            {
                const float Fade=DetailBandFadeAt(X,Y,Z,D,SurfaceSDF,CaveInfluence,FloorInfluence);
                if(Fade>0.0f)
                {Detail=DetailNoise(X,Y,Z,SurfaceSDF,CaveInfluence,FloorInfluence)*Fade;D+=Detail;}
            }
            if(bCriticalAir)D=std::min(D,-1.0f/16.0f);
            if(OutFamily)*OutFamily=Family;
            if(OutRaw)*OutRaw=D-Detail;
            if(OutDetail)*OutDetail=Detail;
            if(OutFloorSurface)*OutFloorSurface=bFloorSurface;
            return D;
        };
        uint64 TotalHits=0,TotalAxisLocked=0,TotalNearZero=0,FloorTaggedHits=0,FloorAxisLocked=0;
        uint64 FamilyTransitionPairs=0,StableFamilyPairs=0;
        uint64 ThinDetachedSheets=0;
        double ThinnestDetachedSheet=std::numeric_limits<double>::max();
        double FamilyTransitionGradientSq=0.0,StableFamilyGradientSq=0.0;
        for(float YawOffset:YawOffsets)for(float PitchOffset:PitchOffsets)
        {
            const double YawRad=Yaw+YawOffset*0.017453292519943295;
            const double PitchRad=Pitch+PitchOffset*0.017453292519943295;
            const double CP=std::cos(PitchRad),DX=std::sin(YawRad)*CP,
                DY=std::cos(YawRad)*CP,DZ=std::sin(PitchRad);
            double PreviousT=0.0;
            float PreviousD=Distance(OX,OY,OZ+0.001);
            bool bFound=false;
            double PrimaryHitT=-1.0;
            for(double T=0.5;T<=256.0;T+=0.5)
            {
                const double X=OX+DX*T,Y=OY+DY*T,Z=OZ+DZ*T;
                const float D=Distance(X,Y,Z);
                if(PreviousD<0.0f&&D>=0.0f)
                {
                    double Lo=PreviousT,Hi=T;
                    for(int32 It=0;It<10;++It)
                    {
                        const double Mid=(Lo+Hi)*0.5;
                        if(Distance(OX+DX*Mid,OY+DY*Mid,OZ+DZ*Mid)<0.0f)Lo=Mid;else Hi=Mid;
                    }
                    const double HitT=(Lo+Hi)*0.5;
                    PrimaryHitT=HitT;
                    const double XH=OX+DX*HitT,YH=OY+DY*HitT,ZH=OZ+DZ*HitT;
                    ESurfaceFamily Family=ESurfaceFamily::Strata;float Raw=0.0f,Detail=0.0f;
                    bool bFloorSurface=false;
                    Distance(XH,YH,ZH,&Family,&Raw,&Detail,&bFloorSurface);
                    const double GX=Distance(XH+0.5,YH,ZH)-Distance(XH-0.5,YH,ZH);
                    const double GY=Distance(XH,YH+0.5,ZH)-Distance(XH,YH-0.5,ZH);
                    const double GZ=Distance(XH,YH,ZH+0.5)-Distance(XH,YH,ZH-0.5);
                    const double GL=std::sqrt(GX*GX+GY*GY+GZ*GZ);
                    const double AxisDot=GL>1.0e-8?std::max({std::abs(GX),std::abs(GY),std::abs(GZ)})/GL:1.0;
                    for(int32 Axis=0;Axis<3;++Axis)
                    {
                        double AP[3]={XH,YH,ZH},AM[3]={XH,YH,ZH};
                        AP[Axis]+=0.5;AM[Axis]-=0.5;
                        ESurfaceFamily FamilyP=ESurfaceFamily::Strata,FamilyM=ESurfaceFamily::Strata;
                        float DetailP=0.0f,DetailM=0.0f;bool bFloorP=false,bFloorM=false;
                        Distance(AP[0],AP[1],AP[2],&FamilyP,nullptr,&DetailP,&bFloorP);
                        Distance(AM[0],AM[1],AM[2],&FamilyM,nullptr,&DetailM,&bFloorM);
                        const double Gradient=static_cast<double>(DetailP-DetailM);
                        if(FamilyP!=FamilyM||bFloorP!=bFloorM)
                        {FamilyTransitionGradientSq+=Gradient*Gradient;++FamilyTransitionPairs;}
                        else
                        {StableFamilyGradientSq+=Gradient*Gradient;++StableFamilyPairs;}
                    }
                    constexpr double ScanRadius=16.0;
                    constexpr int32 ScanSteps=128;
                    double ScanT=std::max(0.0,HitT-ScanRadius);
                    ESurfaceFamily ScanFamily=ESurfaceFamily::Strata;
                    bool bScanFloor=false;float ScanRaw=0.0f,ScanDetail=0.0f;
                    Distance(OX+DX*ScanT,OY+DY*ScanT,OZ+DZ*ScanT,
                        &ScanFamily,&ScanRaw,&ScanDetail,&bScanFloor);
                    for(int32 Step=1;Step<=ScanSteps;++Step)
                    {
                        const double T=std::min(HitT+ScanRadius,ScanT+0.25);
                        ESurfaceFamily NextFamily=ESurfaceFamily::Strata;
                        bool bNextFloor=false;float NextRaw=0.0f,NextDetail=0.0f;
                        Distance(OX+DX*T,OY+DY*T,OZ+DZ*T,
                            &NextFamily,&NextRaw,&NextDetail,&bNextFloor);
                        if(std::abs(ScanRaw)<=5.0f&&std::abs(NextRaw)<=5.0f)
                        {
                            const double Gradient=static_cast<double>(NextDetail-ScanDetail);
                            if(ScanFamily!=NextFamily||bScanFloor!=bNextFloor)
                            {FamilyTransitionGradientSq+=Gradient*Gradient;++FamilyTransitionPairs;}
                            else
                            {StableFamilyGradientSq+=Gradient*Gradient;++StableFamilyPairs;}
                        }
                        ScanT=T;ScanFamily=NextFamily;bScanFloor=bNextFloor;
                        ScanRaw=NextRaw;ScanDetail=NextDetail;
                    }
                    FFamilyStats& S=FamilyStats[static_cast<size_t>(Family)];
                    ++S.Hits;++TotalHits;S.Sum+=Detail;S.SumSq+=Detail*Detail;
                    if(std::abs(Detail)<0.05f){++S.NearZero;++TotalNearZero;}
                    if(AxisDot>0.985){++S.AxisLocked;++TotalAxisLocked;}
                    if(bFloorSurface){++FloorTaggedHits;if(AxisDot>0.985)++FloorAxisLocked;}
                    if(YawOffset==0.0f&&PitchOffset==0.0f)
                    {
                        const char* Name=Family==ESurfaceFamily::Strata?"strata":
                            Family==ESurfaceFamily::Cave?"cave":"solid_overlay";
                        std::printf("[SurfaceRay] t=%.2f vox hit=(%.2f,%.2f,%.2f) family=%s floor_surface=%s raw_sdf=%.4f detail=%.4f axis_dot=%.4f\n",
                            HitT,XH,YH,ZH,Name,bFloorSurface?"yes":"no",Raw,Detail,AxisDot);
                    }
                    bFound=true;break;
                }
                PreviousT=T;PreviousD=D;
            }
            if(PrimaryHitT>=0.0)
            {
                constexpr double SheetStep=0.25;
                constexpr double MinimumAirGap=3.0;
                constexpr double MaximumSheetThickness=1.0;
                double T=PrimaryHitT+SheetStep;
                float Prev=Distance(OX+DX*T,OY+DY*T,OZ+DZ*T);
                bool bInSolid=Prev>=0.0f;
                double SolidStart=T,LastSolidExit=PrimaryHitT;
                bool bSeparatedFromPrevious=false;
                for(T+=SheetStep;T<=256.0;T+=SheetStep)
                {
                    const float D=Distance(OX+DX*T,OY+DY*T,OZ+DZ*T);
                    if(!bInSolid&&Prev<0.0f&&D>=0.0f)
                    {
                        SolidStart=T;
                        bSeparatedFromPrevious=(T-LastSolidExit)>=MinimumAirGap;
                        bInSolid=true;
                    }
                    else if(bInSolid&&Prev>=0.0f&&D<0.0f)
                    {
                        const double Thickness=T-SolidStart;
                        if(bSeparatedFromPrevious&&Thickness<MaximumSheetThickness)
                        {
                            ++ThinDetachedSheets;
                            ThinnestDetachedSheet=std::min(ThinnestDetachedSheet,Thickness);
                            if(ThinDetachedSheets<=8)
                                std::printf("[ThinSheetCandidate] yaw_offset=%.0f pitch_offset=%.0f entry=%.2f vox thickness=%.3f vox air_gap=%.2f vox\n",
                                    YawOffset,PitchOffset,SolidStart,Thickness,SolidStart-LastSolidExit);
                        }
                        LastSolidExit=T;
                        bSeparatedFromPrevious=false;
                        bInSolid=false;
                    }
                    Prev=D;
                }
            }
            (void)bFound;
        }
        const double NearZeroShare=TotalHits?static_cast<double>(TotalNearZero)/TotalHits:1.0;
        const double AxisShare=TotalHits?static_cast<double>(TotalAxisLocked)/TotalHits:1.0;
        const double FloorAxisShare=FloorTaggedHits?static_cast<double>(FloorAxisLocked)/FloorTaggedHits:1.0;
        const double TransitionRMS=FamilyTransitionPairs
            ?std::sqrt(FamilyTransitionGradientSq/FamilyTransitionPairs):0.0;
        const double StableRMS=StableFamilyPairs
            ?std::sqrt(StableFamilyGradientSq/StableFamilyPairs):0.0;
        const double TransitionRatio=StableRMS>1.0e-6?TransitionRMS/StableRMS:1000000.0;
        double Sum=0.0,SumSq=0.0;
        for(const FFamilyStats& S:FamilyStats){Sum+=S.Sum;SumSq+=S.SumSq;}
        const double Mean=TotalHits?Sum/TotalHits:0.0;
        const double Variance=TotalHits?std::max(0.0,SumSq/TotalHits-Mean*Mean):0.0;
        const double SurfaceDetailRMS=std::sqrt(Variance);
        const double RelativeJoinStep=SurfaceDetailRMS>1.0e-6?TransitionRMS/SurfaceDetailRMS:1000000.0;
        const bool bDetailPass=TotalHits>=24&&std::sqrt(Variance)>=0.15&&NearZeroShare<0.50;
        constexpr double GridStairThreshold=0.40;
        const bool bJoinCoverage=FloorTaggedHits>=8
            &&(FamilyTransitionPairs==0||(FamilyTransitionPairs>=8&&RelativeJoinStep<GridStairThreshold));
        const bool bGridPass=TotalHits>=24&&AxisShare<0.60
            &&bJoinCoverage;
        const bool bThinSheetPass=ThinDetachedSheets==0;
        std::printf("[FlatSurfaceAudit] rays=%llu detail_stddev=%.4f vox near_zero_share=%.1f%% threshold=50.0%% status=%s\n",
            static_cast<unsigned long long>(TotalHits),std::sqrt(Variance),NearZeroShare*100.0,
            bDetailPass?"PASS":"FAIL");
        std::printf("[GridStairAudit] rays=%llu axis_locked_normals=%.1f%% threshold=60.0%% floor_tagged_hits=%llu floor_axis_locked=%.1f%% family_transition_pairs=%llu transition_rms=%.4f vox stable_rms=%.4f vox surface_detail_rms=%.4f vox relative_join_step=%.3f threshold=%.3f legacy_join_ratio=%.3f status=%s\n",
            static_cast<unsigned long long>(TotalHits),AxisShare*100.0,
            static_cast<unsigned long long>(FloorTaggedHits),FloorAxisShare*100.0,
            static_cast<unsigned long long>(FamilyTransitionPairs),TransitionRMS,StableRMS,
            SurfaceDetailRMS,RelativeJoinStep,GridStairThreshold,TransitionRatio,bGridPass?"PASS":"FAIL");
        std::printf("[ThinSheetAudit] rays=%llu isolated_sheets=%llu thickness_lt=1.000 vox air_gap_ge=3.000 vox thinnest=%.3f vox status=%s\n",
            static_cast<unsigned long long>(TotalHits),static_cast<unsigned long long>(ThinDetachedSheets),
            ThinDetachedSheets?ThinnestDetachedSheet:0.0,bThinSheetPass?"PASS":"FAIL");
        for(size_t I=0;I<FamilyStats.size();++I)
        {
            const char* Name=I==0?"strata":I==1?"cave":"solid_overlay";
            const FFamilyStats& S=FamilyStats[I];
            const double M=S.Hits?S.Sum/S.Hits:0.0;
            const double V=S.Hits?std::sqrt(std::max(0.0,S.SumSq/S.Hits-M*M)):0.0;
            std::printf("[SurfaceFamilyAudit] family=%s rays=%llu detail_mean=%.4f detail_stddev=%.4f near_zero=%llu axis_locked=%llu\n",
                Name,static_cast<unsigned long long>(S.Hits),M,V,
                static_cast<unsigned long long>(S.NearZero),static_cast<unsigned long long>(S.AxisLocked));
        }
        const bool bFloorProfilePass=Stats.MaxUnmarkedFloorSlopeDelta<=0.75
            &&Stats.MaxFloorCurveDeviation<=2.0;
        return bDetailPass&&bGridPass&&bThinSheetPass&&bFloorProfilePass;
    }

private:
    float ApplyEditBrushes(float Density,double X,double Y,double Z) const
    {
        if(EditBrushes.empty()||X<EditBoundsMin.X||X>EditBoundsMax.X
            ||Y<EditBoundsMin.Y||Y>EditBoundsMax.Y||Z<EditBoundsMin.Z||Z>EditBoundsMax.Z)
            return Density;
        const double P[3]={X,Y,Z};
        for(const FEditBrush& Brush:EditBrushes)
        {
            double Distance=0.0;
            if(Brush.Shape==FEditBrush::EShape::Sphere)
            {
                const double DX=X-Brush.A[0],DY=Y-Brush.A[1],DZ=Z-Brush.A[2];
                Distance=DeterministicSqrt(DX*DX+DY*DY+DZ*DZ)-Brush.Radius;
            }
            else if(Brush.Shape==FEditBrush::EShape::Capsule)
            {
                double AB[3]{},AP[3]{},LengthSq=0.0,Dot=0.0;
                for(int32 Axis=0;Axis<3;++Axis)
                {AB[Axis]=Brush.B[Axis]-Brush.A[Axis];AP[Axis]=P[Axis]-Brush.A[Axis];LengthSq+=AB[Axis]*AB[Axis];Dot+=AB[Axis]*AP[Axis];}
                const double T=LengthSq>1.0e-12?std::clamp(Dot/LengthSq,0.0,1.0):0.0;
                double D2=0.0;
                for(int32 Axis=0;Axis<3;++Axis)
                {const double D=AP[Axis]-T*AB[Axis];D2+=D*D;}
                Distance=DeterministicSqrt(D2)-Brush.Radius;
            }
            else
            {
                double OutsideSq=0.0,MaxQ=-std::numeric_limits<double>::infinity();
                for(int32 Axis=0;Axis<3;++Axis)
                {
                    const double Q=std::abs(P[Axis]-Brush.A[Axis])-Brush.HalfExtent[Axis];
                    const double Outside=std::max(Q,0.0);
                    OutsideSq+=Outside*Outside;MaxQ=std::max(MaxQ,Q);
                }
                Distance=DeterministicSqrt(OutsideSq)+std::min(MaxQ,0.0)-Brush.Rounding;
            }
            if(Brush.Operation==FEditBrush::EOperation::Carve)
                Density=std::max(Density,static_cast<float>(-Distance));
            else
                Density=std::min(Density,static_cast<float>(Distance));
        }
        return Density;
    }

    struct FAxisBins { int32 MinA=0,MinB=0,CountA=0,CountB=0; std::vector<std::vector<uint32>> Items; };
    struct FLayerField
    {
        int32 Lo=0,Hi=-1;
        int32 StrateIndex=INDEX_NONE;
        FStrateGenerationParams Params;
    };
    enum class ESurfaceFamily : uint8 { Strata, Cave, SolidOverlay };

    static bool NearInteger(float V,int32& Out)
    {
        const float R=std::round(V);
        if(std::abs(V-R)>0.0001f||R<static_cast<float>(MIN_int32)||R>static_cast<float>(MAX_int32))return false;
        Out=static_cast<int32>(R);return true;
    }

    static double DeterministicSqrt(double Value)
    {
        if(!(Value>0.0))return 0.0;
        return _mm_cvtsd_f64(_mm_sqrt_sd(_mm_setzero_pd(),_mm_set_sd(Value)));
    }

    bool Contains(int32 X,int32 Y,int32 Z) const
    {return X>=Min[0]&&X<=Max[0]&&Y>=Min[1]&&Y<=Max[1]&&Z>=Min[2]&&Z<=Max[2];}

    static int32 FloorDiv(int32 A,int32 B)
    {int32 Q=A/B,R=A%B;return Q-((R!=0&&R<0)?1:0);}

    static int64 RoundDiv(int64 N,int64 D)
    {return N>=0?(N+D/2)/D:-((-N+D/2)/D);}
    static float SmoothMinDistance(float A,float B,float Width)
    {
        if(Width<=0.0f)return std::min(A,B);
        const float H=std::max(Width-std::abs(A-B),0.0f)/Width;
        return std::min(A,B)-(H*H*H)*(Width/6.0f);
    }
    static float SmoothMaxDistance(float A,float B,float Width)
    {return -SmoothMinDistance(-A,-B,Width);}
    static float SmoothUnit(float V)
    {
        V=std::clamp(V,0.0f,1.0f);
        return V*V*(3.0f-2.0f*V);
    }
    static float DetailBandFade(float SignedDistance)
    {
        // Keep the larger candidate mask conservative around smooth family unions;
        // visible relief still fades to zero beyond its 4.125-voxel support radius.
        constexpr float Edge=DetailSupremum+0.125f;
        constexpr float Width=1.5f;
        return SmoothUnit((Edge-std::abs(SignedDistance))/Width);
    }
    static constexpr float DetailBandInnerRadius=1.5f;
    static float DetailBandFade(float SignedDistance,float Edge)
    {
        const float Distance=std::abs(SignedDistance);
        if(Distance<=DetailBandInnerRadius)return 1.0f;
        const float Width=Edge-DetailBandInnerRadius;
        return Width>0.0f?SmoothUnit((Edge-Distance)/Width):0.0f;
    }
    static float LocalDetailBandEdge(float MacroValue,float SurfaceSDF,
        float CaveInfluence,float FloorInfluence,const FLayerField* Field)
    {
        float OwnerBound=0.0f;
        if(Field)
        {
            const FStrateGenerationParams& Params=Field->Params;
            const float Depth=Params.SurfaceRoughness*2.0f;
            if(Depth>0.0f&&std::abs(SurfaceSDF)<Depth&&Params.RoughnessFrequency>0.0f)
            {
                const float Fade=1.0f-std::abs(SurfaceSDF)/Depth;
                // Perlin is bounded by 2; the squared ridged fold by 2.25. This
                // uses the largest type bound and remains valid under domain warp.
                OwnerBound=3.15f*VOXEL_NOISE_SCALE*Params.SurfaceRoughness*Fade*Fade;
            }
        }
        const float ReliefBound=std::abs(MacroValue)
            +(0.80f+0.25f*MacroValue)*(0.85f+0.18f);
        const float Contrast=1.35f-0.80f*std::clamp(CaveInfluence,0.0f,1.0f)
            +1.85f*std::clamp(FloorInfluence,0.0f,1.0f);
        const float LocalSupremum=std::min(3.9f,(ReliefBound+0.5f*OwnerBound)*Contrast);
        return std::clamp(LocalSupremum+0.125f,DetailBandInnerRadius,DetailSupremum+0.125f);
    }

    float DetailBandFadeAt(double X,double Y,double Z,float SignedDistance,float SurfaceSDF,
        float CaveInfluence,float FloorInfluence) const
    {
        if(DetailVariant!=EDetailVariant::Coherent
            ||std::abs(SignedDistance)<=DetailBandInnerRadius)
            return DetailBandFade(SignedDistance);
        const float MacroValue=FStandaloneResearchTunnels::ValueNoise(
            static_cast<float>(X),static_cast<float>(Y),static_cast<float>(Z),32,
            static_cast<uint32>(Owner->Seed)^0xB5297A4Du);
        const FLayerField* Field=FindLayerField(Z);
        return DetailBandFade(SignedDistance,LocalDetailBandEdge(MacroValue,SurfaceSDF,
            CaveInfluence,FloorInfluence,Field));
    }

    void BuildLayerBounds()
    {
        LayerRanges.clear();LayerFields.clear();
        const auto& Layout=Owner->Manager->GetLayout();
        int32 Top=MIN_int32,Bottom=MAX_int32;
        for(const FStrateSlot& Slot:Layout)
        {
            Top=std::max(Top,Slot.TopChunkZ);Bottom=std::min(Bottom,Slot.BottomChunkZ);
            if(!Slot.Definition)continue;
            const int32 MidChunkZ=(Slot.TopChunkZ+Slot.BottomChunkZ)/2;
            FLayerField Field;
            Field.Lo=Slot.BottomChunkZ*CHUNK_SIZE;
            Field.Hi=(Slot.TopChunkZ+1)*CHUNK_SIZE;
            Field.StrateIndex=Slot.StrateIndex;
            Field.Params=Owner->Manager->GetGenerationParams(FIntVector(0,0,MidChunkZ));
            LayerFields.push_back(std::move(Field));
        }
        if(Top>=Bottom)
        {
            LayerRanges.push_back({Bottom*CHUNK_SIZE,(Top+1)*CHUNK_SIZE});
        }
    }

    static std::array<double,3> V3(const FVector& V){return {V.X,V.Y,V.Z};}
    static double Dot(const std::array<double,3>& A,const std::array<double,3>& B)
    {return A[0]*B[0]+A[1]*B[1]+A[2]*B[2];}
    static std::array<double,3> Sub(const FVector& A,const FVector& B)
    {return {A.X-B.X,A.Y-B.Y,A.Z-B.Z};}
    static std::array<double,3> Normalize(std::array<double,3> V)
    {
        const double L=std::sqrt(std::max(1.0e-24,Dot(V,V)));
        for(double& C:V)C/=L;return V;
    }
    static std::array<double,3> Cross(const std::array<double,3>& A,const std::array<double,3>& B)
    {return {A[1]*B[2]-A[2]*B[1],A[2]*B[0]-A[0]*B[2],A[0]*B[1]-A[1]*B[0]};}

    void AddEllipsoid(const FVector& Center,const std::array<double,3>& Radii,
                      const std::array<double,3>& Axis0={1,0,0},
                      const std::array<double,3>& Axis1={0,1,0},
                      const std::array<double,3>& Axis2={0,0,1},bool bCarve=true)
    {
        FPrimitive P;const auto C=V3(Center);
        for(int32 A=0;A<3;++A)P.CenterSub[A]=static_cast<int32>(std::llround(C[A]*SubPerVoxel));
        const std::array<std::array<double,3>,3> Axes={Axis0,Axis1,Axis2};
        int32 RSub[3];
        for(int32 K=0;K<3;++K)RSub[K]=std::max(16,static_cast<int32>(std::llround(std::max(1.0,Radii[K])*SubPerVoxel)));
        for(int32 I=0;I<3;++I)for(int32 J=0;J<3;++J)
        {
            int64 Sum=0;
            for(int32 K=0;K<3;++K)
            {
                const int64 RI=static_cast<int64>(std::llround(Axes[K][I]*32768.0));
                const int64 RJ=static_cast<int64>(std::llround(Axes[K][J]*32768.0));
                const int64 Den=static_cast<int64>(32768)*32768*RSub[K]*RSub[K];
                Sum+=RoundDiv(RI*RJ*QScale,Den);
            }
            P.Q[I][J]=Sum;
        }
        for(int32 A=0;A<3;++A)
        {
            double E=0.0;for(int32 K=0;K<3;++K)E+=std::abs(Axes[K][A])*Radii[K];
            P.Min[A]=static_cast<int32>(std::floor(C[A]-E))-1;
            P.Max[A]=static_cast<int32>(std::ceil(C[A]+E))+1;
        }
        P.MinRadius=static_cast<float>(std::min({Radii[0],Radii[1],Radii[2]}));
        P.BlendRadius=ActiveBlendRadius;
        P.bCarve=bCarve;Primitives.push_back(P);
    }

    void AddRoundedBox(const FVector& Center,const std::array<double,3>& HalfExtent,
                       double Rounding,bool bCarve=true)
    {
        const std::array<double,3> Outer={HalfExtent[0]+Rounding,HalfExtent[1]+Rounding,HalfExtent[2]+Rounding};
        AddEllipsoid(Center,Outer,{1,0,0},{0,1,0},{0,0,1},bCarve);
        FPrimitive& P=Primitives.back();P.bRoundedBox=true;P.Rounding=Rounding;
        for(int32 A=0;A<3;++A)P.RoundedHalfExtent[A]=HalfExtent[A];
    }

    static void ApplyFloorClip(FPrimitive& P,double FloorCenterZ,double SlopeX=0.0,double SlopeY=0.0,
                               double FloorBlendRadius=0.0)
    {
        P.bHasFloorClip=true;P.FloorCenterZ=FloorCenterZ-0.5;
        P.FloorSlopeX=SlopeX;P.FloorSlopeY=SlopeY;P.FloorBlendRadius=FloorBlendRadius;
    }

    void AddCapsule(const FVector& A,const FVector& B,float Radius,bool bCarve=true,bool bCritical=false,
                    bool bFloorClip=false,double FloorA=0.0,double FloorB=0.0)
    {
        AddLegacyEllipsoidCapsule(A,B,Radius,bCarve,bCritical,bFloorClip,FloorA,FloorB);
    }

    static FVector ApproximateUnwarp(const FVector& Target,const FStrateGenerationParams& Params,
                                     uint32 Seed)
    {
        FVector WorldPoint=Target;
        const float ZScale=Params.VerticalScale>0.0f?Params.VerticalScale:1.0f;
        if(ZScale!=1.0f)WorldPoint.Z*=ZScale;
        if(Params.CaveWarpStrength<=0.0f)return WorldPoint;
        for(int32 Iteration=0;Iteration<8;++Iteration)
        {
            const FVector Mapped=VoxelCaveMorphology::ApplyCaveWarp(WorldPoint,Params,Seed);
            const FVector Error=Target-Mapped;
            if(!Error.ContainsNaN()&&VoxelMath::IsFinite(Error.X)
                &&VoxelMath::IsFinite(Error.Y)&&VoxelMath::IsFinite(Error.Z))
            {
                WorldPoint.X+=Error.X;WorldPoint.Y+=Error.Y;WorldPoint.Z+=Error.Z*ZScale;
            }
        }
        return WorldPoint;
    }

    void AddLegacyEllipsoidCapsule(const FVector& A,const FVector& B,float Radius,bool bCarve,
                                   bool bCritical,bool bFloorClip,double FloorA,double FloorB)
    {
        const auto D=Sub(B,A);const double Length=std::sqrt(Dot(D,D));
        const double PreviousBlendRadius=ActiveBlendRadius;ActiveBlendRadius=0.0;
        if(!(Length>1.0e-4))
        {
            const size_t Before=Primitives.size();
            AddEllipsoid(A,{Radius,Radius,Radius},{1,0,0},{0,1,0},{0,0,1},bCarve);
            if(bFloorClip)ApplyFloorClip(Primitives.back(),std::min(FloorA,FloorB),0.0,0.0,
                0.0);
            if(bCritical)for(size_t I=Before;I<Primitives.size();++I)Primitives[I].bCritical=true;
            ActiveBlendRadius=PreviousBlendRadius;return;
        }
        const auto U=Normalize(D);
        const std::array<double,3> Ref=std::abs(U[2])<0.90?std::array<double,3>{0,0,1}:std::array<double,3>{0,1,0};
        const auto V=Normalize(Cross(U,Ref));const auto W=Normalize(Cross(U,V));
        const double SegmentStep=std::max(8.0,static_cast<double>(Radius)*2.0);
        const int32 Segments=std::max(1,static_cast<int32>(std::ceil(Length/SegmentStep)));
        for(int32 I=0;I<Segments;++I)
        {
            const double T0=static_cast<double>(I)/Segments,T1=static_cast<double>(I+1)/Segments;
            const FVector P0(A.X+D[0]*T0,A.Y+D[1]*T0,A.Z+D[2]*T0);
            const FVector P1(A.X+D[0]*T1,A.Y+D[1]*T1,A.Z+D[2]*T1);
            const FVector C((P0.X+P1.X)*0.5,(P0.Y+P1.Y)*0.5,(P0.Z+P1.Z)*0.5);
            const double Half=Length*(T1-T0)*0.5+Radius;
            const size_t Before=Primitives.size();AddEllipsoid(C,{Half,Radius,Radius},U,V,W,bCarve);
            if(bFloorClip)
            {
                const double Mid=(T0+T1)*0.5;
                const double Floor=(1.0-Mid)*FloorA+Mid*FloorB;
                const double DX=P1.X-P0.X,DY=P1.Y-P0.Y,HorizontalSq=DX*DX+DY*DY;
                const double FloorStep=(FloorB-FloorA)*(T1-T0);
                if(HorizontalSq>1.0e-6)
                    ApplyFloorClip(Primitives.back(),Floor,FloorStep*DX/HorizontalSq,
                        FloorStep*DY/HorizontalSq,0.0);
                else ApplyFloorClip(Primitives.back(),std::min(FloorA,FloorB),0.0,0.0,
                    0.0);
            }
            if(bCritical)for(size_t J=Before;J<Primitives.size();++J)Primitives[J].bCritical=true;
        }
        ActiveBlendRadius=PreviousBlendRadius;
    }

    void AddTaperedCapsule(const FVector& A,const FVector& B,float RadiusA,float RadiusB,
                           bool bCarve=true,bool bCritical=false,bool bFloorClip=false,
                           double FloorA=0.0,double FloorB=0.0,
                           double ExtendStart=0.0,double ExtendEnd=0.0)
    {
        RadiusA=std::max(1.0f,RadiusA);RadiusB=std::max(1.0f,RadiusB);
        const auto D=Sub(B,A);const double Length=std::sqrt(Dot(D,D));
        if(!(Length>1.0e-4))
        {
            AddEllipsoid(A,{RadiusA,RadiusA,RadiusA},{1,0,0},{0,1,0},{0,0,1},bCarve);
            Primitives.back().bCritical=bCritical;
            if(bFloorClip)ApplyFloorClip(Primitives.back(),std::min(FloorA,FloorB),0.0,0.0,
                ActiveFloorBlendRadius);
            return;
        }
        ExtendStart=std::clamp(ExtendStart,0.0,Length*0.4);
        ExtendEnd=std::clamp(ExtendEnd,0.0,Length*0.4);
        const double RadiusSlope=(static_cast<double>(RadiusB)-RadiusA)/Length;
        const float ExtendedRadiusA=static_cast<float>(std::max(1.0,RadiusA-RadiusSlope*ExtendStart));
        const float ExtendedRadiusB=static_cast<float>(std::max(1.0,RadiusB+RadiusSlope*ExtendEnd));
        const double Start[3]={A.X-D[0]*(ExtendStart/Length),A.Y-D[1]*(ExtendStart/Length),A.Z-D[2]*(ExtendStart/Length)};
        const double End[3]={B.X+D[0]*(ExtendEnd/Length),B.Y+D[1]*(ExtendEnd/Length),B.Z+D[2]*(ExtendEnd/Length)};
        FPrimitive P;
        P.bCapsule=true;P.bCarve=bCarve;P.bCritical=bCritical;
        P.CapsuleStartRadius=ExtendedRadiusA;P.CapsuleEndRadius=ExtendedRadiusB;
        P.MinRadius=std::min(ExtendedRadiusA,ExtendedRadiusB);P.BlendRadius=ActiveBlendRadius;
        const double R=std::max(ExtendedRadiusA,ExtendedRadiusB);
        double AxisLengthSquared=0.0;
        for(int32 Axis=0;Axis<3;++Axis)
        {
            P.CapsuleStartSub[Axis]=static_cast<int32>(std::llround(Start[Axis]*SubPerVoxel));
            P.CapsuleEndSub[Axis]=static_cast<int32>(std::llround(End[Axis]*SubPerVoxel));
            P.CapsuleStartWorld[Axis]=P.CapsuleStartSub[Axis]/static_cast<double>(SubPerVoxel);
            const double CapsuleEndWorld=P.CapsuleEndSub[Axis]/static_cast<double>(SubPerVoxel);
            P.CapsuleAxisDelta[Axis]=CapsuleEndWorld-P.CapsuleStartWorld[Axis];
            AxisLengthSquared+=P.CapsuleAxisDelta[Axis]*P.CapsuleAxisDelta[Axis];
            P.CenterSub[Axis]=static_cast<int32>(std::llround((Start[Axis]+End[Axis])*0.5*SubPerVoxel));
            P.Min[Axis]=static_cast<int32>(std::floor(std::min(Start[Axis],End[Axis])-R))-1;
            P.Max[Axis]=static_cast<int32>(std::ceil(std::max(Start[Axis],End[Axis])+R))+1;
        }
        P.CapsuleAxisLengthSquared=AxisLengthSquared;
        P.CapsuleAxisLength=std::sqrt(std::max(0.0,AxisLengthSquared));
        if(bFloorClip)
        {
            const double FloorSlope=(FloorB-FloorA)/Length;
            const double ExtendedFloorA=FloorA-FloorSlope*ExtendStart;
            const double ExtendedFloorB=FloorB+FloorSlope*ExtendEnd;
            const double ExtendedD[3]={End[0]-Start[0],End[1]-Start[1],End[2]-Start[2]};
            const double HorizontalSq=ExtendedD[0]*ExtendedD[0]+ExtendedD[1]*ExtendedD[1];
            if(HorizontalSq>1.0e-6)
                ApplyFloorClip(P,(ExtendedFloorA+ExtendedFloorB)*0.5,
                    (ExtendedFloorB-ExtendedFloorA)*ExtendedD[0]/HorizontalSq,
                    (ExtendedFloorB-ExtendedFloorA)*ExtendedD[1]/HorizontalSq,ActiveFloorBlendRadius);
            else ApplyFloorClip(P,std::min(ExtendedFloorA,ExtendedFloorB),0.0,0.0,ActiveFloorBlendRadius);
        }
        Primitives.push_back(P);
    }

    static double ContinuousJoinExtension(float RadiusA,float RadiusB,double SegmentLength)
    {
        return std::min(std::max(0.0,SegmentLength*0.35),
            JoinBlendRadius+0.5*static_cast<double>(std::max(RadiusA,RadiusB)));
    }

    static float CapsuleRadiusAt(const FPrimitive& P,const FVector& Point)
    {
        if(!P.bCapsule)return P.MinRadius;
        const double Start[3]={P.CapsuleStartSub[0]/static_cast<double>(SubPerVoxel),
            P.CapsuleStartSub[1]/static_cast<double>(SubPerVoxel),
            P.CapsuleStartSub[2]/static_cast<double>(SubPerVoxel)};
        const double End[3]={P.CapsuleEndSub[0]/static_cast<double>(SubPerVoxel),
            P.CapsuleEndSub[1]/static_cast<double>(SubPerVoxel),
            P.CapsuleEndSub[2]/static_cast<double>(SubPerVoxel)};
        const double V[3]={Point.X-Start[0],Point.Y-Start[1],Point.Z-Start[2]};
        const double D[3]={End[0]-Start[0],End[1]-Start[1],End[2]-Start[2]};
        const double D2=D[0]*D[0]+D[1]*D[1]+D[2]*D[2];
        const double T=D2>1.0e-12?(V[0]*D[0]+V[1]*D[1]+V[2]*D[2])/D2:0.0;
        return CapsuleRadiusAtT(P,T,std::sqrt(std::max(0.0,D2)));
    }

    static void FitCapsuleRadiiAtEndpoints(FPrimitive& P,const FVector& PointA,float RadiusA,
                                            const FVector& PointB,float RadiusB)
    {
        const double Start[3]={P.CapsuleStartSub[0]/static_cast<double>(SubPerVoxel),
            P.CapsuleStartSub[1]/static_cast<double>(SubPerVoxel),
            P.CapsuleStartSub[2]/static_cast<double>(SubPerVoxel)};
        const double End[3]={P.CapsuleEndSub[0]/static_cast<double>(SubPerVoxel),
            P.CapsuleEndSub[1]/static_cast<double>(SubPerVoxel),
            P.CapsuleEndSub[2]/static_cast<double>(SubPerVoxel)};
        const double D[3]={End[0]-Start[0],End[1]-Start[1],End[2]-Start[2]};
        const double D2=D[0]*D[0]+D[1]*D[1]+D[2]*D[2];
        if(D2<=1.0e-12)return;
        const auto Parameter=[&](const FVector& Point)
        {
            const double V[3]={Point.X-Start[0],Point.Y-Start[1],Point.Z-Start[2]};
            return (V[0]*D[0]+V[1]*D[1]+V[2]*D[2])/D2;
        };
        const double TA=Parameter(PointA),TB=Parameter(PointB);
        if(std::abs(TB-TA)<=1.0e-9)return;
        const double Slope=(RadiusB-RadiusA)/(TB-TA);
        P.CapsuleStartRadius=static_cast<float>(RadiusA-Slope*TA);
        P.CapsuleEndRadius=static_cast<float>(P.CapsuleStartRadius+Slope);
        P.MinRadius=std::min(P.CapsuleStartRadius,P.CapsuleEndRadius);
    }

    static void SetCapsuleRadiusProfile(FPrimitive& P,const FVector& PointA,float RadiusA,
        double SlopeA,const FVector& PointB,float RadiusB,double SlopeB)
    {
        const double Start[3]={P.CapsuleStartSub[0]/static_cast<double>(SubPerVoxel),
            P.CapsuleStartSub[1]/static_cast<double>(SubPerVoxel),
            P.CapsuleStartSub[2]/static_cast<double>(SubPerVoxel)};
        const double End[3]={P.CapsuleEndSub[0]/static_cast<double>(SubPerVoxel),
            P.CapsuleEndSub[1]/static_cast<double>(SubPerVoxel),
            P.CapsuleEndSub[2]/static_cast<double>(SubPerVoxel)};
        const double D[3]={End[0]-Start[0],End[1]-Start[1],End[2]-Start[2]};
        const double D2=D[0]*D[0]+D[1]*D[1]+D[2]*D[2];
        if(D2<=1.0e-12)return;
        const auto Parameter=[&](const FVector& Point)
        {
            const double V[3]={Point.X-Start[0],Point.Y-Start[1],Point.Z-Start[2]};
            return std::clamp((V[0]*D[0]+V[1]*D[1]+V[2]*D[2])/D2,0.0,1.0);
        };
        P.CapsuleProfileStartT=Parameter(PointA);
        P.CapsuleProfileEndT=Parameter(PointB);
        P.CapsuleProfileStartRadius=RadiusA;
        P.CapsuleProfileEndRadius=RadiusB;
        P.CapsuleRadiusSlopeStart=SlopeA;
        P.CapsuleRadiusSlopeEnd=SlopeB;
        P.bSmoothCapsuleRadius=P.CapsuleProfileEndT-P.CapsuleProfileStartT>1.0e-6;
    }

    static float SmoothedControlRadius(const TArray<float>& Radii,int32 Index,float Fallback,
                                        float Minimum)
    {
        const int32 Count=Radii.Num();
        const float Current=std::max(Minimum,Radii.IsValidIndex(Index)?Radii[Index]:Fallback);
        if(Index<=0||Index+1>=Count)return Current;
        const auto At=[&](int32 I)
        {return std::max(Minimum,Radii[std::clamp(I,0,Count-1)]);};
        return (At(Index-4)+8.0f*At(Index-3)+28.0f*At(Index-2)+56.0f*At(Index-1)
            +70.0f*Current+56.0f*At(Index+1)+28.0f*At(Index+2)+8.0f*At(Index+3)
            +At(Index+4))*(1.0f/256.0f);
    }

    static double ControlRadiusSlopeAt(const TArray<FVector>& Points,const TArray<float>& Radii,
        int32 Index,float Fallback,float Minimum)
    {
        if(Points.Num()<2||Index<0||Index>=Points.Num())return 0.0;
        const auto RadiusAt=[&](int32 I)
        {return static_cast<double>(SmoothedControlRadius(Radii,I,Fallback,Minimum));};
        const auto LengthAt=[&](int32 I)
        {
            const FVector& A=Points[I];const FVector& B=Points[I+1];
            const double DX=B.X-A.X,DY=B.Y-A.Y,DZ=B.Z-A.Z;
            return std::sqrt(DX*DX+DY*DY+DZ*DZ);
        };
        if(Index==0)
        {
            const double Length=LengthAt(0);
            return Length>1.0e-9?(RadiusAt(1)-RadiusAt(0))/Length:0.0;
        }
        if(Index+1==Points.Num())
        {
            const double Length=LengthAt(Index-1);
            return Length>1.0e-9?(RadiusAt(Index)-RadiusAt(Index-1))/Length:0.0;
        }
        const double PrevLength=LengthAt(Index-1),NextLength=LengthAt(Index);
        if(PrevLength<=1.0e-9||NextLength<=1.0e-9)return 0.0;
        const double PrevSlope=(RadiusAt(Index)-RadiusAt(Index-1))/PrevLength;
        const double NextSlope=(RadiusAt(Index+1)-RadiusAt(Index))/NextLength;
        if(PrevSlope*NextSlope<=0.0)return 0.0;
        const double W1=2.0*NextLength+PrevLength,W2=NextLength+2.0*PrevLength;
        const double Harmonic=(W1+W2)/(W1/PrevSlope+W2/NextSlope);
        const double Limit=3.0*std::min(std::abs(PrevSlope),std::abs(NextSlope));
        return std::copysign(std::min(std::abs(Harmonic),Limit),Harmonic);
    }

    static double EvaluateFloorProfile(double A,double B,double TangentA,double TangentB,
        double Length,double StartBlend,double EndBlend,double S)
    {
        const double Secant=(B-A)/Length;
        if(S<0.0)return A+TangentA*S;
        if(S>Length)return B+TangentB*(S-Length);
        double Result=std::lerp(A,B,S/Length);
        const auto Correction=[](double DeltaSlope,double Span,double T)
        {
            if(!(Span>0.0)||T<0.0||T>=1.0)return 0.0;
            return DeltaSlope*Span*(T-2.0*T*T+T*T*T);
        };
        if(S<StartBlend)
            Result+=Correction(TangentA-Secant,StartBlend,S/StartBlend);
        if(Length-S<EndBlend)
            Result-=Correction(TangentB-Secant,EndBlend,(Length-S)/EndBlend);
        return Result;
    }

    static double FloorSurfaceRelief(double X,double Y,uint32 Seed)
    {
        // A shared XY field keeps independently authored floor clips from ending on
        // separate flat planes. Its peak displacement stays under 1.2 voxels.
        const float Broad=FStandaloneResearchTunnels::ValueNoise(
            static_cast<float>(X),static_cast<float>(Y),0.0f,32,Seed^0x7F4A7C15u);
        const float Mid=FStandaloneResearchTunnels::ValueNoise(
            static_cast<float>(X),static_cast<float>(Y),0.0f,16,Seed^0x94D049BBu);
        return 0.80*Broad+0.35*Mid;
    }

    static std::vector<double> BuildFloorTangents(const TArray<FVector>& Points,
        const std::vector<float>& Floor,const std::vector<uint8_t>& MarkedSegments)
    {
        const size_t Count=Floor.size();std::vector<double> Tangent(Count,0.0);
        if(Count<2||Count!=static_cast<size_t>(Points.Num())||MarkedSegments.size()+1!=Count)return Tangent;
        std::vector<double> Length(Count-1,0.0),Slope(Count-1,0.0);
        for(size_t I=0;I+1<Count;++I)
        {
            const FVector& A=Points[static_cast<int32>(I)];const FVector& B=Points[static_cast<int32>(I+1)];
            Length[I]=std::hypot(B.X-A.X,B.Y-A.Y);
            if(Length[I]>1.0e-6)Slope[I]=(Floor[I+1]-Floor[I])/Length[I];
        }
        for(size_t I=0;I<Count;++I)
        {
            const bool bHasPrev=I>0&&Length[I-1]>1.0e-6;
            const bool bHasNext=I+1<Count&&Length[I]>1.0e-6;
            if(I==0){Tangent[I]=MarkedSegments[0]?0.0:Slope[0];continue;}
            if(I+1==Count){Tangent[I]=MarkedSegments[I-1]?0.0:Slope[I-1];continue;}
            const bool bPrevMarked=MarkedSegments[I-1]!=0,bNextMarked=MarkedSegments[I]!=0;
            if(bPrevMarked||bNextMarked)
            {
                if(!bPrevMarked&&bHasPrev)Tangent[I]=Slope[I-1];
                else if(!bNextMarked&&bHasNext)Tangent[I]=Slope[I];
                continue;
            }
            if(!bHasPrev||!bHasNext||Slope[I-1]*Slope[I]<=0.0)continue;
            const double W0=2.0*Length[I]+Length[I-1],W1=Length[I]+2.0*Length[I-1];
            const double Denominator=W0/Slope[I-1]+W1/Slope[I];
            if(std::abs(Denominator)>1.0e-12)
                Tangent[I]=(W0+W1)/Denominator;
        }
        return Tangent;
    }

    static void SetCurvedFloorClip(FPrimitive& P,const FVector& A,const FVector& B,
        double FloorA,double FloorB,double TangentA,double TangentB)
    {
        if(!P.bHasFloorClip)return;
        const double DX=B.X-A.X,DY=B.Y-A.Y,Length=std::hypot(DX,DY);
        if(Length<=1.0e-6)return;
        P.bCurvedFloorClip=true;P.FloorCurveOriginX=A.X;P.FloorCurveOriginY=A.Y;
        P.FloorCurveDeltaX=DX;P.FloorCurveDeltaY=DY;P.FloorCurveLength=Length;
        P.FloorCurveStartZ=FloorA;P.FloorCurveEndZ=FloorB;
        P.FloorCurveTangentStart=TangentA;P.FloorCurveTangentEnd=TangentB;
        P.FloorCurveStartBlend=std::min(8.0,Length*0.35);
        P.FloorCurveEndBlend=std::min(8.0,Length*0.35);
        for(int32 I=1;I<64;++I)
        {
            const double T=I/64.0;
            const double Curve=EvaluateFloorProfile(FloorA,FloorB,TangentA,TangentB,
                Length,P.FloorCurveStartBlend,P.FloorCurveEndBlend,Length*T);
            const double Linear=std::lerp(FloorA,FloorB,T);
            P.FloorCurveMaxDeviation=std::max(P.FloorCurveMaxDeviation,std::abs(Curve-Linear));
        }
    }

    double FloorAt(const FPrimitive& P,double X,double Y) const
    {
        double Base=0.0;
        if(P.bCurvedFloorClip)
        {
            const double T=((X-P.FloorCurveOriginX)*P.FloorCurveDeltaX
                +(Y-P.FloorCurveOriginY)*P.FloorCurveDeltaY)
                /(P.FloorCurveLength*P.FloorCurveLength);
            Base=EvaluateFloorProfile(P.FloorCurveStartZ,P.FloorCurveEndZ,
                P.FloorCurveTangentStart,P.FloorCurveTangentEnd,P.FloorCurveLength,
                P.FloorCurveStartBlend,P.FloorCurveEndBlend,T*P.FloorCurveLength);
        }
        else
            Base=P.FloorCenterZ+P.FloorSlopeX*(X-P.CenterSub[0]/16.0)
                +P.FloorSlopeY*(Y-P.CenterSub[1]/16.0);
        return Base;
    }

    double FloorSlopeAtEndpoint(const FPrimitive& P,bool bEnd) const
    {
        const double DX=P.bCurvedFloorClip?P.FloorCurveDeltaX
            :P.CapsuleEndSub[0]-P.CapsuleStartSub[0];
        const double DY=P.bCurvedFloorClip?P.FloorCurveDeltaY
            :P.CapsuleEndSub[1]-P.CapsuleStartSub[1];
        const double Length=P.bCurvedFloorClip?P.FloorCurveLength
            :std::hypot(DX,DY)/SubPerVoxel;
        if(Length<=1.0e-9)return 0.0;
        const double OriginX=P.bCurvedFloorClip?P.FloorCurveOriginX
            :P.CapsuleStartSub[0]/static_cast<double>(SubPerVoxel);
        const double OriginY=P.bCurvedFloorClip?P.FloorCurveOriginY
            :P.CapsuleStartSub[1]/static_cast<double>(SubPerVoxel);
        const double X=OriginX+(bEnd?DX/(P.bCurvedFloorClip?1.0:SubPerVoxel):0.0);
        const double Y=OriginY+(bEnd?DY/(P.bCurvedFloorClip?1.0:SubPerVoxel):0.0);
        constexpr double Epsilon=0.125;
        const double UX=DX/(Length*(P.bCurvedFloorClip?1.0:SubPerVoxel));
        const double UY=DY/(Length*(P.bCurvedFloorClip?1.0:SubPerVoxel));
        return (FloorAt(P,X+UX*Epsilon,Y+UY*Epsilon)
            -FloorAt(P,X-UX*Epsilon,Y-UY*Epsilon))/(2.0*Epsilon);
    }

    void BuildOwnerPrimitives()
    {
        const auto& Layout=Owner->Manager->GetLayout();
        const float MinX=static_cast<float>(Min[0]-256),MinY=static_cast<float>(Min[1]-256);
        const float MaxX=static_cast<float>(Max[0]+256),MaxY=static_cast<float>(Max[1]+256);
        for(const FStrateSlot& Slot:Layout)
        {
            if(!Slot.Definition)continue;
            const FIntVector MidChunk(0,0,(Slot.TopChunkZ+Slot.BottomChunkZ)/2);
            const FStrateGenerationParams Params=Owner->Manager->GetGenerationParams(MidChunk);
            ActiveBlendRadius=0.0;
            ActiveFloorBlendRadius=0.0;
            FChunkSDFCache Cache;
            VoxelCaveMorphology::BuildChunkCache(Cache,MinX,MinY,MaxX,MaxY,Params,
                static_cast<uint32>(Owner->Seed),
                Slot.StrateIndex,&Slot.Definition->TerrainOperations,ERoomGraphBuildSite::Unknown);
            for(const FCachedRoom& R:Cache.Rooms)
            {
                ++Stats.RoomCount;
                const size_t RoomPrimitiveStart=Primitives.size();
                const FVector RoomCenter=ApproximateUnwarp(R.Center,Params,
                    static_cast<uint32>(Owner->Seed));
                if(R.ShapeType==2)
                    AddCapsule(ApproximateUnwarp(R.ShapeA,Params,static_cast<uint32>(Owner->Seed)),
                        ApproximateUnwarp(R.ShapeB,Params,static_cast<uint32>(Owner->Seed)),
                        std::max(1.0f,R.ShapeR),true);
                else if(R.ShapeType==1)
                {
                    AddRoundedBox(RoomCenter,{std::abs(R.ShapeA.X),std::abs(R.ShapeA.Y),std::abs(R.ShapeA.Z)},
                        std::max(0.0f,R.ShapeR),true);
                }
                else
                    AddEllipsoid(RoomCenter,{std::abs(R.ShapeA.X),std::abs(R.ShapeA.Y),std::abs(R.ShapeA.Z)},
                        {1,0,0},{0,1,0},{0,0,1},true);
                if(R.FloorCutZ>-FLT_MAX*0.5f)
                {
                    const double Relief=std::max(0.0f,R.FloorReliefStrength);
                    const double FloorZ=ApproximateUnwarp(
                        FVector(R.Center.X,R.Center.Y,R.FloorCutZ-Relief),Params,
                        static_cast<uint32>(Owner->Seed)).Z;
                    for(size_t I=RoomPrimitiveStart;I<Primitives.size();++I)
                        if(Primitives[I].bCarve)ApplyFloorClip(Primitives[I],FloorZ,0.0,0.0,
                            FloorClipBlendRadius);
                }
                // The packed signs still union room intervals exactly. Blend only
                // the analytic cave field where neighboring authored rooms meet;
                // this rounds the room-to-room edge without touching tunnel radii.
                for(size_t I=RoomPrimitiveStart;I<Primitives.size();++I)
                    if(Primitives[I].bCarve)Primitives[I].BlendRadius=RoomUnionBlendRadius;
            }
            // Do not turn automatic room-overlap joins into thin positive ellipsoids.
            // Their 2.5-voxel support plates can sit above the lower room floor and render as
            // disconnected shards. Each room already carries its own clipped floor; the join
            // metadata remains available to the owner generator and is audited here.
            Stats.RoomFloorJoinCount+=static_cast<uint64>(Cache.RoomFloorJoins.Num());
            const double RoomBlendRadius=ActiveBlendRadius;
            ActiveBlendRadius=JoinBlendRadius;
            ActiveFloorBlendRadius=FloorClipBlendRadius;
            for(const FCachedTunnel& T:Cache.Tunnels)
            {
                ++Stats.TunnelCount;
                const bool bHasWorldChain=T.WorldControlPoints.Num()>1
                    &&T.WorldControlRadii.Num()==T.WorldControlPoints.Num();
                const TArray<FVector>& Points=bHasWorldChain?T.WorldControlPoints:T.ControlPoints;
                const TArray<float>& Radii=bHasWorldChain?T.WorldControlRadii:T.ControlRadii;
                const TArray<FTunnelFloorSegmentProfile>& FloorProfiles=bHasWorldChain
                    &&T.WorldFloorProfiles.Num()==Points.Num()-1?T.WorldFloorProfiles:T.FloorProfiles;
                std::vector<float> RawFloor(static_cast<size_t>(Points.Num()));
                std::vector<uint8_t> MarkedFloorSegments(Points.Num()>0?static_cast<size_t>(Points.Num()-1):0,0);
                for(int32 I=0;I<Points.Num();++I)
                    RawFloor[static_cast<size_t>(I)]=static_cast<float>(Points[I].Z
                        -SmoothedControlRadius(Radii,I,T.RadiusA,2.0f));
                if(FloorProfiles.Num()==Points.Num()-1)
                {
                    for(int32 I=0;I+1<Points.Num();++I)
                    {
                        const FTunnelFloorSegmentProfile& Profile=FloorProfiles[I];
                        if(I==0)RawFloor[0]=Profile.StartFloorZ;
                        else RawFloor[static_cast<size_t>(I)]=0.5f*(FloorProfiles[I-1].EndFloorZ+Profile.StartFloorZ);
                        RawFloor[static_cast<size_t>(I+1)]=Profile.EndFloorZ;
                        MarkedFloorSegments[static_cast<size_t>(I)]=Profile.bLedgeTransition?1:0;
                    }
                }
                for(int32 I=0;I<Points.Num();++I)
                {
                    const bool bPreviousUnmarked=I==0
                        ||!MarkedFloorSegments[static_cast<size_t>(I-1)];
                    const bool bNextUnmarked=I+1>=Points.Num()
                        ||!MarkedFloorSegments[static_cast<size_t>(I)];
                    if(bPreviousUnmarked&&bNextUnmarked)
                        RawFloor[static_cast<size_t>(I)]+=static_cast<float>(FloorSurfaceRelief(
                            Points[I].X,Points[I].Y,static_cast<uint32>(Owner->Seed)));
                }
                const std::vector<double> FloorTangents=BuildFloorTangents(
                    Points,RawFloor,MarkedFloorSegments);
                for(int32 Join=1;Join+1<Points.Num();++Join)
                {
                    const FVector A=Points[Join]-Points[Join-1],B=Points[Join+1]-Points[Join];
                    const double AL=std::sqrt(A.X*A.X+A.Y*A.Y+A.Z*A.Z);
                    const double BL=std::sqrt(B.X*B.X+B.Y*B.Y+B.Z*B.Z);
                    if(AL>1.0e-6&&BL>1.0e-6)
                    {
                        const double Dot=(A.X*B.X+A.Y*B.Y+A.Z*B.Z)/(AL*BL);
                        const double Angle=std::acos(std::clamp(Dot,-1.0,1.0))*57.29577951308232;
                        if(Angle>Stats.MaxTunnelBendDegrees)
                        {
                            Stats.MaxTunnelBendDegrees=Angle;Stats.MaxTunnelBendVoxels=Points[Join];
                            const double DX=A.X/AL+B.X/BL,DY=A.Y/AL+B.Y/BL,DZ=A.Z/AL+B.Z/BL;
                            const double DL=std::sqrt(DX*DX+DY*DY+DZ*DZ);
                            if(DL>1.0e-6)Stats.MaxTunnelBendDirection=FVector(DX/DL,DY/DL,DZ/DL);
                        }
                    }
                    const float RPrev=std::max(2.0f,Radii.IsValidIndex(Join-1)?Radii[Join-1]:T.RadiusA);
                    const float RHere=std::max(2.0f,Radii.IsValidIndex(Join)?Radii[Join]:T.RadiusA);
                    const float RNext=std::max(2.0f,Radii.IsValidIndex(Join+1)?Radii[Join+1]:T.RadiusB);
                    Stats.LegacySegmentRadiusStep=std::max(Stats.LegacySegmentRadiusStep,
                        static_cast<double>(std::abs(std::max(RPrev,RHere)-std::max(RHere,RNext))));
                    if(FloorProfiles.IsValidIndex(Join-1)&&FloorProfiles.IsValidIndex(Join))
                    {
                        Stats.MaxCapsuleJoinFloorDelta=std::max(Stats.MaxCapsuleJoinFloorDelta,
                            static_cast<double>(std::abs(FloorProfiles[Join-1].EndFloorZ
                                -FloorProfiles[Join].StartFloorZ)));
                        const FVector PrevD=Points[Join]-Points[Join-1];
                        const FVector NextD=Points[Join+1]-Points[Join];
                        const double PrevXY=std::sqrt(PrevD.X*PrevD.X+PrevD.Y*PrevD.Y);
                        const double NextXY=std::sqrt(NextD.X*NextD.X+NextD.Y*NextD.Y);
                        const double RawPrevSlope=PrevXY>1.0e-6
                            ?(FloorProfiles[Join-1].EndFloorZ-FloorProfiles[Join-1].StartFloorZ)/PrevXY:0.0;
                        const double RawNextSlope=NextXY>1.0e-6
                            ?(FloorProfiles[Join].EndFloorZ-FloorProfiles[Join].StartFloorZ)/NextXY:0.0;
                        const double RawSlopeDelta=std::abs(RawPrevSlope-RawNextSlope);
                        if(FloorProfiles[Join-1].bLedgeTransition||FloorProfiles[Join].bLedgeTransition)
                        {
                            ++Stats.MarkedFloorSlopeJoins;
                            Stats.MaxMarkedFloorSlopeDelta=std::max(Stats.MaxMarkedFloorSlopeDelta,RawSlopeDelta);
                        }
                        else
                        {
                            Stats.MaxRawUnmarkedFloorSlopeDelta=std::max(Stats.MaxRawUnmarkedFloorSlopeDelta,RawSlopeDelta);
                        }
                    }
                    ++Stats.CapsuleJoinCount;
                }
                if(Points.Num()>1)
                {
                    const size_t SegmentPrimitiveStart=Primitives.size();
                    for(int32 I=0;I+1<Points.Num();++I)
                    {
                        const float R0=SmoothedControlRadius(Radii,I,T.RadiusA,2.0f);
                        const float R1=SmoothedControlRadius(Radii,I+1,T.RadiusB,2.0f);
                        const double RadiusSlope0=ControlRadiusSlopeAt(Points,Radii,I,T.RadiusA,2.0f);
                        const double RadiusSlope1=ControlRadiusSlopeAt(Points,Radii,I+1,T.RadiusB,2.0f);
                        Stats.SmoothedProfileRadiusStep=std::max(Stats.SmoothedProfileRadiusStep,
                            static_cast<double>(std::abs(R1-R0)));
                        const FVector Segment=Points[I+1]-Points[I];
                        const double SegmentLength=std::sqrt(Segment.X*Segment.X+Segment.Y*Segment.Y+Segment.Z*Segment.Z);
                        const double Extension=ContinuousJoinExtension(R0,R1,SegmentLength);
                        const float ZA=RawFloor[static_cast<size_t>(I)];
                        const float ZB=RawFloor[static_cast<size_t>(I+1)];
                        AddTaperedCapsule(Points[I],Points[I+1],R0,R1,true,false,true,ZA,ZB,
                            I>0?Extension:0.0,I+2<Points.Num()?Extension:0.0);
                        if(!MarkedFloorSegments[static_cast<size_t>(I)])
                        {
                            SetCurvedFloorClip(Primitives.back(),Points[I],Points[I+1],ZA,ZB,
                                FloorTangents[static_cast<size_t>(I)],FloorTangents[static_cast<size_t>(I+1)]);
                            ++Stats.CurvedFloorSegments;
                            Stats.MaxFloorCurveDeviation=std::max(Stats.MaxFloorCurveDeviation,
                                Primitives.back().FloorCurveMaxDeviation);
                        }
                        FitCapsuleRadiiAtEndpoints(Primitives.back(),Points[I],R0,Points[I+1],R1);
                        SetCapsuleRadiusProfile(Primitives.back(),Points[I],R0,RadiusSlope0,
                            Points[I+1],R1,RadiusSlope1);
                    }
                    for(int32 Join=1;Join+1<Points.Num();++Join)
                    {
                        const FPrimitive& Previous=Primitives[SegmentPrimitiveStart+static_cast<size_t>(Join-1)];
                        const FPrimitive& Next=Primitives[SegmentPrimitiveStart+static_cast<size_t>(Join)];
                        const float PreviousEnd=CapsuleRadiusAt(Previous,Points[Join]);
                        const float NextStart=CapsuleRadiusAt(Next,Points[Join]);
                        Stats.MaxCapsuleJoinRadiusDelta=std::max(Stats.MaxCapsuleJoinRadiusDelta,
                            static_cast<double>(std::abs(PreviousEnd-NextStart)));
                        Stats.MaxCapsuleJoinRadiusSlopeDelta=std::max(
                            Stats.MaxCapsuleJoinRadiusSlopeDelta,
                            std::abs(Previous.CapsuleRadiusSlopeEnd-Next.CapsuleRadiusSlopeStart));
                        if(!MarkedFloorSegments[static_cast<size_t>(Join-1)]
                            &&!MarkedFloorSegments[static_cast<size_t>(Join)])
                        {
                            const double TangentDelta=Previous.bCurvedFloorClip&&Next.bCurvedFloorClip
                                ?std::abs(Previous.FloorCurveTangentEnd-Next.FloorCurveTangentStart)
                                :std::abs(FloorSlopeAtEndpoint(Previous,true)-FloorSlopeAtEndpoint(Next,false));
                            const double DirectionalSlopeChange=std::abs(FloorSlopeAtEndpoint(Previous,true)
                                -FloorSlopeAtEndpoint(Next,false));
                            Stats.MaxCapsuleJoinFloorSlopeDelta=std::max(
                                Stats.MaxCapsuleJoinFloorSlopeDelta,TangentDelta);
                            Stats.MaxUnmarkedFloorSlopeDelta=std::max(
                                Stats.MaxUnmarkedFloorSlopeDelta,TangentDelta);
                            Stats.MaxUnmarkedFloorDirectionSlopeDelta=std::max(
                                Stats.MaxUnmarkedFloorDirectionSlopeDelta,DirectionalSlopeChange);
                        }
                    }
                }
                else AddCapsule(T.EndpointA,T.EndpointB,std::max(T.RadiusA,T.RadiusB),true);
            }
            ActiveBlendRadius=RoomBlendRadius;
            ActiveFloorBlendRadius=0.0;
            for(const FCachedPit& P:Cache.Pits)
                AddCapsule(FVector(P.CenterX,P.CenterY,P.TopZ-P.Depth),
                    FVector(P.CenterX,P.CenterY,P.TopZ),std::max(2.0f,P.Radius),true);
            for(const FCachedChimney& C:Cache.Chimneys)
                AddCapsule(FVector(C.CenterX,C.CenterY,C.BottomZ),
                    FVector(C.CenterX,C.CenterY,C.BottomZ+C.Height),std::max(2.0f,C.Radius),true);
            for(const FCachedColumn& C:Cache.Columns)
            {
                FVector Center(C.CenterX,C.CenterY,(LayerRanges[0].Lo+LayerRanges[0].Hi)*0.5);
                const double HalfZ=std::max(8.0,(LayerRanges[0].Hi-LayerRanges[0].Lo)*0.48);
                AddEllipsoid(Center,{std::max(2.0f,C.Radius),std::max(2.0f,C.Radius),HalfZ},
                    {1,0,0},{0,1,0},{0,0,1},false);
            }
        }
        const double PassageJoinBlend=ActiveBlendRadius;ActiveBlendRadius=JoinBlendRadius;
        ActiveFloorBlendRadius=FloorClipBlendRadius;
        for(const FVoxelPassage& P:Owner->Manager->GetPassages())
        {
            ++Stats.PassageCount;
            const TArray<FVector>& Points=P.ControlPoints;
            const TArray<float>& Radii=P.ControlRadii;
            const size_t PassagePrimitiveStart=Primitives.size();
            std::vector<float> RawFloor(static_cast<size_t>(Points.Num()));
            std::vector<uint8_t> MarkedFloorSegments(Points.Num()>0?static_cast<size_t>(Points.Num()-1):0,0);
            for(int32 I=0;I<Points.Num();++I)
                RawFloor[static_cast<size_t>(I)]=static_cast<float>(Points[I].Z
                    -SmoothedControlRadius(Radii,I,P.Radius,6.0f));
            if(P.NativeFloorProfileZ.Num()==Points.Num())
                for(int32 I=0;I<Points.Num();++I)RawFloor[static_cast<size_t>(I)]=P.NativeFloorProfileZ[I];
            for(int32 I=0;I<Points.Num();++I)
                RawFloor[static_cast<size_t>(I)]+=static_cast<float>(FloorSurfaceRelief(
                    Points[I].X,Points[I].Y,static_cast<uint32>(Owner->Seed)));
            const std::vector<double> FloorTangents=BuildFloorTangents(
                Points,RawFloor,MarkedFloorSegments);
            for(int32 Join=1;Join+1<Points.Num();++Join)
            {
                const float RPrev=std::max(6.0f,Radii.IsValidIndex(Join-1)?Radii[Join-1]:P.Radius);
                const float RHere=std::max(6.0f,Radii.IsValidIndex(Join)?Radii[Join]:P.Radius);
                const float RNext=std::max(6.0f,Radii.IsValidIndex(Join+1)?Radii[Join+1]:P.Radius);
                Stats.LegacySegmentRadiusStep=std::max(Stats.LegacySegmentRadiusStep,
                    static_cast<double>(std::abs(std::max(RPrev,RHere)-std::max(RHere,RNext))));
                if(P.NativeFloorProfileZ.IsValidIndex(Join-1)&&P.NativeFloorProfileZ.IsValidIndex(Join))
                {
                    const FVector PrevD=Points[Join]-Points[Join-1];
                    const FVector NextD=Points[Join+1]-Points[Join];
                    const double PrevXY=std::sqrt(PrevD.X*PrevD.X+PrevD.Y*PrevD.Y);
                    const double NextXY=std::sqrt(NextD.X*NextD.X+NextD.Y*NextD.Y);
                    const double PrevSlope=PrevXY>1.0e-6
                        ?(P.NativeFloorProfileZ[Join]-P.NativeFloorProfileZ[Join-1])/PrevXY:0.0;
                    const double NextSlope=NextXY>1.0e-6
                        ?(P.NativeFloorProfileZ[Join+1]-P.NativeFloorProfileZ[Join])/NextXY:0.0;
                    const double RawSlopeDelta=std::abs(PrevSlope-NextSlope);
                    Stats.MaxRawUnmarkedFloorSlopeDelta=std::max(Stats.MaxRawUnmarkedFloorSlopeDelta,RawSlopeDelta);
                }
                ++Stats.CapsuleJoinCount;
            }
            if(Points.Num()>1)
            {
                for(int32 I=0;I+1<Points.Num();++I)
                {
                    const float RadiusA=SmoothedControlRadius(Radii,I,P.Radius,6.0f);
                    const float RadiusB=SmoothedControlRadius(Radii,I+1,P.Radius,6.0f);
                    const double RadiusSlopeA=ControlRadiusSlopeAt(Points,Radii,I,P.Radius,6.0f);
                    const double RadiusSlopeB=ControlRadiusSlopeAt(Points,Radii,I+1,P.Radius,6.0f);
                    Stats.SmoothedProfileRadiusStep=std::max(Stats.SmoothedProfileRadiusStep,
                        static_cast<double>(std::abs(RadiusB-RadiusA)));
                    const FVector Segment=Points[I+1]-Points[I];
                    const double SegmentLength=std::sqrt(Segment.X*Segment.X+Segment.Y*Segment.Y+Segment.Z*Segment.Z);
                    const float ZA=RawFloor[static_cast<size_t>(I)];
                    const float ZB=RawFloor[static_cast<size_t>(I+1)];
                    const double Extension=ContinuousJoinExtension(RadiusA,RadiusB,SegmentLength);
                    const double CoreExtension=ContinuousJoinExtension(std::min(2.5f,RadiusA*0.45f),
                        std::min(2.5f,RadiusB*0.45f),SegmentLength);
                    const double ExtendStart=I>0?Extension:0.0,ExtendEnd=I+2<Points.Num()?Extension:0.0;
                    const double CoreExtendStart=I>0?CoreExtension:0.0,CoreExtendEnd=I+2<Points.Num()?CoreExtension:0.0;
                    AddTaperedCapsule(Points[I],Points[I+1],RadiusA,RadiusB,true,false,true,ZA,ZB,
                        ExtendStart,ExtendEnd);
                    SetCurvedFloorClip(Primitives.back(),Points[I],Points[I+1],ZA,ZB,
                        FloorTangents[static_cast<size_t>(I)],FloorTangents[static_cast<size_t>(I+1)]);
                    ++Stats.CurvedFloorSegments;
                    Stats.MaxFloorCurveDeviation=std::max(Stats.MaxFloorCurveDeviation,
                        Primitives.back().FloorCurveMaxDeviation);
                    FitCapsuleRadiiAtEndpoints(Primitives.back(),Points[I],RadiusA,Points[I+1],RadiusB);
                    SetCapsuleRadiusProfile(Primitives.back(),Points[I],RadiusA,RadiusSlopeA,
                        Points[I+1],RadiusB,RadiusSlopeB);
                    const float CoreRadiusA=std::min(2.5f,RadiusA*0.45f);
                    const float CoreRadiusB=std::min(2.5f,RadiusB*0.45f);
                    const double CoreSlopeA=CoreRadiusA<2.5f-1.0e-4?RadiusSlopeA*0.45:0.0;
                    const double CoreSlopeB=CoreRadiusB<2.5f-1.0e-4?RadiusSlopeB*0.45:0.0;
                    AddTaperedCapsule(Points[I],Points[I+1],CoreRadiusA,
                        CoreRadiusB,true,true,true,ZA,ZB,
                        CoreExtendStart,CoreExtendEnd);
                    SetCurvedFloorClip(Primitives.back(),Points[I],Points[I+1],ZA,ZB,
                        FloorTangents[static_cast<size_t>(I)],FloorTangents[static_cast<size_t>(I+1)]);
                    ++Stats.CurvedFloorSegments;
                    Stats.MaxFloorCurveDeviation=std::max(Stats.MaxFloorCurveDeviation,
                        Primitives.back().FloorCurveMaxDeviation);
                    FitCapsuleRadiiAtEndpoints(Primitives.back(),Points[I],CoreRadiusA,
                        Points[I+1],CoreRadiusB);
                    SetCapsuleRadiusProfile(Primitives.back(),Points[I],CoreRadiusA,CoreSlopeA,
                        Points[I+1],CoreRadiusB,CoreSlopeB);
                }
                for(int32 Join=1;Join+1<Points.Num();++Join)
                {
                    const size_t PreviousBase=PassagePrimitiveStart+static_cast<size_t>(Join-1)*2;
                    const size_t NextBase=PassagePrimitiveStart+static_cast<size_t>(Join)*2;
                    for(size_t Tube=0;Tube<2;++Tube)
                    {
                        const FPrimitive& Previous=Primitives[PreviousBase+Tube];
                        const FPrimitive& Next=Primitives[NextBase+Tube];
                        const float PreviousEnd=CapsuleRadiusAt(Previous,Points[Join]);
                        const float NextStart=CapsuleRadiusAt(Next,Points[Join]);
                        Stats.MaxCapsuleJoinRadiusDelta=std::max(Stats.MaxCapsuleJoinRadiusDelta,
                            static_cast<double>(std::abs(PreviousEnd-NextStart)));
                        Stats.MaxCapsuleJoinRadiusSlopeDelta=std::max(
                            Stats.MaxCapsuleJoinRadiusSlopeDelta,
                            std::abs(Previous.CapsuleRadiusSlopeEnd-Next.CapsuleRadiusSlopeStart));
                        if(Tube==0)
                        {
                            const double TangentDelta=Previous.bCurvedFloorClip&&Next.bCurvedFloorClip
                                ?std::abs(Previous.FloorCurveTangentEnd-Next.FloorCurveTangentStart)
                                :std::abs(FloorSlopeAtEndpoint(Previous,true)-FloorSlopeAtEndpoint(Next,false));
                            const double DirectionalSlopeChange=std::abs(FloorSlopeAtEndpoint(Previous,true)
                                -FloorSlopeAtEndpoint(Next,false));
                            Stats.MaxCapsuleJoinFloorSlopeDelta=std::max(
                                Stats.MaxCapsuleJoinFloorSlopeDelta,TangentDelta);
                            Stats.MaxUnmarkedFloorSlopeDelta=std::max(
                                Stats.MaxUnmarkedFloorSlopeDelta,TangentDelta);
                            Stats.MaxUnmarkedFloorDirectionSlopeDelta=std::max(
                                Stats.MaxUnmarkedFloorDirectionSlopeDelta,DirectionalSlopeChange);
                        }
                    }
                }
            }
            else AddCapsule(P.UpperPoint,P.LowerPoint,std::max(6.0f,P.Radius),true,true);
            AddLanding(P.UpperLanding,true);
            AddLanding(P.LowerLanding,true);
        }
        ActiveBlendRadius=PassageJoinBlend;ActiveFloorBlendRadius=0.0;
        if(Owner->Settings->bOpenSurfaceEntry)
        {
            const double Radius=std::max(8.0f,Owner->Settings->OriginSpineRadius);
            const double Top=LayerRanges[0].Hi;
            const FVector A(0,0,Top),B(0,0,std::max<double>(LayerRanges[0].Lo,Top-256.0));
            AddCapsule(A,B,static_cast<float>(Radius),true);
        }
    }

    void AddLanding(const FVoxelPassageLanding& L,bool bCritical=false)
    {
        const double HalfZ=std::max(2.0,(L.CeilingZ-L.FloorZ)*0.5);
        const FVector Center(L.StandingPoint.X,L.StandingPoint.Y,(L.CeilingZ+L.FloorZ)*0.5);
        const double Width=std::max(6.0f,L.HalfWidth);
        AddRoundedBox(Center,{Width,Width,HalfZ},2.0,true);
        if(bCritical)
        {
            const size_t Before=Primitives.size();
            AddEllipsoid(Center,{2.5,2.5,HalfZ},{1,0,0},{0,1,0},{0,0,1},true);
            for(size_t I=Before;I<Primitives.size();++I)Primitives[I].bCritical=true;
        }
        const FVector Floor(Center.X,Center.Y,L.FloorZ-1.25);
        AddEllipsoid(Floor,{Width,Width,1.25},{1,0,0},{0,1,0},{0,0,1},false);
    }

    void BuildRoundTunnelPrimitives()
    {
        const double PreviousBlendRadius=ActiveBlendRadius;ActiveBlendRadius=0.0;
        constexpr double Pitch=144.0,RadiusBase=28.0,HalfHeight=26.0;
        const int32 CX0=FloorDiv(Min[0],144)-1,CY0=FloorDiv(Min[1],144)-1;
        const int32 CX1=FloorDiv(Max[0],144)+1,CY1=FloorDiv(Max[1],144)+1;
        const uint32 Seed=static_cast<uint32>(Owner->Seed);
        auto Unit=[](uint32 H){return static_cast<double>(FStandaloneResearchTunnels::Unit(H));};
        for(int32 CY=CY0;CY<=CY1;++CY)for(int32 CX=CX0;CX<=CX1;++CX)
        for(int32 Axis=0;Axis<2;++Axis)
        {
            const uint32 H=FStandaloneResearchTunnels::EdgeHash(CX,CY,Axis,Seed);
            if(Unit(H^0x68E31DA4u)>0.56)continue;
            const int32 NX=CX+(Axis==0?1:0),NY=CY+(Axis==1?1:0);
            auto Cell=[&](int32 X,int32 Y,double& OX,double& OY)
            {
                const uint32 CH=FStandaloneResearchTunnels::EdgeHash(X,Y,3,Seed);
                OX=(X+0.5)*Pitch+(Unit(CH^0xA511E9B3u)-0.5)*Pitch*0.18;
                OY=(Y+0.5)*Pitch+(Unit(CH^0x63D83595u)-0.5)*Pitch*0.18;
            };
            double AX,AY,BX,BY;Cell(CX,CY,AX,AY);Cell(NX,NY,BX,BY);
            const double DX=BX-AX,DY=BY-AY,L2=DX*DX+DY*DY,L=std::sqrt(std::max(1.0e-8,L2));
            const double NXrm=-DY/L,NYrm=DX/L;
            const double BA=(Unit(H^0xB5297A4Du)*2.0-1.0)*Pitch*0.16;
            const double BB=(Unit(H^0x1B56C4E9u)*2.0-1.0)*Pitch*0.16;
            const std::array<double,3> P0={AX,AY,FStandaloneResearchTunnels::WorldCenterZVoxels};
            const std::array<double,3> C1={AX+DX/3.0+NXrm*BA,AY+DY/3.0+NYrm*BA,-80.0};
            const std::array<double,3> C2={AX+DX*2.0/3.0+NXrm*BB,AY+DY*2.0/3.0+NYrm*BB,-80.0};
            const std::array<double,3> P3={BX,BY,FStandaloneResearchTunnels::WorldCenterZVoxels};
            const uint32 PhaseHash=H^0xA511E9B3u;
            const double Phase=Unit(PhaseHash)*6.2831853071795864769;
            auto Bezier=[&](double T)
            {
                const double O=1.0-T;
                std::array<double,3> P{};
                for(int32 A=0;A<3;++A)
                    P[A]=O*O*O*P0[A]+3*O*O*T*C1[A]+3*O*T*T*C2[A]+T*T*T*P3[A];
                const double R=RadiusBase*(0.82+0.36*Unit(H^0x63D83595u))
                    *(1.0+0.22*VoxelMath::DetSin(static_cast<float>(6.2831853071795864769*T+Phase)));
                const double VW=4.0*VoxelMath::DetSin(static_cast<float>(3.141592653589793*T))
                    *VoxelMath::DetSin(static_cast<float>(6.2831853071795864769*T+Phase));
                P[2]+=VW;
                return std::pair{P,R};
            };
            constexpr int32 Segments=12;
            auto Last=Bezier(0.0);
            for(int32 I=1;I<=Segments;++I)
            {
                const auto Next=Bezier(static_cast<double>(I)/Segments);
                const FVector A(Last.first[0],Last.first[1],Last.first[2]);
                const FVector B(Next.first[0],Next.first[1],Next.first[2]);
                const FVector Segment=B-A;
                const double SegmentLength=std::sqrt(Segment.X*Segment.X+Segment.Y*Segment.Y+Segment.Z*Segment.Z);
                const double Extension=ContinuousJoinExtension(static_cast<float>(Last.second),
                    static_cast<float>(Next.second),SegmentLength);
                AddTaperedCapsule(A,B,static_cast<float>(Last.second),static_cast<float>(Next.second),
                    true,false,false,0.0,0.0,I>1?Extension:0.0,I<Segments?Extension:0.0);
                Last=Next;
            }
        }
        ActiveBlendRadius=PreviousBlendRadius;
    }

    void BuildSpatialBins(int32 BinResolution,bool bBuildAxisBins=true)
    {
        SpatialBins.clear();
        CriticalAirBins.clear();
        ActiveSpatialBinSize=BinResolution;
        // Primitive bounds overlap multiple cells. Reserve from the measured bins-per-primitive
        // density for each grid instead of repeatedly rehashing from a 2x primitive guess.
          const size_t BinsPerPrimitive=bBuildAxisBins?96u:24u;
        SpatialBins.reserve(std::max<size_t>(Primitives.size()*BinsPerPrimitive,1u));
        if(!bBuildAxisBins)
        {
            CriticalAirBins.reserve(Primitives.size()/16+1);
        }
        if(bBuildAxisBins)for(int32 Axis=0;Axis<3;++Axis)
        {
            FAxisBins& B=AxisBins[Axis];
            const int32 A=(Axis+1)%3,C=(Axis+2)%3;
            B.MinA=FloorDiv(Min[A],LineBinSize);B.MinB=FloorDiv(Min[C],LineBinSize);
            const int32 MaxA=FloorDiv(Max[A],LineBinSize),MaxB=FloorDiv(Max[C],LineBinSize);
            B.CountA=MaxA-B.MinA+1;B.CountB=MaxB-B.MinB+1;
            B.Items.resize(static_cast<size_t>(B.CountA)*B.CountB);
            for(uint32 Id=0;Id<Primitives.size();++Id)
            {
                const FPrimitive& P=Primitives[Id];
                const int32 A0=std::max(B.MinA,FloorDiv(P.Min[A],LineBinSize));
                const int32 A1=std::min(MaxA,FloorDiv(P.Max[A],LineBinSize));
                const int32 B0=std::max(B.MinB,FloorDiv(P.Min[C],LineBinSize));
                const int32 B1=std::min(MaxB,FloorDiv(P.Max[C],LineBinSize));
                if(A0>A1||B0>B1)continue;
                for(int32 IB=B0;IB<=B1;++IB)for(int32 IA=A0;IA<=A1;++IA)
                    B.Items[static_cast<size_t>(IB-B.MinB)*B.CountA+(IA-B.MinA)].push_back(Id);
            }
        }
        for(uint32 Id=0;Id<Primitives.size();++Id)
        {
            const FPrimitive& P=Primitives[Id];
            const bool bCriticalAir=P.bCarve&&P.bCritical;
            const int32 X0=FloorDiv(P.Min[0]-DetailBand,BinResolution),X1=FloorDiv(P.Max[0]+DetailBand,BinResolution);
            const int32 Y0=FloorDiv(P.Min[1]-DetailBand,BinResolution),Y1=FloorDiv(P.Max[1]+DetailBand,BinResolution);
            const int32 Z0=FloorDiv(P.Min[2]-DetailBand,BinResolution),Z1=FloorDiv(P.Max[2]+DetailBand,BinResolution);
            if(static_cast<int64>(X1-X0+1)*(Y1-Y0+1)*(Z1-Z0+1)>100000)continue;
            for(int32 Z=Z0;Z<=Z1;++Z)for(int32 Y=Y0;Y<=Y1;++Y)for(int32 X=X0;X<=X1;++X)
            {
                const FBinKey Key{X,Y,Z};
                SpatialBins[Key].push_back(Id);
                if(bCriticalAir)CriticalAirBins[Key].push_back(Id);
            }
        }
    }

    void BuildFarParentShapes()
    {
        constexpr int32 ParentPitch=2048;
        struct FParent { int32 GX=0,GY=0,Strate=0; FVector Center; double Rx=0,Ry=0,Rz=0; };
        const int32 GX0=FloorDiv(ClipMin[0],ParentPitch)-1,GY0=FloorDiv(ClipMin[1],ParentPitch)-1;
        const int32 GX1=FloorDiv(ClipMax[0],ParentPitch)+1,GY1=FloorDiv(ClipMax[1],ParentPitch)+1;
        std::vector<FParent> Parents;
        std::unordered_map<uint64,size_t> ByCell;
        const auto& Layout=Owner->Manager->GetLayout();
        const uint32 Seed=static_cast<uint32>(Owner->Seed);
        const auto Key=[](int32 GX,int32 GY,int32 Strate)
        {return (static_cast<uint64>(static_cast<uint32>(GX+524288)&0xFFFFFu)<<44)
            |(static_cast<uint64>(static_cast<uint32>(GY+524288)&0xFFFFFu)<<24)
            |(static_cast<uint64>(static_cast<uint32>(Strate)&0xFFFFu)<<8);};
        auto Hash=[](uint32 X)
        {X^=X>>16;X*=0x7FEB352Du;X^=X>>15;X*=0x846CA68Bu;X^=X>>16;return X;};
        auto Unit=[&](uint32 H){return static_cast<double>(Hash(H)>>8)*(1.0/16777216.0);};
        for(const FStrateSlot& Slot:Layout)
        {
            const double ZBottom=static_cast<double>(Slot.BottomChunkZ)*CHUNK_SIZE;
            const double ZTop=static_cast<double>(Slot.TopChunkZ+1)*CHUNK_SIZE;
            for(int32 GY=GY0;GY<=GY1;++GY)for(int32 GX=GX0;GX<=GX1;++GX)
            {
                uint32 H=Hash(Seed^static_cast<uint32>(GX)*0x9E3779B9u
                    ^static_cast<uint32>(GY)*0x85EBCA6Bu
                    ^static_cast<uint32>(Slot.StrateIndex)*0xC2B2AE35u);
                const double X=GX*static_cast<double>(ParentPitch)+ParentPitch*0.5
                    +(Unit(H^0xA511E9B3u)-0.5)*192.0;
                const double Y=GY*static_cast<double>(ParentPitch)+ParentPitch*0.5
                    +(Unit(H^0x63D83595u)-0.5)*192.0;
                const double Rx=300.0+Unit(H^0xB5297A4Du)*100.0;
                const double Ry=280.0+Unit(H^0x1B56C4E9u)*120.0;
                const double Rz=100.0+Unit(H^0x68E31DA4u)*60.0;
                if(X+Rx>=Min[0]&&X-Rx<=Max[0]&&Y+Ry>=Min[1]&&Y-Ry<=Max[1])continue;
                const double Z=ZBottom+(ZTop-ZBottom)*(0.35+0.30*Unit(H^0x27D4EB2Fu));
                FParent P{GX,GY,Slot.StrateIndex,FVector(X,Y,Z),Rx,Ry,Rz};
                const size_t Index=Parents.size();Parents.push_back(P);
                ByCell.emplace(Key(GX,GY,Slot.StrateIndex),Index);
                const std::array<double,3> Axis0={1,0,0},Axis1={0,1,0},Axis2={0,0,1};
                AddEllipsoid(P.Center,{P.Rx,P.Ry,P.Rz},Axis0,Axis1,Axis2,true);
            }
        }
        for(const FParent& P:Parents)
        {
            for(const auto [DX,DY]:std::array<std::pair<int32,int32>,2>{{{1,0},{0,1}}})
            {
                const auto It=ByCell.find(Key(P.GX+DX,P.GY+DY,P.Strate));
                if(It==ByCell.end())continue;
                const FVector A(P.Center.X,P.Center.Y,P.Center.Z);
                const FVector B=Parents[It->second].Center;
                const int32 MinX=static_cast<int32>(std::floor(std::min(A.X,B.X)-64.0));
                const int32 MaxX=static_cast<int32>(std::ceil(std::max(A.X,B.X)+64.0));
                const int32 MinY=static_cast<int32>(std::floor(std::min(A.Y,B.Y)-64.0));
                const int32 MaxY=static_cast<int32>(std::ceil(std::max(A.Y,B.Y)+64.0));
                if(MaxX>=Min[0]&&MinX<=Max[0]&&MaxY>=Min[1]&&MinY<=Max[1])continue;
                AddTaperedCapsule(A,B,64.0f,64.0f,true);
                Primitives.back().BlendRadius=0.0;
            }
        }
        if(!Parents.empty())
            std::printf("[FarParents] pitch=%d rooms=%zu links=%llu clip=(%d,%d)..(%d,%d)\n",
                ParentPitch,Parents.size(),static_cast<unsigned long long>(Primitives.size()),
                ClipMin[0],ClipMin[1],ClipMax[0],ClipMax[1]);
    }

    const std::vector<uint32>& GetLineCandidates(int32 Axis,int32 A,int32 B) const
    {
        static const std::vector<uint32> Empty;
        const FAxisBins& G=AxisBins[Axis];
        const int32 IA=FloorDiv(A,LineBinSize)-G.MinA,IB=FloorDiv(B,LineBinSize)-G.MinB;
        if(IA<0||IB<0||IA>=G.CountA||IB>=G.CountB)return Empty;
        return G.Items[static_cast<size_t>(IB)*G.CountA+IA];
    }

    FBinKey PointBin(double X,double Y,double Z) const
    {return {FloorDiv(static_cast<int32>(std::floor(X)),ActiveSpatialBinSize),FloorDiv(static_cast<int32>(std::floor(Y)),ActiveSpatialBinSize),FloorDiv(static_cast<int32>(std::floor(Z)),ActiveSpatialBinSize)};}

    bool ShapeContains(const FPrimitive& P,int32 X,int32 Y,int32 Z) const
    {
        bool bInside=false;
        if(P.bCapsule)
        {
            const double Position[3]={static_cast<double>(X),static_cast<double>(Y),static_cast<double>(Z)};
            double V[3],D[3];double D2=0.0,VD=0.0,V2=0.0;
            for(int32 A=0;A<3;++A)
            {
                const double Start=P.CapsuleStartSub[A]/static_cast<double>(SubPerVoxel);
                const double End=P.CapsuleEndSub[A]/static_cast<double>(SubPerVoxel);
                V[A]=Position[A]-Start;D[A]=End-Start;
                D2+=D[A]*D[A];VD+=V[A]*D[A];V2+=V[A]*V[A];
            }
            const double R0=P.CapsuleStartRadius,R1=P.CapsuleEndRadius,DR=R1-R0;
            const double QA=D2-DR*DR,QB=-2.0*(VD+R0*DR),QC=V2-R0*R0;
            double Minimum=std::min(QC,QA+QB+QC);
            if(QA>1.0e-12)
            {
                const double T=std::clamp(-QB/(2.0*QA),0.0,1.0);
                Minimum=std::min(Minimum,(QA*T+QB)*T+QC);
            }
            bInside=Minimum<=0.0;
        }
        else if(P.bRoundedBox)
        {
            const double D[3]={std::abs(X-P.CenterSub[0]/16.0)-P.RoundedHalfExtent[0],
                std::abs(Y-P.CenterSub[1]/16.0)-P.RoundedHalfExtent[1],
                std::abs(Z-P.CenterSub[2]/16.0)-P.RoundedHalfExtent[2]};
            const double Ox=std::max(D[0],0.0),Oy=std::max(D[1],0.0),Oz=std::max(D[2],0.0);
            const double Outside=std::sqrt(Ox*Ox+Oy*Oy+Oz*Oz);
            const double Inside=std::min(std::max({D[0],D[1],D[2]}),0.0);
            bInside=Outside+Inside-P.Rounding<=0.0;
        }
        else
        {
            const int32 World[3]={X,Y,Z};int64 D[3];
            for(int32 A=0;A<3;++A)D[A]=static_cast<int64>(World[A])*SubPerVoxel-P.CenterSub[A];
            int64 Q=0;
            for(int32 I=0;I<3;++I)for(int32 J=0;J<3;++J)Q+=P.Q[I][J]*D[I]*D[J];
            bInside=Q<=static_cast<int64>(QScale);
        }
        if(!bInside)return false;
        if(P.bHasFloorClip)
        {
            const double Floor=FloorAt(P,X,Y);
            if(P.FloorBlendRadius<=0.0)return Z>=Floor;
            return SmoothMaxDistance(ShapeDistance(P,X,Y,Z),static_cast<float>(Floor-Z),
                static_cast<float>(P.FloorBlendRadius))<=0.0f;
        }
        return true;
    }

    static float CapsuleRadiusAtT(const FPrimitive& P,double T,double AxisLength)
    {
        T=std::clamp(T,0.0,1.0);
        if(!P.bSmoothCapsuleRadius)
            return static_cast<float>(P.CapsuleStartRadius
                +(P.CapsuleEndRadius-P.CapsuleStartRadius)*T);
        const double T0=P.CapsuleProfileStartT,T1=P.CapsuleProfileEndT;
        if(T<T0)return static_cast<float>(P.CapsuleProfileStartRadius
            +P.CapsuleRadiusSlopeStart*AxisLength*(T-T0));
        if(T>T1)return static_cast<float>(P.CapsuleProfileEndRadius
            +P.CapsuleRadiusSlopeEnd*AxisLength*(T-T1));
        const double SpanT=std::max(1.0e-9,T1-T0);
        const double U=std::clamp((T-T0)/SpanT,0.0,1.0),U2=U*U,U3=U2*U;
        const double SpanDistance=AxisLength*SpanT;
        const double H00=2.0*U3-3.0*U2+1.0;
        const double H10=U3-2.0*U2+U;
        const double H01=-2.0*U3+3.0*U2;
        const double H11=U3-U2;
        return static_cast<float>(H00*P.CapsuleProfileStartRadius
            +H10*SpanDistance*P.CapsuleRadiusSlopeStart
            +H01*P.CapsuleProfileEndRadius
            +H11*SpanDistance*P.CapsuleRadiusSlopeEnd);
    }

    float ShapeDistance(const FPrimitive& P,double X,double Y,double Z) const
    {
        if(P.bCapsule)
        {
            const double Position[3]={X,Y,Z};
            double V[3],D[3],D2=0.0,VD=0.0;
            for(int32 A=0;A<3;++A)
            {
                V[A]=Position[A]-P.CapsuleStartWorld[A];
                D[A]=P.CapsuleAxisDelta[A];
                VD+=V[A]*D[A];
            }
            D2=P.CapsuleAxisLengthSquared;
            const double T=D2>1.0e-12?std::clamp(VD/D2,0.0,1.0):0.0;
            double DistanceSquared=0.0;
            for(int32 A=0;A<3;++A)
            {const double E=V[A]-T*D[A];DistanceSquared+=E*E;}
            const double Radius=CapsuleRadiusAtT(P,T,P.CapsuleAxisLength);
            return static_cast<float>(std::sqrt(DistanceSquared)-Radius);
        }
        if(P.bRoundedBox)
        {
            const double D[3]={std::abs(X-P.CenterSub[0]/16.0)-P.RoundedHalfExtent[0],
                std::abs(Y-P.CenterSub[1]/16.0)-P.RoundedHalfExtent[1],
                std::abs(Z-P.CenterSub[2]/16.0)-P.RoundedHalfExtent[2]};
            const double Ox=std::max(D[0],0.0),Oy=std::max(D[1],0.0),Oz=std::max(D[2],0.0);
            const double Outside=std::sqrt(Ox*Ox+Oy*Oy+Oz*Oz);
            const double Inside=std::min(std::max({D[0],D[1],D[2]}),0.0);
            return static_cast<float>(Outside+Inside-P.Rounding);
        }
        const double D[3]={X-P.CenterSub[0]/16.0,Y-P.CenterSub[1]/16.0,Z-P.CenterSub[2]/16.0};
        double Q=0.0;
        for(int32 I=0;I<3;++I)for(int32 J=0;J<3;++J)Q+=static_cast<double>(P.Q[I][J])/QScale*D[I]*D[J]*256.0;
        return static_cast<float>((std::sqrt(std::max(0.0,Q))-1.0)*P.MinRadius);
    }

    bool SolveLine(const FPrimitive& P,int32 Axis,const int32 Fixed[3],FInterval& Out) const
    {
        const int32 A=(Axis+1)%3,B=(Axis+2)%3;
        const int32 Lo=P.Min[Axis],Hi=P.Max[Axis];
        if(Fixed[A]<P.Min[A]||Fixed[A]>P.Max[A]||Fixed[B]<P.Min[B]||Fixed[B]>P.Max[B])return false;
        int32 L=Lo,H=Hi;
        if(!P.bRoundedBox&&!P.bCapsule)
        {
            double CoefA=static_cast<double>(P.Q[Axis][Axis]);
            double CoefB=0.0,CoefC=-static_cast<double>(QScale);
            int64 D[3];D[Axis]=0;
            for(int32 J=0;J<3;++J)if(J!=Axis)D[J]=static_cast<int64>(Fixed[J])*SubPerVoxel-P.CenterSub[J];
            for(int32 J=0;J<3;++J)if(J!=Axis)
            {
                CoefB+=2.0*P.Q[Axis][J]*D[J];
                for(int32 K=0;K<3;++K)if(K!=Axis)CoefC+=static_cast<double>(P.Q[J][K])*D[J]*D[K];
            }
            if(!(CoefA>0.0))return false;
            const double Disc=CoefB*CoefB-4.0*CoefA*CoefC;
            if(Disc<0.0)return false;
            const double S=std::sqrt(Disc);
            const double U0=(-CoefB-S)/(2.0*CoefA),U1=(-CoefB+S)/(2.0*CoefA);
            const double T0=(U0+P.CenterSub[Axis])/SubPerVoxel;
            const double T1=(U1+P.CenterSub[Axis])/SubPerVoxel;
            L=std::max(Lo,static_cast<int32>(std::floor(T0))-1);
            H=std::min(Hi,static_cast<int32>(std::ceil(T1))+1);
        }
        if(L>H)return false;
        if(P.bCarve&&P.bHasFloorClip&&(!P.bCurvedFloorClip||Axis==2))
        {
            double X0=Fixed[0],Y0=Fixed[1],Z0=Fixed[2];
            if(Axis==0)X0=0.0;else if(Axis==1)Y0=0.0;else Z0=0.0;
            const double Floor0=FloorAt(P,X0,Y0);
            const double SmoothAllowance=P.FloorBlendRadius>0.0?P.FloorBlendRadius/6.0:0.0;
            const double Margin0=Z0-(Floor0-SmoothAllowance);
            const double MarginSlope=Axis==0?-P.FloorSlopeX:Axis==1?-P.FloorSlopeY:1.0;
            if(std::abs(MarginSlope)<1.0e-12)
            {if(Margin0<0.0)return false;}
            else
            {
                const double Boundary=-Margin0/MarginSlope;
                if(MarginSlope>0.0)L=std::max(L,static_cast<int32>(std::ceil(Boundary)));
                else H=std::min(H,static_cast<int32>(std::floor(Boundary)));
            }
        }
        if(L>H)return false;
        int32 Mid=L+(H-L)/2;
        if(P.bCapsule)
        {
            double V2=0.0,VD=0.0,D2=0.0;
            for(int32 J=0;J<3;++J)if(J!=Axis)
            {
                const double Start=P.CapsuleStartSub[J]/static_cast<double>(SubPerVoxel);
                const double End=P.CapsuleEndSub[J]/static_cast<double>(SubPerVoxel);
                const double V=Fixed[J]-Start,D=End-Start;
                V2+=V*V;VD+=V*D;D2+=D*D;
            }
            const double R0=P.CapsuleStartRadius,DR=P.CapsuleEndRadius-R0;
            const double QA=D2-DR*DR,QB=-2.0*(VD+R0*DR),QC=V2-R0*R0;
            const auto Value=[&](double T){return (QA*T+QB)*T+QC;};
            double TBest=Value(1.0)<Value(0.0)?1.0:0.0;
            if(QA>1.0e-12)
            {
                const double T=std::clamp(-QB/(2.0*QA),0.0,1.0);
                if(Value(T)<Value(TBest))TBest=T;
            }
            if(Value(TBest)>0.0)return false;
            const double Start=P.CapsuleStartSub[Axis]/static_cast<double>(SubPerVoxel);
            const double End=P.CapsuleEndSub[Axis]/static_cast<double>(SubPerVoxel);
            Mid=std::clamp(static_cast<int32>(std::llround(Start+(End-Start)*TBest)),L,H);
        }
        if(!ShapeContains(P,Axis==0?Mid:Fixed[0],Axis==1?Mid:Fixed[1],Axis==2?Mid:Fixed[2]))
        {
            bool Found=false;
            const int32 SearchLimit=P.bCapsule?std::min(H-L,8):2;
            for(int32 Delta=1;Delta<=SearchLimit&&!Found;++Delta)
            {
                if(Mid-Delta>=L&&ShapeContains(P,Axis==0?Mid-Delta:Fixed[0],Axis==1?Mid-Delta:Fixed[1],Axis==2?Mid-Delta:Fixed[2])){Mid-=Delta;Found=true;}
                else if(Mid+Delta<=H&&ShapeContains(P,Axis==0?Mid+Delta:Fixed[0],Axis==1?Mid+Delta:Fixed[1],Axis==2?Mid+Delta:Fixed[2])){Mid+=Delta;Found=true;}
            }
            if(!Found&&P.bCapsule)
            {
                for(int32 Sample=L;Sample<=H;Sample+=2)
                    if(ShapeContains(P,Axis==0?Sample:Fixed[0],Axis==1?Sample:Fixed[1],Axis==2?Sample:Fixed[2]))
                    {Mid=Sample;Found=true;break;}
                if(!Found&&((H-L)&1)
                    &&ShapeContains(P,Axis==0?H:Fixed[0],Axis==1?H:Fixed[1],Axis==2?H:Fixed[2]))
                {Mid=H;Found=true;}
            }
            if(!Found)return false;
        }
        int32 Left=L,Right=Mid;
        while(Left<Right)
        {
            const int32 M=Left+(Right-Left)/2;
            const bool Inside=ShapeContains(P,Axis==0?M:Fixed[0],Axis==1?M:Fixed[1],Axis==2?M:Fixed[2]);
            if(Inside)Right=M;else Left=M+1;
        }
        const int32 First=Left;Left=Mid;Right=H;
        while(Left<Right)
        {
            const int32 M=Left+(Right-Left+1)/2;
            const bool Inside=ShapeContains(P,Axis==0?M:Fixed[0],Axis==1?M:Fixed[1],Axis==2?M:Fixed[2]);
            if(Inside)Left=M;else Right=M-1;
        }
        Out={First,Left};return true;
    }

    static void MergeIntervals(std::vector<FInterval>& Ranges)
    {
        if(Ranges.size()<2)return;
        std::sort(Ranges.begin(),Ranges.end(),[](const FInterval& A,const FInterval& B){return A.Lo==B.Lo?A.Hi<B.Hi:A.Lo<B.Lo;});
        size_t Out=0;
        for(size_t I=1;I<Ranges.size();++I)
        {
            if(Ranges[I].Lo<=Ranges[Out].Hi+1)Ranges[Out].Hi=std::max(Ranges[Out].Hi,Ranges[I].Hi);
            else Ranges[++Out]=Ranges[I];
        }
        Ranges.resize(Out+1);
    }

    static void FillBits(std::vector<uint64>& Words,size_t Base,int32 Lo,int32 Hi,bool bSet)
    {
        if(Lo>Hi)return;
        const size_t First=static_cast<size_t>(Lo)>>6,Last=static_cast<size_t>(Hi)>>6;
        for(size_t W=First;W<=Last;++W)
        {
            const uint32 Bit0=W==First?(static_cast<uint32>(Lo)&63u):0u;
            const uint32 Bit1=W==Last?(static_cast<uint32>(Hi)&63u):63u;
            const uint64 Mask=(~uint64{0}<<Bit0)&(~uint64{0}>>(63u-Bit1));
            if(bSet)Words[Base+W]|=Mask;else Words[Base+W]&=~Mask;
        }
    }

    void BuildPackedRockIntervals(FParallelTiming& Timing)
    {
        const int32 NX=Size[0],NY=Size[1],NZ=Size[2];
        ParallelLines(static_cast<size_t>(NX)*NY,[&](size_t Line)
        {
            const int32 X=Min[0]+static_cast<int32>(Line%NX),Y=Min[1]+static_cast<int32>(Line/NX);
            std::vector<FInterval> Carve,Solid,ProtectedAir;
            Carve.reserve(24);Solid.reserve(8);ProtectedAir.reserve(4);
            const int32 Fixed[3]={X,Y,0};
            for(uint32 Id:GetLineCandidates(2,X,Y))
            {
                const FPrimitive& P=Primitives[Id];FInterval R;
                if(SolveLine(P,2,Fixed,R))
                {
                    R.Lo=std::max(R.Lo,Min[2])-Min[2];R.Hi=std::min(R.Hi,Max[2])-Min[2];
                    if(R.Lo<=R.Hi)
                    {
                        if(!P.bCarve)Solid.push_back(R);
                        else if(P.bCritical)ProtectedAir.push_back(R);
                        else Carve.push_back(R);
                    }
                }
            }
            MergeIntervals(Carve);MergeIntervals(Solid);MergeIntervals(ProtectedAir);
            const size_t Base=Line*WordsPerZ;
            const int32 ZLo=LayerRanges[0].Lo, ZHi=LayerRanges[0].Hi;
            const int32 BaseLo=std::max(ZLo,Min[2])-Min[2],BaseHi=std::min(ZHi,Max[2])-Min[2];
            if(BaseLo<=BaseHi)FillBits(RockWords,Base,BaseLo,BaseHi,true);
            for(const FInterval& R:Carve)FillBits(RockWords,Base,R.Lo,R.Hi,false);
            for(const FInterval& R:Solid)FillBits(RockWords,Base,R.Lo,R.Hi,true);
            // Critical passage air is the final authored override, including at landing floors.
            for(const FInterval& R:ProtectedAir)FillBits(RockWords,Base,R.Lo,R.Hi,false);
        },&Timing);
    }

    void BuildCandidateBands(FParallelTiming& AxisTiming,FParallelTiming& MergeTiming,
                             double& MergeSeconds)
    {
        const int32 N[3]={Size[0],Size[1],Size[2]};
        for(int32 Axis=0;Axis<3;++Axis)
        {
            const int32 A=(Axis+1)%3,B=(Axis+2)%3;
            const size_t Lines=static_cast<size_t>(Size[A])*Size[B];
            const size_t Words=(static_cast<size_t>(N[Axis])+63u)/64u;
            CandidateAxisWords[Axis].assign(Lines*Words,0);
            ParallelLines(Lines,[&,Axis,A,B,Words](size_t Line)
            {
                const int32 CA=Min[A]+static_cast<int32>(Line%Size[A]);
                const int32 CB=Min[B]+static_cast<int32>(Line/Size[A]);
                int32 Fixed[3]{};Fixed[A]=CA;Fixed[B]=CB;
                for(uint32 Id:GetLineCandidates(Axis,CA,CB))
                {
                    const FPrimitive& P=Primitives[Id];if(P.bCritical)continue;
                    FInterval R;if(!SolveLine(P,Axis,Fixed,R))continue;
                    const int32 Lo=std::max(Min[Axis],R.Lo-DetailBand);
                    const int32 Hi=std::min(Max[Axis],R.Lo+DetailBand);
                    const int32 Lo2=std::max(Min[Axis],R.Hi-DetailBand);
                    const int32 Hi2=std::min(Max[Axis],R.Hi+DetailBand);
                    const size_t Base=Line*Words;
                    FillBits(CandidateAxisWords[Axis],Base,Lo-Min[Axis],Hi-Min[Axis],true);
                    FillBits(CandidateAxisWords[Axis],Base,Lo2-Min[Axis],Hi2-Min[Axis],true);
                }
            },&AxisTiming);
        }
        const auto MergeStart=std::chrono::steady_clock::now();
        CandidateWords=CandidateAxisWords[2];
        for(int32 Axis=0;Axis<2;++Axis)
        {
            const int32 A=(Axis+1)%3,B=(Axis+2)%3;
            const size_t Words=(static_cast<size_t>(Size[Axis])+63u)/64u;
            ParallelLines(static_cast<size_t>(Size[A])*Size[B],[&,Axis,A,B,Words](size_t Line)
            {
                const int32 CA=Min[A]+static_cast<int32>(Line%Size[A]);
                const int32 CB=Min[B]+static_cast<int32>(Line/Size[A]);
                for(size_t W=0;W<Words;++W)
                {
                    uint64 Bits=CandidateAxisWords[Axis][Line*Words+W];
                    while(Bits)
                    {
                        const uint32 Bit=std::countr_zero(Bits);Bits&=Bits-1;
                        const int32 T=Min[Axis]+static_cast<int32>(W*64+Bit);
                        if(T>Max[Axis])continue;
                        int32 Coord[3]{};Coord[Axis]=T;Coord[A]=CA;Coord[B]=CB;
                        const size_t Column=static_cast<size_t>(Coord[1]-Min[1])*Size[0]
                            +static_cast<size_t>(Coord[0]-Min[0]);
                        const size_t Zi=static_cast<size_t>(Coord[2]-Min[2]);
                        std::atomic_ref<uint64>(CandidateWords[Column*WordsPerZ+(Zi>>6)])
                            .fetch_or(uint64{1}<<(Zi&63u),std::memory_order_relaxed);
                    }
                }
            },&MergeTiming);
        }
        const int32 LayerLo=LayerRanges.empty()?1:LayerRanges[0].Lo;
        const int32 LayerHi=LayerRanges.empty()?0:LayerRanges[0].Hi;
        if(LayerLo<=LayerHi)
        {
            for(size_t Line=0;Line<static_cast<size_t>(Size[0])*Size[1];++Line)
            {
                const size_t Base=Line*WordsPerZ;
                const int32 Lo0=std::max(Min[2],LayerLo-DetailBand)-Min[2];
                const int32 Hi0=std::min(Max[2],LayerLo+DetailBand)-Min[2];
                const int32 Lo1=std::max(Min[2],LayerHi-DetailBand)-Min[2];
                const int32 Hi1=std::min(Max[2],LayerHi+DetailBand)-Min[2];
                FillBits(CandidateWords,Base,Lo0,Hi0,true);
                FillBits(CandidateWords,Base,Lo1,Hi1,true);
            }
        }
        const double MergeElapsed=std::chrono::duration<double>(
            std::chrono::steady_clock::now()-MergeStart).count();
        MergeSeconds=std::max(0.0,MergeElapsed-MergeTiming.WallSeconds);
    }

    void ApplyDetailBand(FParallelTiming& Timing)
    {
        using Clock=std::chrono::steady_clock;
        const bool bProfile=ProfileDetailBand.load(std::memory_order_acquire);
        std::array<FDetailWorkerProfile,16> WorkerProfile{};
        std::array<double,16> WorkerActiveSeconds{};
        std::atomic<uint64> CandidateCount{0},ProxyCount{0},FlipCount{0};
        ParallelLines(static_cast<size_t>(Size[0])*Size[1],[&](size_t Line,int32 Worker)
        {
            const int32 X=Min[0]+static_cast<int32>(Line%Size[0]);
            const int32 Y=Min[1]+static_cast<int32>(Line/Size[0]);
            FDetailNoiseLineContext NoiseLine(static_cast<float>(X),static_cast<float>(Y),
                static_cast<uint32>(Owner->Seed));
            FDetailWorkerProfile& Profile=WorkerProfile[static_cast<size_t>(Worker)];
            if(bProfile)++Profile.Lines;
            uint64 LocalCandidates=0,LocalProxies=0,LocalFlips=0;
            const size_t WordBase=Line*WordsPerZ;
            for(size_t W=0;W<WordsPerZ;++W)
            {
                uint64 Bits=CandidateWords[WordBase+W];
                uint64 ValidBits=0;
                while(Bits)
                {
                    const uint32 Bit=std::countr_zero(Bits);Bits&=Bits-1;
                    const size_t Zi=W*64+Bit;if(Zi>=static_cast<size_t>(Size[2]))continue;
                    const int32 Z=Min[2]+static_cast<int32>(Zi);
                    ++LocalCandidates;
                    const bool bTimeDistance=bProfile&&((Profile.Candidates++&255u)==0);
                    const auto DistanceStart=bTimeDistance?Clock::now():Clock::time_point{};
                    ESurfaceFamily SurfaceFamily=ESurfaceFamily::Strata;float SurfaceSDF=100000.0f;
                    bool bFloorSurface=false,bCriticalAir=false;float CaveInfluence=0.0f,FloorInfluence=0.0f;
                    const float SDF=EffectiveRockDistance(X,Y,Z,&SurfaceFamily,&SurfaceSDF,
                        &bFloorSurface,&CaveInfluence,&FloorInfluence,&bCriticalAir);
                    if(bTimeDistance)
                    {
                        Profile.DistanceSampleNanoseconds+=static_cast<uint64>(
                            std::chrono::duration_cast<std::chrono::nanoseconds>(Clock::now()-DistanceStart).count());
                        ++Profile.DistanceSamples;
                    }
                    if(std::abs(SDF)>DetailSupremum+0.125f)continue;
                    float DetailFade=DetailBandFade(SDF);
                    const FLayerField* DetailField=nullptr;
                    bool bDetailFieldKnown=false;
                    if(DetailVariant==EDetailVariant::Coherent
                        &&std::abs(SDF)>DetailBandInnerRadius)
                    {
                        const float MacroValue=NoiseLine.Macro.Sample(static_cast<float>(Z));
                        DetailField=FindLayerField(Z);
                        bDetailFieldKnown=true;
                        DetailFade=DetailBandFade(SDF,LocalDetailBandEdge(MacroValue,SurfaceSDF,
                            CaveInfluence,FloorInfluence,DetailField));
                    }
                    if(DetailFade<=0.0f&&bProfile)++Profile.LocalBandSkips;
                    const bool bTimeDetail=DetailFade>0.0f&&bProfile
                        &&((Profile.DetailEvaluations&63u)==0);
                    const auto DetailStart=bTimeDetail?Clock::now():Clock::time_point{};
                    float Displaced=SDF;
                    if(DetailFade>0.0f)
                    {
                        if(bProfile)++Profile.DetailEvaluations;
                        Displaced+=DetailNoiseOnLine(X,Y,Z,SurfaceSDF,CaveInfluence,
                            FloorInfluence,NoiseLine,bTimeDetail?&Profile:nullptr,
                            DetailField,bDetailFieldKnown)*DetailFade;
                    }
                    if(bTimeDetail)
                    {
                        Profile.DetailSampleNanoseconds+=static_cast<uint64>(
                            std::chrono::duration_cast<std::chrono::nanoseconds>(Clock::now()-DetailStart).count());
                        ++Profile.DetailSamples;
                    }
                    if(bCriticalAir)Displaced=std::min(Displaced,-1.0f/16.0f);
                    const bool bNewRock=Displaced>=0.0f;
                    const bool bOldRock=((RockWords[WordBase+(Zi>>6)]>>(Zi&63u))&1u)!=0;
                    if(bNewRock!=bOldRock)
                    {
                        RockWords[WordBase+(Zi>>6)]^=uint64{1}<<(Zi&63u);
                        ++LocalFlips;
                        if(bProfile)++Profile.Flips;
                    }
                    const size_t Node=(Zi*static_cast<size_t>(Size[1])+static_cast<size_t>(Y-Min[1]))
                        *Size[0]+static_cast<size_t>(X-Min[0]);
                    StoreSurfaceDistance(Node,-Displaced);
                    ValidBits|=uint64{1}<<Bit;
                    ++LocalProxies;
                }
                CandidateWords[WordBase+W]=ValidBits;
            }
            CandidateCount.fetch_add(LocalCandidates,std::memory_order_relaxed);
            ProxyCount.fetch_add(LocalProxies,std::memory_order_relaxed);
            FlipCount.fetch_add(LocalFlips,std::memory_order_relaxed);
        },&Timing,&WorkerActiveSeconds);
        Stats.SurfaceCandidates=CandidateCount.load(std::memory_order_relaxed);
        Stats.SurfaceProxyNodes=ProxyCount.load(std::memory_order_relaxed);
        Stats.DetailFlips=FlipCount.load(std::memory_order_relaxed);
        if(bProfile)
        {
            uint64 DistanceCalls=0,DetailCalls=0,DistanceSamples=0,DistanceSampleNs=0;
            uint64 DetailSamples=0,DetailSampleNs=0,ReliefSampleNs=0,LocalBandSkips=0;
            uint64 LayerLookupSampleNs=0,OwnerSampleNs=0,OwnerFieldCalls=0,OwnerActiveCalls=0;
            uint64 OwnerCacheHits=0,OwnerCacheMisses=0;
            for(int32 Worker=0;Worker<WorkerCount;++Worker)
            {
                const FDetailWorkerProfile& P=WorkerProfile[static_cast<size_t>(Worker)];
                DistanceCalls+=P.Candidates;DetailCalls+=P.DetailEvaluations;
                DistanceSamples+=P.DistanceSamples;DistanceSampleNs+=P.DistanceSampleNanoseconds;
                DetailSamples+=P.DetailSamples;DetailSampleNs+=P.DetailSampleNanoseconds;
                LocalBandSkips+=P.LocalBandSkips;
                ReliefSampleNs+=P.ReliefSampleNanoseconds;
                LayerLookupSampleNs+=P.LayerLookupSampleNanoseconds;
                OwnerSampleNs+=P.OwnerRoughnessSampleNanoseconds;
                OwnerFieldCalls+=P.OwnerRoughnessFieldCalls;
                OwnerActiveCalls+=P.OwnerRoughnessActiveCalls;
                OwnerCacheHits+=P.OwnerNoiseCacheHits;
                OwnerCacheMisses+=P.OwnerNoiseCacheMisses;
            }
            const double DistanceAverageNs=DistanceSamples?static_cast<double>(DistanceSampleNs)/DistanceSamples:0.0;
            const double DetailAverageNs=DetailSamples?static_cast<double>(DetailSampleNs)/DetailSamples:0.0;
            const double CandidateWorkerNs=Stats.SurfaceCandidates
                ?Timing.WorkerSeconds*1.0e9/static_cast<double>(Stats.SurfaceCandidates):0.0;
            std::printf("[DetailBandProfile] candidates=%llu distance_evaluations=%llu detail_noise_evaluations=%llu proxies=%llu flips=%llu candidate_worker_ns=%.1f distance_sample_ns=%.1f distance_estimated_worker_ms=%.1f detail_sample_ns=%.1f detail_estimated_worker_ms=%.1f worker_sum_ms=%.3f wall_ms=%.3f sampling=(distance_1_of_256,detail_1_of_64)\n",
                static_cast<unsigned long long>(Stats.SurfaceCandidates),
                static_cast<unsigned long long>(DistanceCalls),static_cast<unsigned long long>(DetailCalls),
                static_cast<unsigned long long>(Stats.SurfaceProxyNodes),static_cast<unsigned long long>(Stats.DetailFlips),
                CandidateWorkerNs,DistanceAverageNs,DistanceAverageNs*static_cast<double>(DistanceCalls)/1.0e6,
                DetailAverageNs,DetailAverageNs*static_cast<double>(DetailCalls)/1.0e6,
                Timing.WorkerSeconds*1000.0,Timing.WallSeconds*1000.0);
            std::printf("[DetailBandSubprofile] relief_sample_ns=%.1f relief_estimated_worker_ms=%.1f layer_lookup_sample_ns=%.1f layer_lookup_estimated_worker_ms=%.1f owner_sample_ns=%.1f owner_estimated_worker_ms=%.1f owner_sampled_field_calls=%llu owner_sampled_noise_calls=%llu owner_sampled_core_cache_hits=%llu owner_sampled_core_cache_misses=%llu local_band_skipped=%llu\n",
                DetailSamples?static_cast<double>(ReliefSampleNs)/DetailSamples:0.0,
                DetailSamples?static_cast<double>(ReliefSampleNs)/DetailSamples*DetailCalls/1.0e6:0.0,
                DetailSamples?static_cast<double>(LayerLookupSampleNs)/DetailSamples:0.0,
                DetailSamples?static_cast<double>(LayerLookupSampleNs)/DetailSamples*DetailCalls/1.0e6:0.0,
                DetailSamples?static_cast<double>(OwnerSampleNs)/DetailSamples:0.0,
                DetailSamples?static_cast<double>(OwnerSampleNs)/DetailSamples*DetailCalls/1.0e6:0.0,
                static_cast<unsigned long long>(OwnerFieldCalls),
                static_cast<unsigned long long>(OwnerActiveCalls),
                static_cast<unsigned long long>(OwnerCacheHits),
                static_cast<unsigned long long>(OwnerCacheMisses),
                static_cast<unsigned long long>(LocalBandSkips));
            for(int32 Worker=0;Worker<WorkerCount;++Worker)
            {
                const FDetailWorkerProfile& P=WorkerProfile[static_cast<size_t>(Worker)];
                std::printf("[DetailBandWorker] id=%d active_ms=%.3f lines=%llu candidates=%llu detail_evaluations=%llu flips=%llu\n",
                    Worker,WorkerActiveSeconds[static_cast<size_t>(Worker)]*1000.0,
                    static_cast<unsigned long long>(P.Lines),static_cast<unsigned long long>(P.Candidates),
                    static_cast<unsigned long long>(P.DetailEvaluations),static_cast<unsigned long long>(P.Flips));
            }
        }
    }

    void StoreSurfaceDistance(size_t Node,float Density)
    {
        int32 Q=static_cast<int32>(std::lround(Density*16.0f));
        Q=std::clamp(Q,-127,127);
        if(Q==0)Q=Density<0.0f?-1:1;
        SurfaceDistanceQ[Node]=static_cast<int8>(Q);
    }

    bool IsInsideSolidOverlay(int32 X,int32 Y,int32 Z) const
    {
        const FBinKey K=PointBin(X,Y,Z);const auto It=SpatialBins.find(K);
        if(It==SpatialBins.end())return false;
        for(uint32 Id:It->second)
        {
            const FPrimitive& P=Primitives[Id];if(P.bCarve)continue;
            if(X<P.Min[0]||X>P.Max[0]||Y<P.Min[1]||Y>P.Max[1]||Z<P.Min[2]||Z>P.Max[2])continue;
            if(ShapeContains(P,X,Y,Z))return true;
        }
        return false;
    }

    bool IsInsideCriticalAir(int32 X,int32 Y,int32 Z) const
    {
        const auto It=CriticalAirBins.find(PointBin(X,Y,Z));
        if(It==CriticalAirBins.end())return false;
        for(uint32 Id:It->second)
        {
            const FPrimitive& P=Primitives[Id];if(!P.bCarve||!P.bCritical)continue;
            if(X<P.Min[0]||X>P.Max[0]||Y<P.Min[1]||Y>P.Max[1]||Z<P.Min[2]||Z>P.Max[2])continue;
            if(ShapeContains(P,X,Y,Z))return true;
        }
        return false;
    }

    bool IsInsideSolidOverlay(double X,double Y,double Z) const
    {
        const auto It=SpatialBins.find(PointBin(X,Y,Z));
        if(It==SpatialBins.end())return false;
        for(uint32 Id:It->second)
        {
            const FPrimitive& P=Primitives[Id];if(P.bCarve)continue;
            if(X<P.Min[0]||X>P.Max[0]||Y<P.Min[1]||Y>P.Max[1]||Z<P.Min[2]||Z>P.Max[2])continue;
            if(ShapeDistance(P,X,Y,Z)<=0.0f)return true;
        }
        return false;
    }

    float SolidOverlayDistance(double X,double Y,double Z) const
    {
        const auto It=SpatialBins.find(PointBin(X,Y,Z));
        if(It==SpatialBins.end())return -100000.0f;
        float Best=-100000.0f;
        for(uint32 Id:It->second)
        {
            const FPrimitive& P=Primitives[Id];if(P.bCarve)continue;
            if(X<P.Min[0]||X>P.Max[0]||Y<P.Min[1]||Y>P.Max[1]||Z<P.Min[2]||Z>P.Max[2])continue;
            Best=std::max(Best,-ShapeDistance(P,X,Y,Z));
        }
        return Best;
    }

    bool IsInsideCriticalAir(double X,double Y,double Z) const
    {
        const auto It=CriticalAirBins.find(PointBin(X,Y,Z));
        if(It==CriticalAirBins.end())return false;
        for(uint32 Id:It->second)
        {
            const FPrimitive& P=Primitives[Id];if(!P.bCarve||!P.bCritical)continue;
            if(X<P.Min[0]||X>P.Max[0]||Y<P.Min[1]||Y>P.Max[1]||Z<P.Min[2]||Z>P.Max[2])continue;
            if(ShapeContains(P,X,Y,Z))return true;
        }
        return false;
    }

    float CaveDistance(double X,double Y,double Z,bool* OutFloorSurface=nullptr,
                       float* OutFloorInfluence=nullptr) const
    {
        const auto It=SpatialBins.find(PointBin(X,Y,Z));
        if(It==SpatialBins.end())
        {
            if(OutFloorSurface)*OutFloorSurface=false;
            if(OutFloorInfluence)*OutFloorInfluence=0.0f;
            return 100000.0f;
        }
        float Best=100000.0f;
        float BestFloorInfluence=0.0f;
        for(uint32 Id:It->second)
        {
            const FPrimitive& P=Primitives[Id];if(!P.bCarve)continue;
            if(X<P.Min[0]-DetailBand||X>P.Max[0]+DetailBand||Y<P.Min[1]-DetailBand||Y>P.Max[1]+DetailBand
                ||Z<P.Min[2]-DetailBand||Z>P.Max[2]+DetailBand)continue;
            float Distance=ShapeDistance(P,X,Y,Z);
            float PrimitiveFloorInfluence=0.0f;
            if(P.bHasFloorClip)
            {
                const double Floor=FloorAt(P,X,Y);
                const float FloorT=(static_cast<float>(std::abs(Z-Floor))-static_cast<float>(P.FloorBlendRadius))
                    /static_cast<float>(std::max(DetailBand*2,1));
                PrimitiveFloorInfluence=1.0f-SmoothUnit(FloorT);
                Distance=SmoothMaxDistance(Distance,static_cast<float>(Floor-Z),
                    static_cast<float>(P.FloorBlendRadius));
            }
            const float MetadataWidth=static_cast<float>(std::max(P.BlendRadius,6.0));
            const float PrimitiveWeight=SmoothUnit(0.5f+(Best-Distance)/(2.0f*MetadataWidth));
            const float Combined=SmoothMinDistance(Best,Distance,static_cast<float>(P.BlendRadius));
            BestFloorInfluence=std::lerp(BestFloorInfluence,PrimitiveFloorInfluence,PrimitiveWeight);
            Best=Combined;
        }
        if(OutFloorInfluence)*OutFloorInfluence=BestFloorInfluence;
        if(OutFloorSurface)*OutFloorSurface=BestFloorInfluence>=0.5f;
        return Best;
    }

    float StrataDistance(double Z) const
    {
        if(LayerRanges.empty())return -100000.0f;
        return static_cast<float>(std::min(Z-LayerRanges[0].Lo,LayerRanges[0].Hi-Z));
    }

    const FLayerField* FindLayerField(double Z) const
    {
        for(const FLayerField& Field:LayerFields)
            if(Z>=Field.Lo&&Z<=Field.Hi)return &Field;
        return nullptr;
    }

    float EffectiveRockDistance(double X,double Y,double Z,
                                ESurfaceFamily* OutFamily=nullptr,
                                float* OutSurfaceSDF=nullptr,
                                bool* OutFloorSurface=nullptr,
                                float* OutCaveInfluence=nullptr,
                                float* OutFloorInfluence=nullptr,
                                bool* OutCriticalAir=nullptr) const
    {
        const float Base=StrataDistance(Z);
        float Cave=100000.0f;
        float Overlay=-100000.0f;
        float CaveFloorInfluence=0.0f;
        bool bCriticalAir=false;
        const auto Bin=SpatialBins.find(PointBin(X,Y,Z));
        if(Bin!=SpatialBins.end())
        {
            for(uint32 Id:Bin->second)
            {
                const FPrimitive& P=Primitives[Id];
                if(P.bCarve)
                {
                    if(P.bCritical&&X>=P.Min[0]&&X<=P.Max[0]&&Y>=P.Min[1]&&Y<=P.Max[1]
                        &&Z>=P.Min[2]&&Z<=P.Max[2]&&ShapeContains(P,
                            static_cast<int32>(X),static_cast<int32>(Y),static_cast<int32>(Z)))
                        bCriticalAir=true;
                    if(X<P.Min[0]-DetailBand||X>P.Max[0]+DetailBand
                        ||Y<P.Min[1]-DetailBand||Y>P.Max[1]+DetailBand
                        ||Z<P.Min[2]-DetailBand||Z>P.Max[2]+DetailBand)continue;
                    float Distance=ShapeDistance(P,X,Y,Z);
                    float PrimitiveFloorInfluence=0.0f;
                    if(P.bHasFloorClip)
                    {
                        const double Floor=FloorAt(P,X,Y);
                        const float FloorT=(static_cast<float>(std::abs(Z-Floor))
                            -static_cast<float>(P.FloorBlendRadius))
                            /static_cast<float>(std::max(DetailBand*2,1));
                        PrimitiveFloorInfluence=1.0f-SmoothUnit(FloorT);
                        Distance=SmoothMaxDistance(Distance,static_cast<float>(Floor-Z),
                            static_cast<float>(P.FloorBlendRadius));
                    }
                    const float MetadataWidth=static_cast<float>(std::max(P.BlendRadius,6.0));
                    const float PrimitiveWeight=SmoothUnit(0.5f+(Cave-Distance)/(2.0f*MetadataWidth));
                    const float Combined=SmoothMinDistance(Cave,Distance,static_cast<float>(P.BlendRadius));
                    CaveFloorInfluence=std::lerp(CaveFloorInfluence,PrimitiveFloorInfluence,PrimitiveWeight);
                    Cave=Combined;
                }
                else if(X>=P.Min[0]&&X<=P.Max[0]&&Y>=P.Min[1]&&Y<=P.Max[1]
                    &&Z>=P.Min[2]&&Z<=P.Max[2])
                    Overlay=std::max(Overlay,-ShapeDistance(P,X,Y,Z));
            }
        }
        if(OutCriticalAir)*OutCriticalAir=bCriticalAir;
        const float CarveOrOverlay=std::max(Cave,Overlay);
        const ESurfaceFamily Family=Base<=CarveOrOverlay?ESurfaceFamily::Strata:
            (Cave>=Overlay?ESurfaceFamily::Cave:ESurfaceFamily::SolidOverlay);
        constexpr float FamilyBlendWidth=6.0f;
        const float CaveVsOverlay=SmoothUnit(0.5f+(Cave-Overlay)/(2.0f*FamilyBlendWidth));
        const float CarveInfluence=SmoothUnit(0.5f+(Base-CarveOrOverlay)/(2.0f*FamilyBlendWidth));
        const float CaveInfluence=CarveInfluence*CaveVsOverlay;
        const float FloorInfluence=CaveInfluence*CaveFloorInfluence;
        const float RockDistance=SmoothMinDistance(Base,CarveOrOverlay,StrataCaveBlendRadius);
        if(OutFamily)
            *OutFamily=Family;
        if(OutSurfaceSDF)
            *OutSurfaceSDF=Family==ESurfaceFamily::Strata?Base:
                (Family==ESurfaceFamily::Cave?Cave:Overlay);
        if(OutCaveInfluence)*OutCaveInfluence=CaveInfluence;
        if(OutFloorInfluence)*OutFloorInfluence=FloorInfluence;
        if(OutFloorSurface)*OutFloorSurface=FloorInfluence>=0.5f;
        return RockDistance;
    }

    float CoarseRockDistance(double X,double Y,double Z,bool* OutCriticalAir) const
    {
        const float Base=StrataDistance(Z);
        float Cave=100000.0f,Overlay=-100000.0f;
        bool bCriticalAir=false;
        const auto Bin=SpatialBins.find(PointBin(X,Y,Z));
        if(Bin!=SpatialBins.end())
        {
            for(uint32 Id:Bin->second)
            {
                const FPrimitive& P=Primitives[Id];
                if(P.bCarve)
                {
                    if(P.bCritical&&X>=P.Min[0]&&X<=P.Max[0]&&Y>=P.Min[1]&&Y<=P.Max[1]
                        &&Z>=P.Min[2]&&Z<=P.Max[2]&&ShapeContains(P,
                            static_cast<int32>(X),static_cast<int32>(Y),static_cast<int32>(Z)))
                        bCriticalAir=true;
                    if(X<P.Min[0]-DetailBand||X>P.Max[0]+DetailBand
                        ||Y<P.Min[1]-DetailBand||Y>P.Max[1]+DetailBand
                        ||Z<P.Min[2]-DetailBand||Z>P.Max[2]+DetailBand)continue;
                    float Distance=ShapeDistance(P,X,Y,Z);
                    if(P.bHasFloorClip)
                    {
                        const double Floor=FloorAt(P,X,Y);
                        Distance=SmoothMaxDistance(Distance,static_cast<float>(Floor-Z),
                            static_cast<float>(P.FloorBlendRadius));
                    }
                    Cave=SmoothMinDistance(Cave,Distance,static_cast<float>(P.BlendRadius));
                }
                else if(X>=P.Min[0]&&X<=P.Max[0]&&Y>=P.Min[1]&&Y<=P.Max[1]
                    &&Z>=P.Min[2]&&Z<=P.Max[2])
                    Overlay=std::max(Overlay,-ShapeDistance(P,X,Y,Z));
            }
        }
        if(OutCriticalAir)*OutCriticalAir=bCriticalAir;
        return SmoothMinDistance(Base,std::max(Cave,Overlay),StrataCaveBlendRadius);
    }

    struct FOwnerNoiseLineContext
    {
        const FLayerField* Field=nullptr;
        int32 GridZ0=MIN_int32;
        float LowerCore=0.0f,UpperCore=0.0f;
    };

    float OwnerCaveRoughnessCore(double X,double Y,double Z,
                                 const FStrateGenerationParams& Params) const
    {
        const uint32 Seed=static_cast<uint32>(Owner->Seed);
        const float RF=Params.RoughnessFrequency;
        const float EffectiveZ=(Params.VerticalScale>0.0f&&Params.VerticalScale!=1.0f)
            ?static_cast<float>(Z)/Params.VerticalScale:static_cast<float>(Z);
        FVector3f Main(static_cast<float>(X)*RF+VoxelHash::SeedOffset(Seed,11.3f),
            static_cast<float>(Y)*RF+VoxelHash::SeedOffset(Seed,13.7f),
            EffectiveZ*RF+VoxelHash::SeedOffset(Seed,17.1f));
        FVector3f Fine(static_cast<float>(X)*RF*3.0f+VoxelHash::SeedOffset(Seed,19.1f)+2000.0f,
            static_cast<float>(Y)*RF*3.0f+VoxelHash::SeedOffset(Seed,23.7f)+2500.0f,
            EffectiveZ*RF*3.0f+VoxelHash::SeedOffset(Seed,29.3f)+3000.0f);
        if(Params.DomainWarpStrength>0.0f)
        {
            const float WF=Params.DomainWarpFrequency,WS=Params.DomainWarpStrength;
            const float WX=VoxelNoise::Perlin3D(FVector3f(
                static_cast<float>(X)*WF+VoxelHash::SeedOffset(Seed,5.2f),
                static_cast<float>(Y)*WF+VoxelHash::SeedOffset(Seed,1.3f),
                EffectiveZ*WF+VoxelHash::SeedOffset(Seed,9.7f)))*VOXEL_NOISE_SCALE*WS;
            const float WY=VoxelNoise::Perlin3D(FVector3f(
                static_cast<float>(X)*WF+100.0f+VoxelHash::SeedOffset(Seed,7.7f),
                static_cast<float>(Y)*WF+200.0f+VoxelHash::SeedOffset(Seed,3.1f),
                EffectiveZ*WF+300.0f))*VOXEL_NOISE_SCALE*WS;
            const float WZ=VoxelNoise::Perlin3D(FVector3f(
                static_cast<float>(X)*WF+400.0f,
                static_cast<float>(Y)*WF+500.0f+VoxelHash::SeedOffset(Seed,11.9f),
                EffectiveZ*WF+600.0f+VoxelHash::SeedOffset(Seed,13.3f)))*VOXEL_NOISE_SCALE*WS;
            Main+=FVector3f(WX,WY,WZ);Fine+=FVector3f(WX,WY,WZ);
        }
        float Rough=0.0f,FineValue=0.0f;
        switch(Params.RoughnessNoiseType)
        {
        case EVoxelNoiseType::Ridged:
            Rough=VoxelNoise::Ridged(Main.X,Main.Y,Main.Z,3);
            FineValue=VoxelNoise::Ridged(Fine.X,Fine.Y,Fine.Z,2);break;
        case EVoxelNoiseType::Mixed:
            Rough=0.5f*(VoxelNoise::FBM(Main.X,Main.Y,Main.Z,3)
                +VoxelNoise::Ridged(Main.X,Main.Y,Main.Z,3));
            FineValue=0.5f*(VoxelNoise::FBM(Fine.X,Fine.Y,Fine.Z,2)
                +VoxelNoise::Ridged(Fine.X,Fine.Y,Fine.Z,2));break;
        case EVoxelNoiseType::Cellular:
            Rough=VoxelNoise::Cellular3D(Main);FineValue=VoxelNoise::Cellular3D(Fine);break;
        case EVoxelNoiseType::FBM:
        default:
            Rough=VoxelNoise::FBM(Main.X,Main.Y,Main.Z,3);
            FineValue=VoxelNoise::FBM(Fine.X,Fine.Y,Fine.Z,2);break;
        }
        return (Rough+FineValue*0.4f)*VOXEL_NOISE_SCALE*Params.SurfaceRoughness;
    }

    float OwnerCaveRoughness(double X,double Y,double Z,float SurfaceSDF,
                             const FLayerField& Field,float CaveInfluence,
                             FOwnerNoiseLineContext& LineContext,
                             FDetailWorkerProfile* Profile=nullptr) const
    {
        const FStrateGenerationParams& Params=Field.Params;
        const float Depth=Params.SurfaceRoughness*2.0f;
        if(!(Depth>0.0f)||std::abs(SurfaceSDF)>=Depth||!(Params.RoughnessFrequency>0.0f))return 0.0f;
        if(Profile)++Profile->OwnerRoughnessActiveCalls;
        constexpr int32 ZStep=8;
        const int32 GridZ0=static_cast<int32>(std::floor(Z/static_cast<double>(ZStep)))*ZStep;
        const int32 SampleZ0=std::max(GridZ0,Field.Lo);
        const int32 SampleZ1=std::min(GridZ0+ZStep,Field.Hi);
        float Total=0.0f;
        if(SampleZ1>SampleZ0)
        {
            const bool bCacheHit=LineContext.Field==&Field&&LineContext.GridZ0==GridZ0;
            if(Profile)
            {
                if(bCacheHit)++Profile->OwnerNoiseCacheHits;
                else ++Profile->OwnerNoiseCacheMisses;
            }
            if(!bCacheHit)
            {
                LineContext.Field=&Field;
                LineContext.GridZ0=GridZ0;
                LineContext.LowerCore=OwnerCaveRoughnessCore(X,Y,SampleZ0,Params);
                LineContext.UpperCore=OwnerCaveRoughnessCore(X,Y,SampleZ1,Params);
            }
            const float Alpha=static_cast<float>((Z-SampleZ0)/(SampleZ1-SampleZ0));
            Total=FMath::Lerp(LineContext.LowerCore,LineContext.UpperCore,Alpha);
        }
        else
            Total=OwnerCaveRoughnessCore(X,Y,Z,Params);
        if(SurfaceSDF<0.0f)
            Total=std::lerp(Total,std::min(Total,0.0f),std::clamp(CaveInfluence,0.0f,1.0f));
        const float Fade=1.0f-std::abs(SurfaceSDF)/Depth;
        return Total*Fade*Fade;
    }

    static float SurfaceRelief(float X,float Y,float Z,uint32 Seed)
    {
        // The low band carries the broad rock shapes and also gently varies mid-band
        // amplitude. The quieter 6-voxel octave keeps fine relief without restoring
        // the same tight bubble pattern on every exposed surface.
        const float Macro=FStandaloneResearchTunnels::ValueNoise(X,Y,Z,32,Seed^0xB5297A4Du);
        const float Mid=FStandaloneResearchTunnels::ValueNoise(X,Y,Z,16,Seed^0x63D83595u);
        const float Fine=FStandaloneResearchTunnels::ValueNoise(X,Y,Z,6,Seed^0xA511E9B3u);
        return 1.00f*Macro+(0.80f+0.25f*Macro)*(0.85f*Mid+0.18f*Fine);
    }

    struct FValueNoiseLineContext
    {
        int32 X0=0,Y0=0,Z0=std::numeric_limits<int32>::min();
        int32 Wavelength=1;
        uint32 Seed=0;
        float Inv=1.0f,TX=0.0f,TY=0.0f;
        std::array<float,8> Corners{};

        FValueNoiseLineContext(float X,float Y,int32 InWavelength,uint32 InSeed)
            :Wavelength(InWavelength),Seed(InSeed),Inv(1.0f/static_cast<float>(InWavelength))
        {
            const float FX=X*Inv,FY=Y*Inv;
            X0=FMath::FloorToInt(FX);Y0=FMath::FloorToInt(FY);
            TX=FStandaloneResearchTunnels::Fade(FX-static_cast<float>(X0));
            TY=FStandaloneResearchTunnels::Fade(FY-static_cast<float>(Y0));
        }

        void SetPlane(int32 Z,bool bUpper)
        {
            const size_t I=bUpper?4u:0u;
            Corners[I+0]=FStandaloneResearchTunnels::NoiseValue(X0,Y0,Z,Seed);
            Corners[I+1]=FStandaloneResearchTunnels::NoiseValue(X0+1,Y0,Z,Seed);
            Corners[I+2]=FStandaloneResearchTunnels::NoiseValue(X0,Y0+1,Z,Seed);
            Corners[I+3]=FStandaloneResearchTunnels::NoiseValue(X0+1,Y0+1,Z,Seed);
        }

        float Sample(float Z)
        {
            const float FZ=Z*Inv;
            const int32 NewZ0=FMath::FloorToInt(FZ);
            const float TZ=FStandaloneResearchTunnels::Fade(FZ-static_cast<float>(NewZ0));
            if(NewZ0!=Z0)
            {
                if(NewZ0==Z0+1)
                {
                    Corners[0]=Corners[4];Corners[1]=Corners[5];
                    Corners[2]=Corners[6];Corners[3]=Corners[7];
                    SetPlane(NewZ0+1,true);
                }
                else
                {
                    SetPlane(NewZ0,false);SetPlane(NewZ0+1,true);
                }
                Z0=NewZ0;
            }
            const float Z00=FMath::Lerp(Corners[0],Corners[1],TX);
            const float Z10=FMath::Lerp(Corners[2],Corners[3],TX);
            const float Z01=FMath::Lerp(Corners[4],Corners[5],TX);
            const float Z11=FMath::Lerp(Corners[6],Corners[7],TX);
            return FMath::Lerp(FMath::Lerp(Z00,Z10,TY),FMath::Lerp(Z01,Z11,TY),TZ);
        }
    };

    struct FDetailNoiseLineContext
    {
        FValueNoiseLineContext Macro,Mid,Fine;
        FOwnerNoiseLineContext OwnerNoise;
        FDetailNoiseLineContext(float X,float Y,uint32 Seed)
            :Macro(X,Y,32,Seed^0xB5297A4Du),Mid(X,Y,16,Seed^0x63D83595u),
             Fine(X,Y,6,Seed^0xA511E9B3u)
        {}

        float SurfaceReliefAt(float Z)
        {
            const float MacroValue=Macro.Sample(Z);
            const float MidValue=Mid.Sample(Z);
            const float FineValue=Fine.Sample(Z);
            return 1.00f*MacroValue+(0.80f+0.25f*MacroValue)
                *(0.85f*MidValue+0.18f*FineValue);
        }
    };

    float DetailNoise(double X,double Y,double Z,float SurfaceSDF,
                      float CaveInfluence,float FloorInfluence) const
    {
        if(DetailVariant==EDetailVariant::Smooth)return 0.0f;
        if(DetailVariant==EDetailVariant::PerVoxel)
        {
            const int32 IX=static_cast<int32>(std::floor(X)),IY=static_cast<int32>(std::floor(Y)),IZ=static_cast<int32>(std::floor(Z));
            return FStandaloneResearchTunnels::NoiseValue(IX,IY,IZ,static_cast<uint32>(Owner->Seed))*DetailSupremum;
        }
        const float Relief=SurfaceRelief(static_cast<float>(X),
            static_cast<float>(Y),static_cast<float>(Z),static_cast<uint32>(Owner->Seed));
        const FLayerField* Field=FindLayerField(Z);
        FOwnerNoiseLineContext OwnerNoise;
        const float OwnerRelief=Field?OwnerCaveRoughness(X,Y,Z,SurfaceSDF,*Field,
            CaveInfluence,OwnerNoise):0.0f;
        const float Contrast=1.35f-0.80f*std::clamp(CaveInfluence,0.0f,1.0f)
            +1.85f*std::clamp(FloorInfluence,0.0f,1.0f);
        return std::clamp((Relief+OwnerRelief*0.5f)*Contrast,-3.9f,3.9f);
    }

    float DetailNoiseOnLine(double X,double Y,double Z,float SurfaceSDF,
                            float CaveInfluence,float FloorInfluence,
                            FDetailNoiseLineContext& NoiseLine,
                            FDetailWorkerProfile* Profile=nullptr,
                            const FLayerField* KnownField=nullptr,
                            bool bFieldKnown=false) const
    {
        if(DetailVariant!=EDetailVariant::Coherent)
            return DetailNoise(X,Y,Z,SurfaceSDF,CaveInfluence,FloorInfluence);
        using Clock=std::chrono::steady_clock;
        const auto ReliefStart=Profile?Clock::now():Clock::time_point{};
        const float Relief=NoiseLine.SurfaceReliefAt(static_cast<float>(Z));
        if(Profile)
            Profile->ReliefSampleNanoseconds+=static_cast<uint64>(
                std::chrono::duration_cast<std::chrono::nanoseconds>(Clock::now()-ReliefStart).count());
        const auto LookupStart=Profile?Clock::now():Clock::time_point{};
        const FLayerField* Field=bFieldKnown?KnownField:FindLayerField(Z);
        if(Profile)
            Profile->LayerLookupSampleNanoseconds+=static_cast<uint64>(
                std::chrono::duration_cast<std::chrono::nanoseconds>(Clock::now()-LookupStart).count());
        const auto OwnerStart=Profile?Clock::now():Clock::time_point{};
        if(Profile&&Field)++Profile->OwnerRoughnessFieldCalls;
        const float OwnerRelief=Field?OwnerCaveRoughness(X,Y,Z,SurfaceSDF,*Field,
            CaveInfluence,NoiseLine.OwnerNoise,Profile):0.0f;
        if(Profile)
            Profile->OwnerRoughnessSampleNanoseconds+=static_cast<uint64>(
                std::chrono::duration_cast<std::chrono::nanoseconds>(Clock::now()-OwnerStart).count());
        const float Contrast=1.35f-0.80f*std::clamp(CaveInfluence,0.0f,1.0f)
            +1.85f*std::clamp(FloorInfluence,0.0f,1.0f);
        return std::clamp((Relief+OwnerRelief*0.5f)*Contrast,-3.9f,3.9f);
    }

    bool EvaluateFast(double X,double Y,double Z) const
    {
        ESurfaceFamily SurfaceFamily=ESurfaceFamily::Strata;float SurfaceSDF=100000.0f;
        bool bFloorSurface=false,bCriticalAir=false;float CaveInfluence=0.0f,FloorInfluence=0.0f;
        float RockDistance=EffectiveRockDistance(X,Y,Z,&SurfaceFamily,&SurfaceSDF,
            &bFloorSurface,&CaveInfluence,&FloorInfluence,&bCriticalAir);
        const bool bInsideFine=X>=Min[0]&&X<=Max[0]&&Y>=Min[1]&&Y<=Max[1]&&Z>=Min[2]&&Z<=Max[2];
        if(bInsideFine&&DetailVariant!=EDetailVariant::Smooth
            &&std::abs(RockDistance)<=DetailSupremum+0.125f)
        {
            const float Fade=DetailBandFadeAt(X,Y,Z,RockDistance,SurfaceSDF,
                CaveInfluence,FloorInfluence);
            if(Fade>0.0f)RockDistance+=DetailNoise(X,Y,Z,SurfaceSDF,CaveInfluence,FloorInfluence)*Fade;
        }
        if(bCriticalAir)RockDistance=std::min(RockDistance,-1.0f/16.0f);
        return RockDistance>=0.0f;
    }

    void ValidatePassageRoutes()
    {
        constexpr double SampleStep=0.5;
        for(const FVoxelPassage& P:Owner->Manager->GetPassages())
        {
            const auto& Points=P.ControlPoints;
            if(Points.Num()<2)continue;
            for(int32 I=0;I+1<Points.Num();++I)
            {
                const FVector A=Points[I],B=Points[I+1];
                const double L=std::sqrt((B.X-A.X)*(B.X-A.X)+(B.Y-A.Y)*(B.Y-A.Y)+(B.Z-A.Z)*(B.Z-A.Z));
                const int32 Count=std::max(1,static_cast<int32>(std::ceil(L/SampleStep)));
                for(int32 J=0;J<=Count;++J)
                {
                    const double T=static_cast<double>(J)/Count;
                    const double X=A.X+(B.X-A.X)*T,Y=A.Y+(B.Y-A.Y)*T,Z=A.Z+(B.Z-A.Z)*T;
                    ++Stats.WalkSamples;
                    if(IsRock(static_cast<int32>(std::lround(X)),
                        static_cast<int32>(std::lround(Y)),static_cast<int32>(std::lround(Z))))
                        ++Stats.WalkFailures;
                }
            }
        }
        if(Stats.WalkFailures)
            std::fprintf(stderr,"[Walkability] FAIL %llu/%llu passage centerline samples are solid\n",
                static_cast<unsigned long long>(Stats.WalkFailures),static_cast<unsigned long long>(Stats.WalkSamples));
        else
            std::printf("[Walkability] PASS %llu/%llu passage centerline samples remain open; native landing floors retained\n",
                static_cast<unsigned long long>(Stats.WalkSamples),static_cast<unsigned long long>(Stats.WalkSamples));
    }

    static double DetailAutocorrelation(EDetailVariant Variant)
    {
        constexpr int32 N=24;double Sum=0.0,SumSq=0.0,PairA=0.0,PairB=0.0,PairAA=0.0,PairBB=0.0,PairAB=0.0;uint64 Count=0;
        auto Value=[&](int32 X,int32 Y,int32 Z)->double
        {
            if(Variant==EDetailVariant::PerVoxel)return FStandaloneResearchTunnels::NoiseValue(X,Y,Z,11);
            if(Variant==EDetailVariant::Smooth)return 0.0;
            return SurfaceRelief(static_cast<float>(X),static_cast<float>(Y),static_cast<float>(Z),11);
        };
        for(int32 Z=0;Z<N;++Z)for(int32 Y=0;Y<N;++Y)for(int32 X=0;X<N;++X)
        {const double V=Value(X,Y,Z);Sum+=V;SumSq+=V*V;}
        const double Mean=Sum/(N*N*N),Variance=SumSq/(N*N*N)-Mean*Mean;
        if(Variance<1.0e-9)return 1.0;
        for(int32 Z=0;Z<N;++Z)for(int32 Y=0;Y<N;++Y)for(int32 X=0;X+1<N;++X)
        {
            const double A=Value(X,Y,Z)-Mean,B=Value(X+1,Y,Z)-Mean;
            PairAA+=A*A;PairBB+=B*B;PairAB+=A*B;++Count;
        }
        return Count?PairAB/std::sqrt(std::max(1.0e-12,PairAA*PairBB)):0.0;
    }

    template<class Fn> void ParallelLines(size_t Count,Fn&& Work,FParallelTiming* Timing=nullptr,
        std::array<double,16>* WorkerActiveSeconds=nullptr) const
    {
        const auto WallStart=std::chrono::steady_clock::now();
        if(WorkerCount<=1||Count<512)
        {
            const auto WorkerStart=std::chrono::steady_clock::now();
            for(size_t I=0;I<Count;++I)
            {
                if constexpr(std::is_invocable_v<Fn&,size_t,int32>)Work(I,0);
                else Work(I);
            }
            if(WorkerActiveSeconds)
                (*WorkerActiveSeconds)[0]=std::chrono::duration<double>(
                    std::chrono::steady_clock::now()-WorkerStart).count();
            if(Timing)
            {
                const double WorkerSeconds=std::chrono::duration<double>(std::chrono::steady_clock::now()-WorkerStart).count();
                Timing->WallSeconds+=std::chrono::duration<double>(std::chrono::steady_clock::now()-WallStart).count();
                Timing->WorkerSeconds+=WorkerSeconds;
            }
            return;
        }
        std::atomic<size_t> Next{0};std::atomic<uint64> WorkerNanoseconds{0};
        std::vector<std::thread> Threads;Threads.reserve(WorkerCount);
        for(int32 W=0;W<WorkerCount;++W)
            Threads.emplace_back([&,W]()
            {
                const auto WorkerStart=std::chrono::steady_clock::now();
                for(;;)
                {
                    const size_t I=Next.fetch_add(1,std::memory_order_relaxed);if(I>=Count)break;
                    if constexpr(std::is_invocable_v<Fn&,size_t,int32>)Work(I,W);
                    else Work(I);
                }
                const uint64 WorkerNs=static_cast<uint64>(std::chrono::duration_cast<std::chrono::nanoseconds>(
                    std::chrono::steady_clock::now()-WorkerStart).count());
                if(WorkerActiveSeconds)(*WorkerActiveSeconds)[static_cast<size_t>(W)]=static_cast<double>(WorkerNs)*1.0e-9;
                WorkerNanoseconds.fetch_add(WorkerNs,std::memory_order_relaxed);
            });
        for(std::thread& T:Threads)T.join();
        if(Timing)
        {
            Timing->WallSeconds+=std::chrono::duration<double>(std::chrono::steady_clock::now()-WallStart).count();
            Timing->WorkerSeconds+=static_cast<double>(WorkerNanoseconds.load(std::memory_order_relaxed))*1.0e-9;
        }
    }

    FStandaloneOwnerWorld* Owner=nullptr;
    int32 WorkerCount=6;
    int32 Min[3]{},Max[3]{},ClipMin[3]{},ClipMax[3]{},Size[3]{};
    FIntVector EditBoundsMin=FIntVector(MAX_int32,MAX_int32,MAX_int32);
    FIntVector EditBoundsMax=FIntVector(MIN_int32,MIN_int32,MIN_int32);
    int32 ActiveSpatialBinSize=LineBinSize;
    double ActiveBlendRadius=0.0;
    double ActiveFloorBlendRadius=0.0;
    size_t WordsPerZ=0;
    EDetailVariant DetailVariant=EDetailVariant::Coherent;
    bool Built=false;
    bool FarBuilt=false;
    std::vector<FInterval> LayerRanges;
    std::vector<FLayerField> LayerFields;
    std::vector<FPrimitive> Primitives;
    std::vector<FEditBrush> EditBrushes;
    std::array<FAxisBins,3> AxisBins;
    std::unordered_map<FBinKey,FBinItems,FBinHash> SpatialBins;
    std::unordered_map<FBinKey,FBinItems,FBinHash> CriticalAirBins;
    std::vector<uint64> RockWords;
    std::unique_ptr<int8[]> SurfaceDistanceQ;
    inline static std::atomic<bool> ProfileDetailBand{false};
    std::array<std::vector<uint64>,3> CandidateAxisWords;
    std::vector<uint64> CandidateWords;
    FBuildStats Stats;
};

namespace FStandaloneShapeWorldAccess
{
    inline std::atomic<const FStandaloneShapeWorld*>& ActiveRef()
    {static std::atomic<const FStandaloneShapeWorld*> Active{nullptr};return Active;}
    inline bool& ThreadBypassRef()
    {static thread_local bool Bypass=false;return Bypass;}
    inline void SetActive(const FStandaloneShapeWorld* World)
    {ActiveRef().store(World,std::memory_order_release);}
    inline const FStandaloneShapeWorld* GetActive()
    {return ThreadBypassRef()?nullptr:ActiveRef().load(std::memory_order_acquire);}
    class FScopedBypass
    {
    public:
        FScopedBypass():Previous(ThreadBypassRef()){ThreadBypassRef()=true;}
        ~FScopedBypass(){ThreadBypassRef()=Previous;}
        FScopedBypass(const FScopedBypass&)=delete;
        FScopedBypass& operator=(const FScopedBypass&)=delete;
    private:
        bool Previous;
    };
}