feat(core): ore vein records and a materialAt query, off by default (L50); research R18, R19

Versioned capsule vein records (stable ID, resource ID, endpoints, radius, seeded irregularity) resolved
into the snapshot only when requested, placed beside the route and away from passages, landings and
landmarks (forced-passage control rejected). The bound field table gains materialAt (host + two resource
IDs + Q16 weights) behind a per-tile BVH; mesh export v2 carries per-vertex materials; a measurement-only
yield computes removed resource volume (0 on already-air space). With veins on, every density, mesh,
payload and topology hash is unchanged; demo off: owner hashes identical to L47. fast/medium/full PASS.
Research R18 (cavern-graph planner) and R19 (why FloatingIslands costs ~10x SurfaceWorld).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016dhxM43dnWvUu2B5Ph1iG1

Standalone/Core/src/ShapeWorldInternal.h

@@ -369,6 +369,7 @@ public:
         }
         std::sort(BoundStratumFields.begin(),BoundStratumFields.end(),
             [](const auto& A,const auto& B){return A.BottomChunkZ<B.BottomChunkZ;});
+        BuildOreVeinSpatialIndex();
     }
 
     const FBoundStratumField* FindBoundStratumForChunk(int32 ChunkZ) const
@@ -2972,6 +2973,27 @@ public:
         return AnalyticDensityAt(X,Y,Z);
     }
 
+    std::vector<std::uint32_t> OreVeinTileCandidates(
+        const std::array<std::int64_t,3>& minimum,
+        const std::array<std::int64_t,3>& maximum) const
+    {return BuildOreVeinTileList(minimum,maximum);}
+
+    VoxelForgeCore::MaterialSampleV1 EvaluateMaterialAt(float X,float Y,float Z,
+        std::span<const std::uint32_t> candidates,
+        std::uint64_t* candidateChecks=nullptr) const
+    {
+        if(SharedDensitySource)
+            return SharedDensitySource->MaterialAt(X,Y,Z,candidates,candidateChecks);
+        return MaterialAt(X,Y,Z,candidates,candidateChecks);
+    }
+
+    VoxelForgeCore::MaterialSampleV1 MaterialAt(float X,float Y,float Z) const
+    {
+        if(SharedDensitySource)return SharedDensitySource->MaterialAt(X,Y,Z);
+        const auto candidates=BuildOreVeinPointList(X,Y,Z);
+        return EvaluateMaterialAt(X,Y,Z,candidates);
+    }
+
     std::string DensityDispatchPathAt(float X,float Y,float Z) const
     {
         if(SharedDensitySource)return SharedDensitySource->DensityDispatchPathAt(X,Y,Z);
@@ -3208,6 +3230,25 @@ private:
         std::uint64_t irregularitySeed=0;
         std::uint8_t quarterTurns=0;
     };
+    struct FOreVeinSpatialBound
+    {
+        std::array<std::int64_t,3> minimumSubvoxels{},maximumSubvoxels{};
+    };
+    struct FOreVeinSpatialNode
+    {
+        std::array<std::int64_t,3> minimumSubvoxels{},maximumSubvoxels{};
+        std::size_t first=0,end=0;
+        int left=-1,right=-1;
+    };
+    struct FOreVeinRuntime
+    {
+        const VoxelForgeCore::ResolvedLayoutOreVeinV1* record=nullptr;
+        std::array<std::array<double,3>,2> endpointsVoxels{};
+        double radiusVoxels=0.0,irregularityAmplitudeVoxels=0.0;
+        double irregularityWavelengthVoxels=0.0,blendWidthVoxels=0.0;
+        std::uint64_t irregularitySeed=0;
+        std::uint16_t resourceMaterialId=0;
+    };
     struct FCavernRuntimeFeature
     {
         const VoxelForgeCore::ResolvedCavernRoomV1* room=nullptr;
@@ -3349,6 +3390,258 @@ private:
         if(!StampSpatialOrder.empty())BuildStampSpatialNode(0,StampSpatialOrder.size());
     }
 
