wip(assembler): M1 repair round -- raised block and cliff fixed; native 7x faster than legacy

Luna, M1 repair round: the rectangular raised block and the broken cliff shared one cause, a type-2 solid AABB in the generated plan; it is now a tapered, terrain-following shelf, and a negative-control audit catches the old box behaviour; seam classification check passes. Bound audit 298 boxes / 8,046 samples, 0 misses. Evaluator and 120-tile density/mesh/payload/ topology hashes identical at W1/W8/W16. Cave regression: walkability 5,642/5,642, 64/64 tiles. Legacy mismatch unchanged (0 biome mismatches, 18.14% signs). Speed (matched workloads): 120-tile full detail 0.089 s native vs 0.625 s legacy; 529-tile clipmap 3.858 s vs 11.687 s -- 2.29 s above the ~1.57 s aim; LOD8 dominates (165 unknown tiles, 13.27 M field calls, 10.88 M column-cache misses). My capture review: coast matches legacy closely; surface clean but a dark flat slab on the far horizon; cliff shelf renders as a stack of thin plates with back faces showing (sub-voxel shelf thickness suspected). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MmSoF64BbRFmk9C3p1VsBJ
2026-10-09 01:22:30 +02:00
1 parent fcb16fd265
commit 2c20f9eaea
@@ -62,6 +62,28 @@ namespace VoxelForgeCore
        std::uint32_t skirtTriangleCount = 0;
    };