+    int BuildOreVeinSpatialNode(std::size_t first,std::size_t end)
+    {
+        FOreVeinSpatialNode node;node.first=first;node.end=end;
+        node.minimumSubvoxels.fill(std::numeric_limits<std::int64_t>::max());
+        node.maximumSubvoxels.fill(std::numeric_limits<std::int64_t>::min());
+        for(std::size_t offset=first;offset<end;++offset)
+        {
+            const auto& bound=OreVeinSpatialBounds[OreVeinSpatialOrder[offset]];
+            for(std::size_t axis=0;axis<3;++axis)
+            {
+                node.minimumSubvoxels[axis]=std::min(node.minimumSubvoxels[axis],
+                    bound.minimumSubvoxels[axis]);
+                node.maximumSubvoxels[axis]=std::max(node.maximumSubvoxels[axis],
+                    bound.maximumSubvoxels[axis]);
+            }
+        }
+        const int nodeIndex=static_cast<int>(OreVeinSpatialNodes.size());
+        OreVeinSpatialNodes.push_back(node);
+        if(end-first>8)
+        {
+            std::size_t splitAxis=0;
+            for(std::size_t axis=1;axis<3;++axis)
+                if(static_cast<long double>(node.maximumSubvoxels[axis])
+                        -node.minimumSubvoxels[axis]
+                    >static_cast<long double>(node.maximumSubvoxels[splitAxis])
+                        -node.minimumSubvoxels[splitAxis])splitAxis=axis;
+            const std::size_t middle=first+(end-first)/2;
+            std::nth_element(OreVeinSpatialOrder.begin()+first,
+                OreVeinSpatialOrder.begin()+middle,OreVeinSpatialOrder.begin()+end,
+                [this,splitAxis](std::uint32_t a,std::uint32_t b)
+                {
+                    const auto& left=OreVeinSpatialBounds[a];
+                    const auto& right=OreVeinSpatialBounds[b];
+                    const long double leftCenter=static_cast<long double>(
+                        left.minimumSubvoxels[splitAxis])+left.maximumSubvoxels[splitAxis];
+                    const long double rightCenter=static_cast<long double>(
+                        right.minimumSubvoxels[splitAxis])+right.maximumSubvoxels[splitAxis];
+                    return leftCenter==rightCenter?a<b:leftCenter<rightCenter;
+                });
+            OreVeinSpatialNodes[nodeIndex].left=BuildOreVeinSpatialNode(first,middle);
+            OreVeinSpatialNodes[nodeIndex].right=BuildOreVeinSpatialNode(middle,end);
+        }
+        return nodeIndex;
+    }
+
+    void QueryOreVeinSpatialNode(int nodeIndex,
+        const std::array<std::int64_t,3>& minimum,
+        const std::array<std::int64_t,3>& maximum,
+        std::vector<std::uint32_t>& result) const
+    {
+        if(nodeIndex<0)return;
+        const auto& node=OreVeinSpatialNodes[static_cast<std::size_t>(nodeIndex)];
+        for(std::size_t axis=0;axis<3;++axis)
+            if(node.maximumSubvoxels[axis]<minimum[axis]
+                ||node.minimumSubvoxels[axis]>maximum[axis])return;
+        if(node.left<0)
+        {
+            for(std::size_t offset=node.first;offset<node.end;++offset)
+            {
+                const std::uint32_t index=OreVeinSpatialOrder[offset];
+                const auto& bound=OreVeinSpatialBounds[index];bool overlap=true;
+                for(std::size_t axis=0;axis<3;++axis)
+                    overlap=overlap&&bound.minimumSubvoxels[axis]<=maximum[axis]
+                        &&bound.maximumSubvoxels[axis]>=minimum[axis];
+                if(overlap)result.push_back(index);
+            }
+            return;
+        }
+        QueryOreVeinSpatialNode(node.left,minimum,maximum,result);
+        QueryOreVeinSpatialNode(node.right,minimum,maximum,result);
+    }
+
+    void BuildOreVeinSpatialIndex()
+    {
+        OreVeinSpatialBounds.clear();OreVeinSpatialOrder.clear();OreVeinSpatialNodes.clear();
+        OreVeinRuntime.clear();
+        if(!ResolvedLayout||ResolvedLayout->oreVeins.empty()||!Owner
+            ||Owner->CoreVoxelSizeMicrons==0)return;
+        constexpr long double Q=VoxelForgeCore::WorldAssemblySubvoxelsPerVoxel;
+        const long double scale=Q/Owner->CoreVoxelSizeMicrons;
+        OreVeinSpatialBounds.reserve(ResolvedLayout->oreVeins.size());
+        OreVeinRuntime.reserve(ResolvedLayout->oreVeins.size());
+        for(const auto& vein:ResolvedLayout->oreVeins)
+        {
+            FOreVeinSpatialBound bound;
+            for(std::size_t axis=0;axis<3;++axis)
+            {
+                bound.minimumSubvoxels[axis]=static_cast<std::int64_t>(std::floor(
+                    static_cast<long double>(vein.boundsMinimumMicrometres[axis])*scale));
+                bound.maximumSubvoxels[axis]=static_cast<std::int64_t>(std::ceil(
+                    static_cast<long double>(vein.boundsMaximumMicrometres[axis])*scale));
+            }
+            OreVeinSpatialBounds.push_back(bound);
+            FOreVeinRuntime runtime;runtime.record=&vein;
+            for(std::size_t endpoint=0;endpoint<2;++endpoint)
+                for(std::size_t axis=0;axis<3;++axis)
+                    runtime.endpointsVoxels[endpoint][axis]=
+                        static_cast<double>(vein.endpointsMicrometres[endpoint][axis])
+                            /Owner->CoreVoxelSizeMicrons;
+            runtime.radiusVoxels=static_cast<double>(vein.radiusMicrometres)
+                /Owner->CoreVoxelSizeMicrons;
+            runtime.irregularityAmplitudeVoxels=
+                static_cast<double>(vein.irregularityAmplitudeMicrometres)
+                    /Owner->CoreVoxelSizeMicrons;
+            runtime.irregularityWavelengthVoxels=
+                static_cast<double>(vein.irregularityWavelengthMicrometres)
+                    /Owner->CoreVoxelSizeMicrons;
+            runtime.blendWidthVoxels=static_cast<double>(vein.blendWidthMicrometres)
+                /Owner->CoreVoxelSizeMicrons;
+            runtime.irregularitySeed=vein.irregularitySeed;
+            runtime.resourceMaterialId=vein.resourceMaterialId;
+            OreVeinRuntime.push_back(runtime);
+        }
+        OreVeinSpatialOrder.resize(OreVeinSpatialBounds.size());
+        for(std::size_t i=0;i<OreVeinSpatialOrder.size();++i)
+            OreVeinSpatialOrder[i]=static_cast<std::uint32_t>(i);
+        if(!OreVeinSpatialOrder.empty())BuildOreVeinSpatialNode(0,OreVeinSpatialOrder.size());
+    }
+
+    std::vector<std::uint32_t> BuildOreVeinTileList(
+        const std::array<std::int64_t,3>& minimum,
+        const std::array<std::int64_t,3>& maximum) const
+    {
+        std::vector<std::uint32_t> result;
+        if(OreVeinSpatialBounds.empty())return result;
+        for(std::size_t axis=0;axis<3;++axis)
+            if(minimum[axis]>maximum[axis])
+            {
+                result.resize(OreVeinSpatialBounds.size());
+                for(std::size_t i=0;i<result.size();++i)
+                    result[i]=static_cast<std::uint32_t>(i);
+                return result;
+            }
+        QueryOreVeinSpatialNode(0,minimum,maximum,result);
+        std::sort(result.begin(),result.end());
+        return result;
+    }
+
+    std::vector<std::uint32_t> BuildOreVeinPointList(float X,float Y,float Z) const
+    {
+        constexpr std::int64_t Q=VoxelForgeCore::WorldAssemblySubvoxelsPerVoxel;
+        const auto toQ=[](float value)
+        {
+            const long double scaled=static_cast<long double>(value)*Q;
+            if(!std::isfinite(value)||scaled<std::numeric_limits<std::int64_t>::min()
+                ||scaled>std::numeric_limits<std::int64_t>::max())return std::int64_t{0};
+            return static_cast<std::int64_t>(std::llround(scaled));
+        };
+        const std::array<std::int64_t,3> point={toQ(X),toQ(Y),toQ(Z)};
+        return BuildOreVeinTileList(point,point);
+    }
+
+    VoxelForgeCore::MaterialSampleV1 MaterialAt(float X,float Y,float Z,
+        std::span<const std::uint32_t> candidates,
+        std::uint64_t* candidateChecks=nullptr) const
+    {
+        using namespace VoxelForgeCore;
+        MaterialSampleV1 result;
+        if(!ResolvedLayout||OreVeinRuntime.empty()||candidates.empty()
+            ||!std::isfinite(X)||!std::isfinite(Y)||!std::isfinite(Z))return result;
+        const FBoundStratumField* bound=FindBoundStratumForZ(Z);
+        if(!bound)return result;
+        const std::array<double,3> point={X,Y,Z};
+        const auto noiseAt=[](std::uint64_t seed,std::int64_t x,std::int64_t y,std::int64_t z)
+        {
+            std::uint64_t hash=14695981039346656037ull;
+            for(const std::uint64_t word:{seed,static_cast<std::uint64_t>(x),
+                static_cast<std::uint64_t>(y),static_cast<std::uint64_t>(z)})
+                for(int byte=0;byte<8;++byte)
+                {hash^=static_cast<std::uint8_t>(word>>(byte*8));hash*=1099511628211ull;}
+            hash^=hash>>30;hash*=0xbf58476d1ce4e5b9ull;
+            hash^=hash>>27;hash*=0x94d049bb133111ebull;hash^=hash>>31;
+            return static_cast<double>(static_cast<std::uint32_t>(hash>>32))/2147483647.5-1.0;
+        };
+        const auto fade=[](double value){return value*value*(3.0-2.0*value);};
+        for(const std::uint32_t candidate:candidates)
+        {
+            if(candidateChecks)++*candidateChecks;
+            if(candidate>=OreVeinRuntime.size())continue;
+            const auto& vein=OreVeinRuntime[candidate];
+            if(!vein.record||vein.record->strateIndex!=bound->StrateIndex)continue;
+            const double wavelength=vein.irregularityWavelengthVoxels;
+            double irregularity=0.0;
+            if(vein.irregularityAmplitudeVoxels>0.0&&wavelength>0.0)
+            {
+                std::array<std::int64_t,3> cell{};std::array<double,3> t{};
+                for(std::size_t axis=0;axis<3;++axis)
+                {
+                    const double coordinate=point[axis]/wavelength;
+                    cell[axis]=static_cast<std::int64_t>(std::floor(coordinate));
+                    t[axis]=fade(coordinate-cell[axis]);
+                }
+                double layer[2][2]{};
+                for(int z=0;z<2;++z)for(int y=0;y<2;++y)
+                {
+                    const double a=noiseAt(vein.irregularitySeed,cell[0],cell[1]+y,cell[2]+z);
+                    const double b=noiseAt(vein.irregularitySeed,cell[0]+1,cell[1]+y,cell[2]+z);
+                    layer[z][y]=std::lerp(a,b,t[0]);
+                }
+                irregularity=std::lerp(std::lerp(layer[0][0],layer[0][1],t[1]),
+                    std::lerp(layer[1][0],layer[1][1],t[1]),t[2]);
+            }
+            const auto& a=vein.endpointsVoxels[0];const auto& b=vein.endpointsVoxels[1];
+            double lengthSquared=0.0,dot=0.0;
+            for(std::size_t axis=0;axis<3;++axis)
+            {const double delta=b[axis]-a[axis];lengthSquared+=delta*delta;dot+=(point[axis]-a[axis])*delta;}
+            const double t=lengthSquared>0.0?std::clamp(dot/lengthSquared,0.0,1.0):0.0;
+            double distanceSquared=0.0;
+            for(std::size_t axis=0;axis<3;++axis)
+            {const double delta=point[axis]-(a[axis]+t*(b[axis]-a[axis]));distanceSquared+=delta*delta;}
+            const double distance=std::sqrt(distanceSquared);
+            const double radius=std::max(0.0,vein.radiusVoxels
+                +irregularity*vein.irregularityAmplitudeVoxels);
+            if(distance>=radius+vein.blendWidthVoxels)continue;
+            double weight=1.0;
+            if(distance>radius)
+            {
+                const double u=std::clamp((radius+vein.blendWidthVoxels-distance)
+                    /vein.blendWidthVoxels,0.0,1.0);
+                weight=u*u*(3.0-2.0*u);
+            }
+            const auto weightQ16=static_cast<std::uint16_t>(std::clamp(
+                std::llround(weight*65535.0),0ll,65535ll));
+            if(weightQ16==0)continue;
+            std::size_t slot=2;
+            for(std::size_t index=0;index<2;++index)
+                if(result.resourceMaterialIds[index]==vein.resourceMaterialId)
+                {slot=index;break;}
+            if(slot==2)
+                for(std::size_t index=0;index<2;++index)
+                    if(result.resourceMaterialIds[index]==0){slot=index;break;}
+            if(slot<2)
+            {
+                result.resourceMaterialIds[slot]=vein.resourceMaterialId;
+                const std::uint32_t sum=static_cast<std::uint32_t>(
+                    result.resourceWeightsQ16[slot])+weightQ16;
+                result.resourceWeightsQ16[slot]=static_cast<std::uint16_t>(
+                    std::min<std::uint32_t>(65535,sum));
+            }
+        }
+        const std::uint32_t total=static_cast<std::uint32_t>(result.resourceWeightsQ16[0])
+            +result.resourceWeightsQ16[1];
+        if(total>65535)
+        {
+            result.resourceWeightsQ16[0]=static_cast<std::uint16_t>(
+                static_cast<std::uint64_t>(result.resourceWeightsQ16[0])*65535/total);
+            result.resourceWeightsQ16[1]=static_cast<std::uint16_t>(
+                65535-result.resourceWeightsQ16[0]);
+        }
+        return result;
+    }
+
     std::vector<std::uint32_t> BuildStampTileList(
         const std::array<std::int64_t,3>& minimum,
         const std::array<std::int64_t,3>& maximum) const
@@ -9412,6 +9705,10 @@ private:
     std::vector<FStampSpatialNode> StampSpatialNodes;
     std::vector<std::uint32_t> StampSpatialOrder;
     std::vector<FStampRuntime> StampRuntime;
+    std::vector<FOreVeinSpatialBound> OreVeinSpatialBounds;
+    std::vector<FOreVeinSpatialNode> OreVeinSpatialNodes;
+    std::vector<std::uint32_t> OreVeinSpatialOrder;
+    std::vector<FOreVeinRuntime> OreVeinRuntime;
     std::vector<FCavernRuntimeFeature> CavernRuntime;
     std::vector<VoxelForgeCore::WorldAssemblyFeatureBoundsV1> PassageFeatureBounds;
     std::vector<FPassageTileFeature> PassageTileFeatures;