+    enum class TileClass : std::uint8_t { allSolid, allAir, mixed, unknown };
+
+    struct TileGenerationStats
+    {
+        TileClass tileClass = TileClass::unknown;
+        std::uint64_t solidBoxes = 0;
+        std::uint64_t airBoxes = 0;
+        std::uint64_t mixedBoxes = 0;
+        std::uint64_t unknownBoxes = 0;
+        std::uint64_t fieldCalls = 0;
+        std::uint64_t floatFieldCalls = 0;
+        std::uint64_t latticeFieldCalls = 0;
+        std::uint64_t unclampedFieldCalls = 0;
+        std::uint64_t columnCacheHits = 0;
+        std::uint64_t columnCacheMisses = 0;
+        double classificationMilliseconds = 0.0;
+        double columnComputeMilliseconds = 0.0;
+        double meshingMilliseconds = 0.0;
+        bool tileWasClassified = false;
+        bool meshingSkipped = false;
+    };
+
    struct TileMesh
    {
        TileKey key;
@@ -72,6 +94,7 @@ namespace VoxelForgeCore
        std::vector<std::uint8_t> skirtFaceMasks;
        std::vector<float> debugDensitySamples;
        MeshGroups groups;
+        TileGenerationStats diagnostics;
        std::uint64_t densitySampleCount = 0;
        std::uint64_t meshPayloadAllocatedBytes = 0;
        std::uint64_t densityHash = 0, meshHash = 0, payloadHash = 0, topologyHash = 0;
@@ -103,7 +126,7 @@ namespace VoxelForgeCore
    // L'assembleur v1 encode coordonnees et scalaires en 1/256 de voxel, independamment des tuiles.
    inline constexpr std::int64_t WorldAssemblySubvoxelsPerVoxel = 256;
    inline constexpr std::uint32_t WorldAssemblySchemaVersion = 1;
-    inline constexpr std::uint32_t WorldAssemblyEvaluatorVersion = 2;
+    inline constexpr std::uint32_t WorldAssemblyEvaluatorVersion = 3;

    enum class FieldChannel : std::uint32_t
    {
@@ -177,9 +200,10 @@ namespace VoxelForgeCore

    struct WorldAssemblyFeatureV1
    {
-        // Stable numeric type IDs are part of the plan contract; v1 reserves 1 for cliff planes,
-        // 2 for bounded overhang shelves and 3 for authored cap features.
-        // Les types numeriques stables font partie du plan; v1 reserve 1 aux falaises, 2 aux surplombs bornes et 3 aux caps.
+        // Stable numeric type IDs are part of the plan contract: 1 is a heightfield cliff plane,
+        // 2 is the legacy bounded box exception, 3 is a cap feature, and 4 is a tapered terrain-
+        // following cliff shelf.
+        // Types stables: 1 plan de falaise, 2 ancienne boite bornee, 3 cap, 4 corniche adaptee au terrain.
        std::uint64_t stableId = 0;
        std::uint32_t type = 0;
        std::array<std::int64_t, 3> minimumSubvoxels{};
@@ -251,6 +275,17 @@ namespace VoxelForgeCore
        std::int64_t maximumDensitySubvoxels = 0;
    };

+    enum class HeightfieldBoxClass : std::uint8_t
+    { unknown, allSolid, allAir, mixed };
+
+    struct HeightfieldColumnCacheStats
+    {
+        std::uint64_t pointCalls = 0;
+        std::uint64_t columnHits = 0;
+        std::uint64_t columnMisses = 0;
+        double columnComputeMilliseconds = 0.0;
+    };
+
    struct FieldInfluenceBounds
    {
        bool unboundedXY = true;
@@ -283,6 +318,12 @@ namespace VoxelForgeCore
        bool isValid() const noexcept;
        HeightfieldPointSample point(std::int64_t xSubvoxels, std::int64_t ySubvoxels,
                                     std::int64_t zSubvoxels) const noexcept;
+        HeightfieldBoxClass classifyLatticeBox(std::int64_t minimumXSubvoxels,
+            std::int64_t minimumYSubvoxels, std::int64_t minimumZSubvoxels,
+            std::int64_t maximumXSubvoxels, std::int64_t maximumYSubvoxels,
+            std::int64_t maximumZSubvoxels, std::int32_t stepVoxels) const noexcept;
+        static void resetThreadColumnCacheStats() noexcept;
+        static HeightfieldColumnCacheStats threadColumnCacheStats() noexcept;
        FieldInterval interval(std::int64_t minimumXSubvoxels,
                               std::int64_t minimumYSubvoxels,
                               std::int64_t minimumZSubvoxels,
@@ -356,8 +397,6 @@ namespace VoxelForgeCore
    };

    enum class DetailVariant : std::uint8_t { coherent, smooth, perVoxel };
-    enum class TileClass : std::uint8_t { allSolid, allAir, mixed };
-
    struct DetailBandTileRegion
    {
        Vec3i originVoxel{};
@@ -33,6 +33,21 @@ class FStandaloneShapeWorld
 {
 public:
    enum class EDetailVariant : uint8 { Coherent, Smooth, PerVoxel };
+    struct FBoxClassification
+    {
+        EVoxelTileClass Verdict=EVoxelTileClass::Mixed;
+        bool bNativeSurface=false;
+        bool bUnknown=false;
+    };
+    struct FClassificationThreadStats
+    {
+        uint64 SolidBoxes=0,AirBoxes=0,MixedBoxes=0,UnknownBoxes=0;
+        uint64 NativeFloatFieldCalls=0,NativeLatticeFieldCalls=0,NativeUnclampedFieldCalls=0;
+        double Milliseconds=0.0;
+    };
+    static void ResetThreadClassificationStats() { ThreadClassificationStats()={}; }
+    static FClassificationThreadStats GetThreadClassificationStats()
+    { return ThreadClassificationStats(); }
    struct FFarConnectorSegment
    {
        FVector Start=FVector::ZeroVector;
@@ -923,7 +938,10 @@ public:
    {
        float NativeDensity=0.0f;
        if(NativeSurfaceDensityAt(X,Y,Z,NativeDensity))
+        {
+            ++ThreadClassificationStats().NativeFloatFieldCalls;
            return ApplyEditBrushes(NativeDensity,X,Y,Z);
+        }
        if(UsesLegacySurfaceField(X,Y,Z))
            return LegacySurfaceDensityAt(X,Y,Z);
        int32 IX=0,IY=0,IZ=0;
@@ -939,7 +957,10 @@ public:
        float NativeDensity=0.0f;
        if(NativeSurfaceDensityAt(static_cast<float>(X),static_cast<float>(Y),
            static_cast<float>(Z),NativeDensity))
+        {
+            ++ThreadClassificationStats().NativeLatticeFieldCalls;
            return ApplyEditBrushes(NativeDensity,X,Y,Z);
+        }
        if(UsesLegacySurfaceField(static_cast<float>(X),static_cast<float>(Y),
            static_cast<float>(Z)))
            return LegacySurfaceDensityAt(static_cast<float>(X),static_cast<float>(Y),
@@ -956,7 +977,10 @@ public:
    {
        float NativeDensity=0.0f;
        if(NativeSurfaceDensityAt(X,Y,Z,NativeDensity))
+        {
+            ++ThreadClassificationStats().NativeUnclampedFieldCalls;
            return ApplyEditBrushes(NativeDensity,X,Y,Z);
+        }
        if(UsesLegacySurfaceField(X,Y,Z))
            return LegacySurfaceDensityAt(X,Y,Z);
        return AnalyticDensityAt(X,Y,Z,false);
@@ -1130,12 +1154,17 @@ private:

 public:
    EVoxelTileClass ClassifyBox(const FIntVector& Origin,int32 Step,int32 Cells) const
+    { return ClassifyBoxDetailed(Origin,Step,Cells).Verdict; }
+
+    FBoxClassification ClassifyBoxDetailed(const FIntVector& Origin,int32 Step,int32 Cells) const
    {
+        using Clock=std::chrono::steady_clock;
+        const auto Start=Clock::now();
        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};
        const VoxelForgeCore::WorldAssemblyStratumPlanV1* SurfacePlan=nullptr;
-        bool bSurfaceOverlap=false,bOtherStrateOverlap=false;
+        bool bSurfaceOverlap=false,bOtherStrateOverlap=false,bBadPlan=false;
        if(Owner&&Owner->Manager)
        {
            const int32 FirstChunk=FMath::FloorToInt(static_cast<float>(BMin[2])
@@ -1152,19 +1181,31 @@ public:
                        static_cast<float>((static_cast<double>(ChunkZ*CHUNK_SIZE+CHUNK_SIZE/2))*VOXEL_SIZE));
                    const auto* Candidate=FindWorldAssemblyStratum(StrateIndex,ChunkZ);
                    if(!Candidate||(SurfacePlan&&SurfacePlan!=Candidate))
-                        return EVoxelTileClass::Mixed;
+                        bBadPlan=true;
                    SurfacePlan=Candidate;
                }
                else bOtherStrateOverlap=true;
            }
        }
+        const auto Finish=[&](EVoxelTileClass Verdict,bool bUnknown=false)
+        {
+            FBoxClassification Result{Verdict,bSurfaceOverlap&&!bOtherStrateOverlap,bUnknown};
+            auto& Counts=ThreadClassificationStats();
+            Counts.Milliseconds+=std::chrono::duration<double,std::milli>(Clock::now()-Start).count();
+            if(!Result.bNativeSurface||bUnknown)++Counts.UnknownBoxes;
+            else if(Verdict==EVoxelTileClass::AllSolid)++Counts.SolidBoxes;
+            else if(Verdict==EVoxelTileClass::AllAir)++Counts.AirBoxes;
+            else ++Counts.MixedBoxes;
+            return Result;
+        };
        // A box crossing recipes cannot use one field interval; unknown composition stays Mixed.
        // Une boite qui traverse plusieurs recettes reste Mixed, car l'intervalle d'un seul champ ne suffit pas.
-        if(bSurfaceOverlap&&bOtherStrateOverlap)return EVoxelTileClass::Mixed;
+        if(bBadPlan|| (bSurfaceOverlap&&bOtherStrateOverlap))
+            return Finish(EVoxelTileClass::Mixed,true);
        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;
+            return Finish(EVoxelTileClass::Mixed,true);
        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)
@@ -1175,32 +1216,39 @@ public:
                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;
+                if(Overlap)return Finish(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;
+                return Finish(EVoxelTileClass::Mixed);
        if(bSurfaceOverlap)
        {
-            if(!WorldAssembly||!SurfacePlan)return EVoxelTileClass::Mixed;
+            if(!WorldAssembly||!SurfacePlan)return Finish(EVoxelTileClass::Mixed,true);
            const VoxelForgeCore::Heightfield2D Block(WorldAssembly->recipe,*SurfacePlan);
-            const auto Bounds=Block.interval(
+            const auto Verdict=Block.classifyLatticeBox(
                static_cast<int64>(BMin[0])*VoxelForgeCore::WorldAssemblySubvoxelsPerVoxel,
                static_cast<int64>(BMin[1])*VoxelForgeCore::WorldAssemblySubvoxelsPerVoxel,
                static_cast<int64>(BMin[2])*VoxelForgeCore::WorldAssemblySubvoxelsPerVoxel,
                static_cast<int64>(BMax[0])*VoxelForgeCore::WorldAssemblySubvoxelsPerVoxel,
                static_cast<int64>(BMax[1])*VoxelForgeCore::WorldAssemblySubvoxelsPerVoxel,
-                static_cast<int64>(BMax[2])*VoxelForgeCore::WorldAssemblySubvoxelsPerVoxel);
-            if(!Bounds.known)return EVoxelTileClass::Mixed;
-            if(Bounds.minimumDensitySubvoxels>0)return EVoxelTileClass::AllAir;
-            if(Bounds.maximumDensitySubvoxels<0)return EVoxelTileClass::AllSolid;
-            return EVoxelTileClass::Mixed;
+                static_cast<int64>(BMax[2])*VoxelForgeCore::WorldAssemblySubvoxelsPerVoxel,
+                std::max(Step,1));
+            switch(Verdict)
+            {
+                case VoxelForgeCore::HeightfieldBoxClass::allAir:
+                    return Finish(EVoxelTileClass::AllAir);
+                case VoxelForgeCore::HeightfieldBoxClass::allSolid:
+                    return Finish(EVoxelTileClass::AllSolid);
+                case VoxelForgeCore::HeightfieldBoxClass::mixed:
+                    return Finish(EVoxelTileClass::Mixed);
+                default:return Finish(EVoxelTileClass::Mixed,true);
+            }
        }
        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;
+        return Finish(bRock?EVoxelTileClass::AllSolid:EVoxelTileClass::AllAir);
    }

    bool IsRock(int32 X,int32 Y,int32 Z) const
@@ -1662,6 +1710,12 @@ public:
    }

 private:
+    static FClassificationThreadStats& ThreadClassificationStats()
+    {
+        static thread_local FClassificationThreadStats Stats;
+        return Stats;
+    }
+
    const VoxelForgeCore::WorldAssemblyStratumPlanV1* FindWorldAssemblyStratum(
        int32 StrateIndex,int32 ChunkZ) const
    {
@@ -632,9 +632,59 @@ namespace VoxelForgeCore
        int64 sampleCount=0;std::vector<float> densityValues;
        const FVoxelTileKey nativeKey=ToNative(key);
        const FIntVector origin=nativeKey.OriginVoxels();
+        TileMesh out;out.key=key;
+        const bool bHasWorldAssembly=impl_->worldAssembly!=nullptr;
+        const bool bNeedsFarShell=bHasWorldAssembly&&key.lod>=7
+            &&(tileShape.HasFarWorld()||impl_->shape->HasFarWorld());
+        FStandaloneShapeWorld::ResetThreadClassificationStats();
+        VoxelForgeCore::Heightfield2D::resetThreadColumnCacheStats();
+        bool bNativeSurfaceTile=false;
+        if(bHasWorldAssembly)
+        {
+            const auto Classification=tileShape.ClassifyBoxDetailed(origin,step,cells);
+            bNativeSurfaceTile=Classification.bNativeSurface;
+            if(bNativeSurfaceTile)
+            {
+                out.diagnostics.tileWasClassified=true;
+                switch(Classification.bUnknown?EVoxelTileClass::Mixed:Classification.Verdict)
+                {
+                    case EVoxelTileClass::AllSolid:out.diagnostics.tileClass=TileClass::allSolid;break;
+                    case EVoxelTileClass::AllAir:out.diagnostics.tileClass=TileClass::allAir;break;
+                    case EVoxelTileClass::Mixed:out.diagnostics.tileClass=TileClass::mixed;break;
+                }
+                if((out.diagnostics.tileClass==TileClass::allSolid
+                    ||out.diagnostics.tileClass==TileClass::allAir)
+                    &&!options.hashDensity&&!options.captureDensitySamples&&!bNeedsFarShell)
+                {
+                    out.diagnostics.meshingSkipped=true;
+                    const auto Counts=FStandaloneShapeWorld::GetThreadClassificationStats();
+                    const auto Columns=VoxelForgeCore::Heightfield2D::threadColumnCacheStats();
+                    out.diagnostics.solidBoxes=Counts.SolidBoxes;
+                    out.diagnostics.airBoxes=Counts.AirBoxes;
+                    out.diagnostics.mixedBoxes=Counts.MixedBoxes;
+                    out.diagnostics.unknownBoxes=Counts.UnknownBoxes;
+                    out.diagnostics.fieldCalls=Columns.pointCalls;
+                    out.diagnostics.floatFieldCalls=Counts.NativeFloatFieldCalls;
+                    out.diagnostics.latticeFieldCalls=Counts.NativeLatticeFieldCalls;
+                    out.diagnostics.unclampedFieldCalls=Counts.NativeUnclampedFieldCalls;
+                    out.diagnostics.columnCacheHits=Columns.columnHits;
+                    out.diagnostics.columnCacheMisses=Columns.columnMisses;
+                    out.diagnostics.classificationMilliseconds=Counts.Milliseconds;
+                    out.diagnostics.columnComputeMilliseconds=Columns.columnComputeMilliseconds;
+                    out.densitySampleCount=0;
+                    return out;
+                }
+                // Native SurfaceWorld classifies each 8-cell block against its exact sample
+                // lattice. This lets the mesher avoid filling density grids for proven air/rock.
+                mesher.bUseBlockEarlyOut=true;
+            }
+        }
+        const auto MeshingStart=std::chrono::steady_clock::now();
        FVoxelMeshData mesh=mesher.GenerateMesh(origin,step,cells,nullptr,
            options.bandMinZVoxel,options.bandMaxZVoxel,&sampleCount,nullptr,nullptr,
            (options.hashDensity||options.captureDensitySamples)?&densityValues:nullptr);
+        const double MeshingMilliseconds=std::chrono::duration<double,std::milli>(
+            std::chrono::steady_clock::now()-MeshingStart).count();
        const FStandaloneShapeWorld& farShape=tileShape.HasFarWorld()?tileShape:*impl_->shape;
        if(farShape.HasFarWorld())
        {
@@ -642,7 +692,33 @@ namespace VoxelForgeCore
            AppendFarConnectorShells(farShape,nativeKey,step,world.CoarseTileCells,mesh);
        }

-        TileMesh out;out.key=key;out.densitySampleCount=static_cast<std::uint64_t>(std::max<int64>(0,sampleCount));
+        out.densitySampleCount=static_cast<std::uint64_t>(std::max<int64>(0,sampleCount));
+        if(bHasWorldAssembly)
+        {
+            const auto Counts=FStandaloneShapeWorld::GetThreadClassificationStats();
+            const auto Columns=VoxelForgeCore::Heightfield2D::threadColumnCacheStats();
+            out.diagnostics.fieldCalls=Columns.pointCalls;
+            out.diagnostics.floatFieldCalls=Counts.NativeFloatFieldCalls;
+            out.diagnostics.latticeFieldCalls=Counts.NativeLatticeFieldCalls;
+            out.diagnostics.unclampedFieldCalls=Counts.NativeUnclampedFieldCalls;
+            out.diagnostics.columnCacheHits=Columns.columnHits;
+            out.diagnostics.columnCacheMisses=Columns.columnMisses;
+            out.diagnostics.columnComputeMilliseconds=Columns.columnComputeMilliseconds;
+            out.diagnostics.classificationMilliseconds=Counts.Milliseconds;
+            if(bNativeSurfaceTile)
+            {
+                out.diagnostics.solidBoxes=Counts.SolidBoxes;
+                out.diagnostics.airBoxes=Counts.AirBoxes;
+                out.diagnostics.mixedBoxes=Counts.MixedBoxes;
+                out.diagnostics.unknownBoxes=Counts.UnknownBoxes;
+            }
+            else
+            {
+                out.diagnostics.unknownBoxes=Counts.UnknownBoxes;
+                out.diagnostics.tileClass=TileClass::unknown;
+            }
+            out.diagnostics.meshingMilliseconds=MeshingMilliseconds;
+        }
        out.meshPayloadAllocatedBytes=static_cast<std::uint64_t>(mesh.Vertices.GetAllocatedSize())
            +static_cast<std::uint64_t>(mesh.Normals.GetAllocatedSize())
            +static_cast<std::uint64_t>(mesh.UVs.GetAllocatedSize())
@@ -4,7 +4,9 @@
 #include "VoxelForgeCore/VoxelForgeCore.h"

 #include <algorithm>
+#include <bit>
 #include <charconv>
+#include <chrono>
 #include <cmath>
 #include <limits>
 #include <sstream>
@@ -99,7 +101,8 @@ namespace world_assembly_detail
    inline bool ValidateRecipe(const WorldAssemblyRecipeV1& recipe, std::string& error)
    {
        if (recipe.schemaVersion != WorldAssemblySchemaVersion
-            || (recipe.evaluatorVersion != 1 && recipe.evaluatorVersion != WorldAssemblyEvaluatorVersion)
+            || recipe.evaluatorVersion == 0
+            || recipe.evaluatorVersion > WorldAssemblyEvaluatorVersion
            || recipe.unitsPerVoxel != WorldAssemblySubvoxelsPerVoxel)
        { error = "world assembly recipe version, evaluator or fixed scale is unsupported"; return false; }
        if (recipe.biomes.empty() || recipe.biomes.size() > MaxBiomeCount)
@@ -183,7 +186,7 @@ namespace world_assembly_detail
        ids.reserve(stratum.features.size());
        for (const auto& feature : stratum.features)
        {
-            if (feature.stableId == 0 || feature.type < 1 || feature.type > 3
+            if (feature.stableId == 0 || feature.type < 1 || feature.type > 4
                || std::find(ids.begin(), ids.end(), feature.stableId) != ids.end())
            { error = "world assembly feature IDs and types must be stable, unique and supported"; return false; }
            ids.push_back(feature.stableId);
@@ -194,7 +197,7 @@ namespace world_assembly_detail
                { error = "world assembly feature has an invalid influence bound"; return false; }
            if (!std::all_of(feature.parametersSubvoxels.begin(), feature.parametersSubvoxels.end(), ValidScalar))
            { error = "world assembly feature parameter is outside the fixed-scale range"; return false; }
-            if (feature.type == 2)
+            if (feature.type == 2 || feature.type == 4)
            {
                const std::int64_t stratumMinZ = static_cast<std::int64_t>(stratum.bottomChunkZ)
                    * 32 * WorldAssemblySubvoxelsPerVoxel;
@@ -210,6 +213,12 @@ namespace world_assembly_detail
                if (feature.minimumSubvoxels[2] < stratumMinZ
                    || feature.maximumSubvoxels[2] > stratumMaxZ)
                { error = "world assembly overhang exception escapes its finite stratum Z bounds"; return false; }
+                if (feature.type == 4
+                    && (std::abs(feature.parametersSubvoxels[0])
+                            + std::abs(feature.parametersSubvoxels[1]) != 256
+                        || feature.parametersSubvoxels[2] <= 0
+                        || feature.parametersSubvoxels[3] <= 0))
+                { error = "world assembly cliff shelf requires a unit horizontal axis and positive height/reach"; return false; }
            }
        }
        for (const auto& point : stratum.route)
@@ -954,6 +963,200 @@ namespace world_assembly_detail
        return SmoothQ16(static_cast<std::uint32_t>(std::clamp<std::int64_t>(edge*Q16/feather,0,Q16)));
    }

+    struct HeightfieldColumnCacheEntry
+    {
+        const WorldAssemblyRecipeV1* recipe = nullptr;
+        const WorldAssemblyStratumPlanV1* plan = nullptr;
+        std::uint64_t recipeHash = 0, planHash = 0;
+        std::int64_t x = 0, y = 0;
+        HeightfieldPointSample column{};
+    };
+    // Coarse LOD tiles revisit a wider set of gradient-sample columns than the full-resolution
+    // band. Keep those columns resident for a complete worker tile so the 3D mesher does not
+    // repeat the procedural biome evaluation for every Z sample.
+    inline thread_local std::array<HeightfieldColumnCacheEntry, 65536> HeightfieldColumnCache{};
+    inline thread_local HeightfieldColumnCacheStats HeightfieldColumnStats{};
+
+    inline HeightfieldPointSample EvaluateColumn(const WorldAssemblyRecipeV1& recipe,
+        const WorldAssemblyStratumPlanV1& plan,std::int64_t x,std::int64_t y)
+    {
+        HeightfieldPointSample out;
+        if(recipe.evaluatorVersion>=2)
+        {
+            const BiomeSelectionV2 selection=SelectBiomeV2(recipe,x,y);
+            const auto& primary=recipe.biomes[selection.primary];
+            const auto& neighbor=recipe.biomes[selection.neighbor];
+            const std::int64_t primaryHeight=SurfaceHeightV2(recipe,plan,primary,x,y);
+            const std::int64_t neighborHeight=selection.primary==selection.neighbor?primaryHeight
+                :SurfaceHeightV2(recipe,plan,neighbor,x,y);
+            std::int64_t ground=primaryHeight+((neighborHeight-primaryHeight)
+                *(static_cast<std::int64_t>(Q16)-selection.primaryWeightQ16))/Q16;
+            const std::int64_t primaryCap=CapHeightV2(recipe,plan,primary,x,y);
+            const std::int64_t neighborCap=selection.primary==selection.neighbor?primaryCap
+                :CapHeightV2(recipe,plan,neighbor,x,y);
+            std::int64_t cap=primaryCap+((neighborCap-primaryCap)
+                *(static_cast<std::int64_t>(Q16)-selection.primaryWeightQ16))/Q16;
+            for(const auto& feature:plan.features)
+            {
+                if(feature.type==1)
+                {
+                    const std::uint32_t fade=FeatureFadeQ16(feature,x,y,feature.parametersSubvoxels[1]);
+                    if(fade)ground+=MulQ16(feature.parametersSubvoxels[0],fade);
+                }
+                else if(feature.type==3)
+                {
+                    const std::uint32_t fade=FeatureFadeQ16(feature,x,y,feature.parametersSubvoxels[2]);
+                    const std::int64_t wave=std::max<std::int64_t>(1,feature.parametersSubvoxels[1]);
+                    if(fade)
+                    {
+                        const std::uint32_t featureSeed=static_cast<std::uint32_t>(feature.stableId);
+                        const std::int64_t nx=NoiseCoordinateQ16(x,wave)+SeedOffsetQ16(featureSeed,1700);
+                        const std::int64_t ny=NoiseCoordinateQ16(y,wave)+SeedOffsetQ16(featureSeed,2900);
+                        const std::int32_t noise=FractalNoiseQ16(nx,ny,
+                            SeedOffsetQ16(featureSeed,4100),3);
+                        cap+=MulQ16(MulQ16(noise,feature.parametersSubvoxels[0]),fade);
+                    }
+                }
+            }
+            out.groundHeightSubvoxels=ground;
+            out.capHeightSubvoxels=cap;
+            out.biomeId=primary.stableId;
+            out.neighborBiomeId=neighbor.stableId;
+            out.biomeWeightQ16=selection.primaryWeightQ16;
+            out.waterSurfaceSubvoxels=plan.waterSurfaceSubvoxels;
+            return out;
+        }
+
+        const std::int64_t warpX=MulQ16(recipe.biomeWarpSubvoxels,
+            Noise2D(x,y,recipe.biomeWarpWavelengthSubvoxels,recipe.seed,0x574152505f585f31ull));
+        const std::int64_t warpY=MulQ16(recipe.biomeWarpSubvoxels,
+            Noise2D(x,y,recipe.biomeWarpWavelengthSubvoxels,recipe.seed,0x574152505f595f31ull));
+        const std::int64_t warpedX=x+warpX,warpedY=y+warpY;
+        const std::int64_t cellX=FloorDiv(warpedX,recipe.biomeCellSizeSubvoxels);
+        const std::int64_t cellY=FloorDiv(warpedY,recipe.biomeCellSizeSubvoxels);
+        std::uint64_t nearestDistance=std::numeric_limits<std::uint64_t>::max();
+        std::uint64_t secondDistance=std::numeric_limits<std::uint64_t>::max();
+        std::uint32_t nearestId=0,secondId=0;
+        const std::uint32_t jitterSeed=recipe.seed^0xa341316cu;
+        const std::int64_t cellSize=recipe.biomeCellSizeSubvoxels;
+        for(std::int64_t cy=cellY-2;cy<=cellY+2;++cy)
+        for(std::int64_t cx=cellX-2;cx<=cellX+2;++cx)
+        {
+            const std::uint64_t h=CoordinateHash(jitterSeed,cx,cy,0x42494f4d455f5349ull);
+            const std::int64_t jitterX=static_cast<std::int64_t>(h&0xffffu)*(cellSize/2)/65535-cellSize/4;
+            const std::int64_t jitterY=static_cast<std::int64_t>((h>>16)&0xffffu)*(cellSize/2)/65535-cellSize/4;
+            const std::int64_t siteX=cx*cellSize+cellSize/2+jitterX;
+            const std::int64_t siteY=cy*cellSize+cellSize/2+jitterY;
+            const std::int64_t dx=warpedX-siteX,dy=warpedY-siteY;
+            const std::uint64_t distance=static_cast<std::uint64_t>(dx*dx+dy*dy);
+            const std::uint32_t biomeId=recipe.biomes[CoordinateHash(jitterSeed,cx,cy,
+                0x42494f4d455f4944ull)%recipe.biomes.size()].stableId;
+            if(distance<nearestDistance)
+            {secondDistance=nearestDistance;secondId=nearestId;nearestDistance=distance;nearestId=biomeId;}
+            else if(distance<secondDistance&&biomeId!=nearestId)
+            {secondDistance=distance;secondId=biomeId;}
+        }
+        const auto findBiome=[&](std::uint32_t id)->const WorldAssemblyBiomeProfileV1*
+        {
+            const auto found=std::find_if(recipe.biomes.begin(),recipe.biomes.end(),
+                [id](const auto& biome){return biome.stableId==id;});
+            return found==recipe.biomes.end()?nullptr:&*found;
+        };
+        const auto* primary=findBiome(nearestId);
+        const auto* secondary=findBiome(secondId?secondId:nearestId);
+        if(!primary||!secondary)return out;
+        std::uint32_t primaryWeight=65536;
+        if(recipe.biomeBorderBlendSubvoxels>0&&nearestDistance<secondDistance)
+        {
+            const std::uint64_t difference=secondDistance-nearestDistance;
+            const std::uint64_t denominator=static_cast<std::uint64_t>(2*cellSize)
+                *static_cast<std::uint64_t>(recipe.biomeBorderBlendSubvoxels);
+            primaryWeight=SmoothQ16(RatioQ16(difference,denominator));
+        }
+        const std::int64_t primaryHeight=ProfileHeight(*primary,warpedX,warpedY,recipe.seed);
+        const std::int64_t secondaryHeight=ProfileHeight(*secondary,warpedX,warpedY,recipe.seed);
+        out.groundHeightSubvoxels=(primaryHeight*primaryWeight
+            +secondaryHeight*(65536-primaryWeight))/65536;
+        out.biomeId=primary->stableId;out.neighborBiomeId=secondary->stableId;
+        out.biomeWeightQ16=primaryWeight;
+        out.capHeightSubvoxels=plan.capHeightSubvoxels+MulQ16(plan.capUndulationSubvoxels,
+            Noise2D(x,y,plan.capWavelengthSubvoxels,recipe.seed,0x4341505f554e4431ull));
+        out.waterSurfaceSubvoxels=plan.waterSurfaceSubvoxels;
+        return out;
+    }
+
+    inline HeightfieldPointSample CachedColumn(const WorldAssemblyRecipeV1& recipe,
+        const WorldAssemblyStratumPlanV1& plan,std::int64_t x,std::int64_t y)
+    {
+        const std::uint64_t mixed=(static_cast<std::uint64_t>(x)*0x9E3779B185EBCA87ull)
+            ^std::rotl(static_cast<std::uint64_t>(y)*0xC2B2AE3D27D4EB4Full,23)
+            ^static_cast<std::uint64_t>(reinterpret_cast<std::uintptr_t>(&recipe))
+            ^std::rotl(static_cast<std::uint64_t>(reinterpret_cast<std::uintptr_t>(&plan)),41);
+        auto& entry=HeightfieldColumnCache[static_cast<std::size_t>(mixed)&
+            (HeightfieldColumnCache.size()-1)];
+        if(entry.recipe==&recipe&&entry.plan==&plan&&entry.recipeHash==recipe.contentHash
+            &&entry.planHash==plan.contentHash&&entry.x==x&&entry.y==y)
+        {++HeightfieldColumnStats.columnHits;return entry.column;}
+        const auto start=std::chrono::steady_clock::now();
+        HeightfieldPointSample column=EvaluateColumn(recipe,plan,x,y);
+        HeightfieldColumnStats.columnComputeMilliseconds+=std::chrono::duration<double,std::milli>(
+            std::chrono::steady_clock::now()-start).count();
+        ++HeightfieldColumnStats.columnMisses;
+        entry={&recipe,&plan,recipe.contentHash,plan.contentHash,x,y,column};
+        return column;
+    }
+
+    inline std::int64_t CliffShelfDensity(const WorldAssemblyFeatureV1& feature,
+        const WorldAssemblyRecipeV1& recipe,const WorldAssemblyStratumPlanV1& plan,
+        std::int64_t x,std::int64_t y,std::int64_t z,
+        const HeightfieldPointSample& localColumn)
+    {
+        const auto& lo=feature.minimumSubvoxels;const auto& hi=feature.maximumSubvoxels;
+        if(x<lo[0]||x>hi[0]||y<lo[1]||y>hi[1]||z<lo[2]||z>hi[2])
+            return std::numeric_limits<std::int64_t>::max()/4;
+        constexpr std::int64_t Feather=8*WorldAssemblySubvoxelsPerVoxel;
+        const std::int64_t dirX=feature.parametersSubvoxels[0]/WorldAssemblySubvoxelsPerVoxel;
+        const std::int64_t dirY=feature.parametersSubvoxels[1]/WorldAssemblySubvoxelsPerVoxel;
+        const std::int64_t height=feature.parametersSubvoxels[2];
+        const std::int64_t reach=feature.parametersSubvoxels[3];
+        if((dirX==0&&dirY==0)||height<=0||reach<=0)return std::numeric_limits<std::int64_t>::max()/4;
+        const std::int64_t centerX=(lo[0]+hi[0])/2,centerY=(lo[1]+hi[1])/2;
+        const std::int64_t edgeX=dirX>0?lo[0]+Feather:dirX<0?hi[0]-Feather:centerX;
+        const std::int64_t edgeY=dirY>0?lo[1]+Feather:dirY<0?hi[1]-Feather:centerY;
+        const std::int64_t dx=x-edgeX,dy=y-edgeY;
+        const std::int64_t along=dx*dirX+dy*dirY;
+        const std::int64_t lateral=-dx*dirY+dy*dirX;
+        const std::int64_t halfWidth=std::max<std::int64_t>(0,
+            (dirX!=0?(hi[1]-lo[1]):(hi[0]-lo[0]))/2-Feather);
+        const std::int64_t aboveGround=z-localColumn.groundHeightSubvoxels;
+        if(aboveGround<=0||aboveGround>=height||along<0||along>reach
+            ||std::abs(lateral)>halfWidth)
+            return std::numeric_limits<std::int64_t>::max()/4;
+        const std::uint32_t vertical=SmoothQ16(static_cast<std::uint32_t>(
+            std::clamp<std::int64_t>(aboveGround*Q16/height,0,Q16)));
+        const std::int64_t sourceOffset=MulQ16(reach,vertical);
+        const std::int64_t sourceX=x-dirX*sourceOffset;
+        const std::int64_t sourceY=y-dirY*sourceOffset;
+        const HeightfieldPointSample source=CachedColumn(recipe,plan,sourceX,sourceY);
+        if(source.groundHeightSubvoxels<=localColumn.groundHeightSubvoxels)
+            return std::numeric_limits<std::int64_t>::max()/4;
+        std::int64_t candidate=std::max({localColumn.groundHeightSubvoxels-z,
+            z-source.groundHeightSubvoxels,-along,along-reach,
+            std::abs(lateral)-halfWidth,-aboveGround,aboveGround-height});
+        const std::uint32_t edgeFade=FeatureFadeQ16(feature,x,y,Feather);
+        candidate+=static_cast<std::int64_t>(Q16-edgeFade)*Feather/Q16;
+        return candidate;
+    }
+
+    inline std::int64_t SmoothMinSubvoxels(std::int64_t a,std::int64_t b,
+        std::int64_t width)
+    {
+        if(width<=0)return std::min(a,b);
+        const std::int64_t h=std::max<std::int64_t>(width-std::abs(a-b),0);
+        const std::int64_t smoothQ16=h*Q16/width;
+        return std::min(a,b)-MulQ16(MulQ16(smoothQ16,smoothQ16),width/4);
+    }
+
    inline bool InInterval(const FieldInterval& interval, std::int64_t density)
    { return interval.known && density >= interval.minimumDensitySubvoxels
        && density <= interval.maximumDensitySubvoxels; }
@@ -1143,7 +1346,8 @@ Heightfield2D::Heightfield2D(const WorldAssemblyRecipeV1& recipe,
 bool Heightfield2D::isValid() const noexcept
 {
    if (!recipe_ || !plan_ || recipe_->schemaVersion != WorldAssemblySchemaVersion
-        || (recipe_->evaluatorVersion != 1 && recipe_->evaluatorVersion != WorldAssemblyEvaluatorVersion)
+        || recipe_->evaluatorVersion == 0
+        || recipe_->evaluatorVersion > WorldAssemblyEvaluatorVersion
        || recipe_->unitsPerVoxel != WorldAssemblySubvoxelsPerVoxel
        || recipe_->biomes.empty() || recipe_->biomes.size() > world_assembly_detail::MaxBiomeCount
        || recipe_->biomeCellSizeSubvoxels <= 0 || recipe_->biomeWarpWavelengthSubvoxels <= 0
@@ -1169,125 +1373,88 @@ HeightfieldPointSample Heightfield2D::point(std::int64_t x, std::int64_t y,
        || world_assembly_detail::OutsideRange(y)
        || world_assembly_detail::OutsideRange(z)) return out;
    using namespace world_assembly_detail;
+    ++HeightfieldColumnStats.pointCalls;
+    out=CachedColumn(*recipe_,*plan_,x,y);
+    out.densitySubvoxels=std::min(z-out.groundHeightSubvoxels,
+        out.capHeightSubvoxels-z);
    if(recipe_->evaluatorVersion>=2)
    {
-        const BiomeSelectionV2 selection=SelectBiomeV2(*recipe_,x,y);
-        const auto& primary=recipe_->biomes[selection.primary];
-        const auto& neighbor=recipe_->biomes[selection.neighbor];
-        const std::int64_t primaryHeight=SurfaceHeightV2(*recipe_,*plan_,primary,x,y);
-        const std::int64_t neighborHeight=selection.primary==selection.neighbor?primaryHeight
-            :SurfaceHeightV2(*recipe_,*plan_,neighbor,x,y);
-        std::int64_t ground=primaryHeight+((neighborHeight-primaryHeight)
-            *(static_cast<std::int64_t>(Q16)-selection.primaryWeightQ16))/Q16;
-        const std::int64_t primaryCap=CapHeightV2(*recipe_,*plan_,primary,x,y);
-        const std::int64_t neighborCap=selection.primary==selection.neighbor?primaryCap
-            :CapHeightV2(*recipe_,*plan_,neighbor,x,y);
-        std::int64_t cap=primaryCap+((neighborCap-primaryCap)
-            *(static_cast<std::int64_t>(Q16)-selection.primaryWeightQ16))/Q16;
        for(const auto& feature:plan_->features)
        {
-            if(feature.type==1)
+            if(feature.type==2)
            {
-                const std::uint32_t fade=FeatureFadeQ16(feature,x,y,feature.parametersSubvoxels[1]);
-                if(fade)ground+=MulQ16(feature.parametersSubvoxels[0],fade);
+                const auto& lo=feature.minimumSubvoxels;const auto& hi=feature.maximumSubvoxels;
+                if(x<lo[0]||x>hi[0]||y<lo[1]||y>hi[1]||z<lo[2]||z>hi[2])continue;
+                const std::int64_t depth=std::min({x-lo[0],hi[0]-x,y-lo[1],hi[1]-y,
+                    z-lo[2],hi[2]-z});
+                out.densitySubvoxels=std::min(out.densitySubvoxels,-depth);
            }
-            else if(feature.type==3)
+            else if(feature.type==4&&recipe_->evaluatorVersion>=3)
            {
-                const std::uint32_t fade=FeatureFadeQ16(feature,x,y,feature.parametersSubvoxels[2]);
-                const std::int64_t wave=std::max<std::int64_t>(1,feature.parametersSubvoxels[1]);
-                if(fade)
-                {
-                    const std::uint32_t featureSeed=static_cast<std::uint32_t>(feature.stableId);
-                    const std::int64_t nx=NoiseCoordinateQ16(x,wave)+SeedOffsetQ16(featureSeed,1700);
-                    const std::int64_t ny=NoiseCoordinateQ16(y,wave)+SeedOffsetQ16(featureSeed,2900);
-                    const std::int32_t noise=FractalNoiseQ16(nx,ny,SeedOffsetQ16(featureSeed,4100),3);
-                    cap+=MulQ16(MulQ16(noise,feature.parametersSubvoxels[0]),fade);
-                }
+                const std::int64_t shelf=CliffShelfDensity(feature,*recipe_,*plan_,x,y,z,out);
+                out.densitySubvoxels=SmoothMinSubvoxels(out.densitySubvoxels,shelf,
+                    4*WorldAssemblySubvoxelsPerVoxel);
            }
        }
-        out.groundHeightSubvoxels=ground;
-        out.capHeightSubvoxels=cap;
-        out.biomeId=primary.stableId;
-        out.neighborBiomeId=neighbor.stableId;
-        out.biomeWeightQ16=selection.primaryWeightQ16;
-        out.densitySubvoxels=std::min(z-ground,cap-z);
-        for(const auto& feature:plan_->features)
-        {
-            if(feature.type!=2)continue;
-            const auto& lo=feature.minimumSubvoxels;const auto& hi=feature.maximumSubvoxels;
-            if(x<lo[0]||x>hi[0]||y<lo[1]||y>hi[1]||z<lo[2]||z>hi[2])continue;
-            const std::int64_t depth=std::min({x-lo[0],hi[0]-x,y-lo[1],hi[1]-y,z-lo[2],hi[2]-z});
-            out.densitySubvoxels=std::min(out.densitySubvoxels,-depth);
-        }
-        out.waterSurfaceSubvoxels=plan_->waterSurfaceSubvoxels;
-        if(plan_->waterEnabled && out.densitySubvoxels>=0 && z>ground
-            &&z<plan_->waterSurfaceSubvoxels && ground<plan_->waterSurfaceSubvoxels)
-            out.waterFillSubvoxels=std::min(z-ground,plan_->waterSurfaceSubvoxels-z);
-        return out;
    }
-    const std::int64_t warpX = MulQ16(recipe_->biomeWarpSubvoxels,
-        Noise2D(x, y, recipe_->biomeWarpWavelengthSubvoxels, recipe_->seed, 0x574152505f585f31ull));
-    const std::int64_t warpY = MulQ16(recipe_->biomeWarpSubvoxels,
-        Noise2D(x, y, recipe_->biomeWarpWavelengthSubvoxels, recipe_->seed, 0x574152505f595f31ull));
-    const std::int64_t warpedX = x + warpX, warpedY = y + warpY;
-    const std::int64_t cellX = FloorDiv(warpedX, recipe_->biomeCellSizeSubvoxels);
-    const std::int64_t cellY = FloorDiv(warpedY, recipe_->biomeCellSizeSubvoxels);
-    std::uint64_t nearestDistance = std::numeric_limits<std::uint64_t>::max();
-    std::uint64_t secondDistance = std::numeric_limits<std::uint64_t>::max();
-    std::uint32_t nearestId = 0, secondId = 0;
-    const std::uint32_t jitterSeed = recipe_->seed ^ 0xa341316cu;
-    const std::int64_t cellSize = recipe_->biomeCellSizeSubvoxels;
-    for (std::int64_t cy = cellY - 2; cy <= cellY + 2; ++cy)
-    for (std::int64_t cx = cellX - 2; cx <= cellX + 2; ++cx)
-    {
-        const std::uint64_t h = CoordinateHash(jitterSeed, cx, cy, 0x42494f4d455f5349ull);
-        const std::int64_t jitterX = static_cast<std::int64_t>((h & 0xffffu)) * (cellSize / 2) / 65535 - cellSize / 4;
-        const std::int64_t jitterY = static_cast<std::int64_t>(((h >> 16) & 0xffffu)) * (cellSize / 2) / 65535 - cellSize / 4;
-        const std::int64_t siteX = cx * cellSize + cellSize / 2 + jitterX;
-        const std::int64_t siteY = cy * cellSize + cellSize / 2 + jitterY;
-        const std::int64_t dx = warpedX - siteX, dy = warpedY - siteY;
-        const std::uint64_t distance = static_cast<std::uint64_t>(dx * dx + dy * dy);
-        const std::uint32_t biomeId = recipe_->biomes[CoordinateHash(jitterSeed, cx, cy,
-            0x42494f4d455f4944ull) % recipe_->biomes.size()].stableId;
-        if (distance < nearestDistance)
-        {
-            secondDistance = nearestDistance; secondId = nearestId;
-            nearestDistance = distance; nearestId = biomeId;
-        }
-        else if (distance < secondDistance && biomeId != nearestId)
-        { secondDistance = distance; secondId = biomeId; }
-    }
-    const auto findBiome = [&](std::uint32_t id) -> const WorldAssemblyBiomeProfileV1*
-    {
-        const auto found = std::find_if(recipe_->biomes.begin(), recipe_->biomes.end(),
-            [id](const auto& biome) { return biome.stableId == id; });
-        return found == recipe_->biomes.end() ? nullptr : &*found;
-    };
-    const auto* primary = findBiome(nearestId);
-    const auto* secondary = findBiome(secondId ? secondId : nearestId);
-    if (!primary || !secondary) return out;
-    std::uint32_t primaryWeight = 65536;
-    if (recipe_->biomeBorderBlendSubvoxels > 0 && nearestDistance < secondDistance)
-    {
-        const std::uint64_t difference = secondDistance - nearestDistance;
-        const std::uint64_t denominator = static_cast<std::uint64_t>(2 * cellSize)
-            * static_cast<std::uint64_t>(recipe_->biomeBorderBlendSubvoxels);
-        primaryWeight = SmoothQ16(RatioQ16(difference, denominator));
-    }
-    const std::int64_t primaryHeight = ProfileHeight(*primary, warpedX, warpedY, recipe_->seed);
-    const std::int64_t secondaryHeight = ProfileHeight(*secondary, warpedX, warpedY, recipe_->seed);
-    out.groundHeightSubvoxels = (primaryHeight * primaryWeight
-        + secondaryHeight * (65536 - primaryWeight)) / 65536;
-    out.biomeId = primary->stableId;
-    out.neighborBiomeId = secondary->stableId;
-    out.biomeWeightQ16 = primaryWeight;
-    out.capHeightSubvoxels = plan_->capHeightSubvoxels + MulQ16(plan_->capUndulationSubvoxels,
-        Noise2D(x, y, plan_->capWavelengthSubvoxels, recipe_->seed, 0x4341505f554e4431ull));
-    out.densitySubvoxels = std::min(z - out.groundHeightSubvoxels,
-        out.capHeightSubvoxels - z);
+    if(recipe_->evaluatorVersion>=2&&plan_->waterEnabled&&out.densitySubvoxels>=0
+        &&z>out.groundHeightSubvoxels&&z<plan_->waterSurfaceSubvoxels
+        &&out.groundHeightSubvoxels<plan_->waterSurfaceSubvoxels)
+        out.waterFillSubvoxels=std::min(z-out.groundHeightSubvoxels,
+            plan_->waterSurfaceSubvoxels-z);
    return out;
 }

+void Heightfield2D::resetThreadColumnCacheStats() noexcept
+{ world_assembly_detail::HeightfieldColumnStats={}; }
+
+HeightfieldColumnCacheStats Heightfield2D::threadColumnCacheStats() noexcept
+{ return world_assembly_detail::HeightfieldColumnStats; }
+
+HeightfieldBoxClass Heightfield2D::classifyLatticeBox(std::int64_t minimumX,
+    std::int64_t minimumY,std::int64_t minimumZ,std::int64_t maximumX,
+    std::int64_t maximumY,std::int64_t maximumZ,std::int32_t stepVoxels) const noexcept
+{
+    using namespace world_assembly_detail;
+    if(!isValid()||stepVoxels<=0||minimumX>maximumX||minimumY>maximumY||minimumZ>maximumZ
+        ||OutsideRange(minimumX)||OutsideRange(minimumY)||OutsideRange(minimumZ)
+        ||OutsideRange(maximumX)||OutsideRange(maximumY)||OutsideRange(maximumZ))
+        return HeightfieldBoxClass::unknown;
+    for(const auto& feature:plan_->features)
+    {
+        if(feature.type!=2&&feature.type!=4)continue;
+        if(feature.minimumSubvoxels[0]<=maximumX&&feature.maximumSubvoxels[0]>=minimumX
+            &&feature.minimumSubvoxels[1]<=maximumY&&feature.maximumSubvoxels[1]>=minimumY
+            &&feature.minimumSubvoxels[2]<=maximumZ&&feature.maximumSubvoxels[2]>=minimumZ)
+            return HeightfieldBoxClass::mixed;
+    }
+    const std::int64_t stride=static_cast<std::int64_t>(stepVoxels)
+        *WorldAssemblySubvoxelsPerVoxel;
+    if((maximumX-minimumX)%stride!=0||(maximumY-minimumY)%stride!=0
+        ||(maximumX-minimumX)/stride>2048||(maximumY-minimumY)/stride>2048)
+        return HeightfieldBoxClass::unknown;
+    std::int64_t groundMin=std::numeric_limits<std::int64_t>::max();
+    std::int64_t groundMax=std::numeric_limits<std::int64_t>::min();
+    std::int64_t capMin=std::numeric_limits<std::int64_t>::max();
+    std::int64_t capMax=std::numeric_limits<std::int64_t>::min();
+    for(std::int64_t y=minimumY;;y+=stride)
+    {
+        for(std::int64_t x=minimumX;;x+=stride)
+        {
+            const auto column=point(x,y,0);
+            groundMin=std::min(groundMin,column.groundHeightSubvoxels);
+            groundMax=std::max(groundMax,column.groundHeightSubvoxels);
+            capMin=std::min(capMin,column.capHeightSubvoxels);
+            capMax=std::max(capMax,column.capHeightSubvoxels);
+            if(x==maximumX)break;
+        }
+        if(y==maximumY)break;
+    }
+    if(minimumZ>groundMax&&maximumZ<capMin)return HeightfieldBoxClass::allAir;
+    if(maximumZ<groundMin||minimumZ>capMax)return HeightfieldBoxClass::allSolid;
+    return HeightfieldBoxClass::mixed;
+}
+
 FieldInterval Heightfield2D::interval(std::int64_t minX, std::int64_t minY,
    std::int64_t minZ, std::int64_t maxX, std::int64_t maxY, std::int64_t maxZ) const noexcept
 {
@@ -1361,16 +1528,23 @@ FieldInterval Heightfield2D::interval(std::int64_t minX, std::int64_t minY,
    {
        for(const auto& feature:plan_->features)
        {
-            if(feature.type!=2)continue;
+            if(feature.type!=2&&feature.type!=4)continue;
            bool overlap=true;
            for(std::size_t axis=0;axis<3;++axis)
                overlap=overlap&&feature.minimumSubvoxels[axis]<=std::array<std::int64_t,3>{maxX,maxY,maxZ}[axis]
                    &&feature.maximumSubvoxels[axis]>=std::array<std::int64_t,3>{minX,minY,minZ}[axis];
            if(!overlap)continue;
-            const std::int64_t depth=std::min({(feature.maximumSubvoxels[0]-feature.minimumSubvoxels[0])/2,
-                (feature.maximumSubvoxels[1]-feature.minimumSubvoxels[1])/2,
-                (feature.maximumSubvoxels[2]-feature.minimumSubvoxels[2])/2});
-            out.minimumDensitySubvoxels=std::min(out.minimumDensitySubvoxels,-depth);
+            if(feature.type==2)
+            {
+                const std::int64_t depth=std::min({
+                    (feature.maximumSubvoxels[0]-feature.minimumSubvoxels[0])/2,
+                    (feature.maximumSubvoxels[1]-feature.minimumSubvoxels[1])/2,
+                    (feature.maximumSubvoxels[2]-feature.minimumSubvoxels[2])/2});
+                out.minimumDensitySubvoxels=std::min(out.minimumDensitySubvoxels,-depth);
+            }
+            else
+                out.minimumDensitySubvoxels=std::min(out.minimumDensitySubvoxels,
+                    -feature.parametersSubvoxels[2]-WorldAssemblySubvoxelsPerVoxel);
        }
    }
    return out;
@@ -1416,6 +1590,8 @@ FieldDetailEnvelope Heightfield2D::detailEnvelope() const noexcept
                (feature.maximumSubvoxels[0]-feature.minimumSubvoxels[0])/2,
                (feature.maximumSubvoxels[1]-feature.minimumSubvoxels[1])/2,
                (feature.maximumSubvoxels[2]-feature.minimumSubvoxels[2])/2}));
+        if(feature.type==4)maximum=std::max(maximum,feature.parametersSubvoxels[2]
+            +WorldAssemblySubvoxelsPerVoxel);
        if(feature.type==3)maximum=std::max(maximum,std::abs(feature.parametersSubvoxels[0]));
    }
    out.maximumAbsoluteDisplacementSubvoxels = maximum;
@@ -1540,7 +1716,7 @@ WorldAssemblyBoundsAuditResult AuditHeightfieldBounds(const Heightfield2D& block
    }
    for(const auto& feature:plan->features)
    {
-        if(feature.type!=2)continue;
+        if(feature.type!=2&&feature.type!=4)continue;
        const std::int64_t x=feature.minimumSubvoxels[0]
            +(feature.maximumSubvoxels[0]-feature.minimumSubvoxels[0])/2;
        const std::int64_t y=feature.minimumSubvoxels[1]
@@ -1563,7 +1739,7 @@ WorldAssemblyBoundsAuditResult AuditHeightfieldBounds(const Heightfield2D& block
        ||result.biomeBorderBoxes>0;
    const bool cliffCoverage=!hasCliffProfile||result.cliffBoxes>0;
    const bool exceptionCoverage=std::none_of(plan->features.begin(),plan->features.end(),
-        [](const auto& feature){return feature.type==2;})||result.exceptionBoxes>0;
+        [](const auto& feature){return feature.type==2||feature.type==4;})||result.exceptionBoxes>0;
    result.passed = result.boxesTested >= boxCount && result.negativeCoordinateBoxes>0
        && borderCoverage&&cliffCoverage&&result.capBoxes>0&&exceptionCoverage
        && result.missedBounds == 0