VoxelForge/Standalone/Research/WorldAssemblerM3/Round1/FloatingIslandsAudit.cpp

#include "CoreMinimal.h"
#include "OwnerWorldLoader.h"
#include "VoxelForgeCore/VoxelForgeCore.h"

#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <filesystem>
#include <fstream>
#include <functional>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <limits>
#include <numeric>
#include <sstream>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>

#if defined(_WIN32)
#ifdef TEXT
#undef TEXT
#endif
#ifndef NOMINMAX
#define NOMINMAX
#endif
#include <Windows.h>
#ifdef TEXT
#undef TEXT
#endif
#endif

namespace
{
    using namespace VoxelForgeCore;
    constexpr std::int64_t Q = WorldAssemblySubvoxelsPerVoxel;
    using Clock = std::chrono::steady_clock;
    using Point = std::array<std::int64_t, 3>;

    struct OracleAnchor
    {
        Point surfaceQ{};
        std::uint32_t supportScore = 0;
    };

    struct SampleStats
    {
        std::string hash;
        double legacyMilliseconds = 0.0;
        double nativeMilliseconds = 0.0;
        std::uint64_t signMismatches = 0;
        std::uint64_t legacySolidSamples = 0;
        std::uint64_t nativeSolidSamples = 0;
        std::uint64_t legacyOnlySolidSamples = 0;
        std::uint64_t nativeOnlySolidSamples = 0;
        double legacyChecksum = 0.0;
        double nativeChecksum = 0.0;
    };

    struct Triangle
    {
        std::array<float, 12> values{}; // Three xyz vertices followed by one flat normal.
    };

    struct MeshResult
    {
        std::uint64_t gridSamples = 0;
        std::uint64_t triangles = 0;
        std::uint64_t bytes = 0;
    };

    struct ExistingHashes
    {
        std::string round1Manifest;
        std::string surfaceWorldRound4Manifest;
    };

    std::uint64_t Mix(std::uint64_t value)
    {
        value += 0x9e3779b97f4a7c15ull;
        value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ull;
        value = (value ^ (value >> 27)) * 0x94d049bb133111ebull;
        return value ^ (value >> 31);
    }

    std::uint64_t Next(std::uint64_t& state)
    {
        state += 0x9e3779b97f4a7c15ull;
        return Mix(state);
    }

    std::int64_t VoxQ(double value)
    { return static_cast<std::int64_t>(std::llround(value * static_cast<double>(Q))); }

    std::string Hex64(std::uint64_t value)
    {
        std::ostringstream output;
        output << std::hex << std::setfill('0') << std::setw(16) << value;
        return output.str();
    }

    std::string HashValues(const std::vector<std::int64_t>& values)
    {
        std::uint64_t hash = 14695981039346656037ull;
        for (const auto value : values)
            for (int byte = 0; byte < 8; ++byte)
            {
                hash ^= static_cast<std::uint8_t>(static_cast<std::uint64_t>(value) >> (byte * 8));
                hash *= 1099511628211ull;
            }
        return Hex64(hash);
    }

    ExistingHashes VerifyExistingManifests(const std::filesystem::path& standaloneRoot)
    {
        const auto read = [](const std::filesystem::path& path)
        {
            std::ifstream input(path, std::ios::binary);
            if (!input) throw std::runtime_error("cannot read existing manifest: " + path.string());
            return std::string(std::istreambuf_iterator<char>(input), std::istreambuf_iterator<char>());
        };
        ExistingHashes hashes;
        std::string error;
        const auto round1 = WorldAssemblyManifest::FromJson(read(
            standaloneRoot / "Research/WorldAssemblerM1/round1-manifest.json"), error);
        if (!round1 || !round1->validate(error)
            || Hex64(round1->contentHash) != "c87ae279b0e26f53")
            throw std::runtime_error("existing M1 manifest hash changed or failed validation: " + error);
        hashes.round1Manifest = Hex64(round1->contentHash);
        error.clear();
        const auto surface = WorldAssemblyManifest::FromJson(read(
            standaloneRoot / "Research/WorldAssemblerM1/round4/world_assembly.json"), error);
        if (!surface || !surface->validate(error)
            || Hex64(surface->contentHash) != "8933693d9560ba9f")
            throw std::runtime_error("existing SurfaceWorld round-4 hash changed or failed validation: " + error);
        hashes.surfaceWorldRound4Manifest = Hex64(surface->contentHash);
        return hashes;
    }

    float LegacyDensity(const FStandaloneOwnerWorld& owner, const Point& q)
    {
        constexpr double CoreToLegacyVoxelScale = 2.0;
        return owner.Generator->GetDensityAt(
            static_cast<float>(CoreToLegacyVoxelScale * q[0] / static_cast<double>(Q)),
            static_cast<float>(CoreToLegacyVoxelScale * q[1] / static_cast<double>(Q)),
            static_cast<float>(CoreToLegacyVoxelScale * q[2] / static_cast<double>(Q)));
    }

    OracleAnchor FindLegacyAnchor(const FStandaloneOwnerWorld& owner)
    {
        OracleAnchor best;
        const float highZ = -560.0f;
        const float lowZ = -2480.0f;
        constexpr float ScanStep = 8.0f;
        for (int y = -448; y <= 448; y += 16)
        for (int x = -448; x <= 448; x += 16)
        {
            float airZ = highZ;
            float above = owner.Generator->GetDensityAt(static_cast<float>(x),
                static_cast<float>(y), airZ);
            if (above < 0.0f) continue;
            for (float solidZ = highZ - ScanStep; solidZ >= lowZ; solidZ -= ScanStep)
            {
                const float density = owner.Generator->GetDensityAt(static_cast<float>(x),
                    static_cast<float>(y), solidZ);
                if (density < 0.0f)
                {
                    float outside = airZ, inside = solidZ;
                    for (int iteration = 0; iteration < 12; ++iteration)
                    {
                        const float middle = (outside + inside) * 0.5f;
                        if (owner.Generator->GetDensityAt(static_cast<float>(x),
                                static_cast<float>(y), middle) < 0.0f) inside = middle;
                        else outside = middle;
                    }
                    const std::int64_t topQ = VoxQ(outside);
                    std::uint32_t score = 0;
                    for (int dy = -16; dy <= 16; dy += 8)
                    for (int dx = -16; dx <= 16; dx += 8)
                        if (owner.Generator->GetDensityAt(static_cast<float>(x + dx),
                                static_cast<float>(y + dy),
                                static_cast<float>((topQ - 4 * Q) / static_cast<double>(Q))) < 0.0f)
                            ++score;
                    if (score > best.supportScore
                        || (score == best.supportScore
                            && std::abs(x) + std::abs(y)
                                < std::abs(best.surfaceQ[0] / Q) + std::abs(best.surfaceQ[1] / Q)))
                    {
                        best.surfaceQ = {static_cast<std::int64_t>(x) * Q,
                            static_cast<std::int64_t>(y) * Q, topQ};
                        best.supportScore = score;
                    }
                    break; // The first air-to-solid crossing is the island top.
                }
                airZ = solidZ;
            }
        }
        return best;
    }

    std::vector<Point> MakeSamples(const FloatingIslandInfluenceBounds& bounds,
                                   std::size_t count, std::uint64_t seed)
    {
        std::vector<Point> points(count);
        std::uint64_t state = seed;
        for (auto& point : points)
            for (std::size_t axis = 0; axis < 3; ++axis)
            {
                const std::uint64_t span = static_cast<std::uint64_t>(
                    bounds.maximumSubvoxels[axis] - bounds.minimumSubvoxels[axis] + 1);
                point[axis] = bounds.minimumSubvoxels[axis]
                    + static_cast<std::int64_t>(Next(state) % std::max<std::uint64_t>(span, 1));
            }
        return points;
    }

    SampleStats HashWithWorkers(const FloatingIslandFieldV1& field,
        const FStandaloneOwnerWorld& owner, const std::vector<Point>& points,
        std::uint32_t workers, bool measure)
    {
        std::vector<std::int64_t> values(points.size());
        std::vector<std::thread> threads;
        threads.reserve(workers);
        const auto evaluate = [&](std::uint32_t worker)
        {
            for (std::size_t i = worker; i < points.size(); i += workers)
                values[i] = field.densitySubvoxels(points[i][0], points[i][1], points[i][2]);
        };
        for (std::uint32_t worker = 0; worker < workers; ++worker)
            threads.emplace_back(evaluate, worker);
        for (auto& thread : threads) thread.join();

        SampleStats stats;
        stats.hash = HashValues(values);
        if (!measure) return stats;
        std::vector<float> legacyValues(points.size());
        const auto beginLegacy = Clock::now();
        double legacyChecksum = 0.0;
        for (std::size_t i = 0; i < points.size(); ++i)
        {
            const float value = LegacyDensity(owner, points[i]);
            legacyValues[i] = value;
            legacyChecksum += value;
        }
        const auto endLegacy = Clock::now();
        const auto beginNative = Clock::now();
        double nativeChecksum = 0.0;
        for (const Point& point : points)
            nativeChecksum += static_cast<double>(field.densitySubvoxels(
                point[0], point[1], point[2])) / static_cast<double>(Q);
        const auto endNative = Clock::now();
        for (std::size_t i = 0; i < values.size(); ++i)
        {
            const bool legacySolid = legacyValues[i] < 0.0f;
            const bool nativeSolid = values[i] < 0;
            if (legacySolid) ++stats.legacySolidSamples;
            if (nativeSolid) ++stats.nativeSolidSamples;
            if (legacySolid && !nativeSolid) ++stats.legacyOnlySolidSamples;
            if (nativeSolid && !legacySolid) ++stats.nativeOnlySolidSamples;
            if (legacySolid != nativeSolid) ++stats.signMismatches;
        }
        stats.legacyMilliseconds = std::chrono::duration<double, std::milli>(endLegacy - beginLegacy).count();
        stats.nativeMilliseconds = std::chrono::duration<double, std::milli>(endNative - beginNative).count();
        stats.legacyChecksum = legacyChecksum;
        stats.nativeChecksum = nativeChecksum;
        return stats;
    }

    using D3 = std::array<double, 3>;
    D3 Add(D3 a, D3 b) { return {a[0] + b[0], a[1] + b[1], a[2] + b[2]}; }
    D3 Sub(D3 a, D3 b) { return {a[0] - b[0], a[1] - b[1], a[2] - b[2]}; }
    D3 Mul(D3 a, double s) { return {a[0] * s, a[1] * s, a[2] * s}; }
    double Dot(D3 a, D3 b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; }
    D3 Cross(D3 a, D3 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]}; }
    double Length(D3 a) { return std::sqrt(Dot(a, a)); }
    D3 Normalize(D3 a)
    {
        const double length = Length(a);
        return length > 1e-12 ? Mul(a, 1.0 / length) : D3{0.0, 0.0, 1.0};
    }

    double SampleNative(const FloatingIslandFieldV1& field, const Point& q)
    { return static_cast<double>(field.densitySubvoxels(q[0], q[1], q[2])) / static_cast<double>(Q); }

    double SampleLegacy(const FStandaloneOwnerWorld& owner, const Point& q)
    { return static_cast<double>(LegacyDensity(owner, q)); }

    double TetGradient(const std::array<D3, 4>& points, const std::array<double, 4>& values, D3& output)
    {
        const D3 a = Sub(points[1], points[0]);
        const D3 b = Sub(points[2], points[0]);
        const D3 c = Sub(points[3], points[0]);
        const D3 bc = Cross(b, c), ca = Cross(c, a), ab = Cross(a, b);
        const double determinant = Dot(a, bc);
        if (std::abs(determinant) < 1e-12) { output = {0.0, 0.0, 1.0}; return 0.0; }
        output = Mul(Add(Add(Mul(bc, values[1] - values[0]),
            Mul(ca, values[2] - values[0])), Mul(ab, values[3] - values[0])), 1.0 / determinant);
        return Length(output);
    }

    MeshResult ExtractMesh(const std::filesystem::path& outputPath,
        Point minimumQ, Point maximumQ, std::int64_t stepQ,
        const std::function<double(const Point&)>& density)
    {
        const std::size_t nx = static_cast<std::size_t>((maximumQ[0] - minimumQ[0]) / stepQ + 1);
        const std::size_t ny = static_cast<std::size_t>((maximumQ[1] - minimumQ[1]) / stepQ + 1);
        const std::size_t nz = static_cast<std::size_t>((maximumQ[2] - minimumQ[2]) / stepQ + 1);
        const std::size_t gridCount = nx * ny * nz;
        if (gridCount > 4'000'000) throw std::runtime_error("mesh grid exceeds the 4M-sample cap");
        std::vector<float> grid(gridCount);
        const auto gridIndex = [=](std::size_t x, std::size_t y, std::size_t z)
        { return (z * ny + y) * nx + x; };
        for (std::size_t z = 0; z < nz; ++z)
        for (std::size_t y = 0; y < ny; ++y)
        for (std::size_t x = 0; x < nx; ++x)
        {
            const Point q{minimumQ[0] + static_cast<std::int64_t>(x) * stepQ,
                minimumQ[1] + static_cast<std::int64_t>(y) * stepQ,
                minimumQ[2] + static_cast<std::int64_t>(z) * stepQ};
            grid[gridIndex(x, y, z)] = static_cast<float>(density(q));
        }

        constexpr std::array<std::array<int, 3>, 8> Corner = {{{0,0,0},{1,0,0},{1,1,0},{0,1,0},
            {0,0,1},{1,0,1},{1,1,1},{0,1,1}}};
        constexpr std::array<std::array<int, 4>, 6> Tets = {{{0,5,1,6},{0,1,2,6},{0,2,3,6},
            {0,3,7,6},{0,7,4,6},{0,4,5,6}}};
        constexpr std::array<std::array<int, 2>, 6> TetEdges = {{{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}}};
        std::vector<Triangle> triangles;
        triangles.reserve(8192);
        for (std::size_t z = 0; z + 1 < nz; ++z)
        for (std::size_t y = 0; y + 1 < ny; ++y)
        for (std::size_t x = 0; x + 1 < nx; ++x)
        {
            std::array<D3, 8> cubePoints{};
            std::array<double, 8> cubeValues{};
            double low = std::numeric_limits<double>::infinity();
            double high = -std::numeric_limits<double>::infinity();
            for (std::size_t corner = 0; corner < Corner.size(); ++corner)
            {
                const auto& offset = Corner[corner];
                cubePoints[corner] = {
                    static_cast<double>(minimumQ[0] + static_cast<std::int64_t>(x + offset[0]) * stepQ) / Q,
                    static_cast<double>(minimumQ[1] + static_cast<std::int64_t>(y + offset[1]) * stepQ) / Q,
                    static_cast<double>(minimumQ[2] + static_cast<std::int64_t>(z + offset[2]) * stepQ) / Q};
                cubeValues[corner] = grid[gridIndex(x + offset[0], y + offset[1], z + offset[2])];
                low = std::min(low, cubeValues[corner]);
                high = std::max(high, cubeValues[corner]);
            }
            if (low >= 0.0 || high < 0.0) continue;

            for (const auto& tet : Tets)
            {
                std::array<D3, 4> points{};
                std::array<double, 4> values{};
                for (std::size_t i = 0; i < 4; ++i)
                { points[i] = cubePoints[tet[i]]; values[i] = cubeValues[tet[i]]; }
                std::array<D3, 4> polygon{};
                std::size_t polygonCount = 0;
                for (const auto& edge : TetEdges)
                {
                    const int ia = edge[0], ib = edge[1];
                    if ((values[ia] < 0.0) == (values[ib] < 0.0)) continue;
                    const double t = values[ia] / (values[ia] - values[ib]);
                    const D3 point = Add(points[ia], Mul(Sub(points[ib], points[ia]), t));
                    bool duplicate = false;
                    for (std::size_t i = 0; i < polygonCount; ++i)
                        if (Length(Sub(polygon[i], point)) < 1e-7) { duplicate = true; break; }
                    if (!duplicate && polygonCount < polygon.size()) polygon[polygonCount++] = point;
                }
                if (polygonCount < 3) continue;
                D3 normal{};
                TetGradient(points, values, normal);
                normal = Normalize(normal);
                D3 center{};
                for (std::size_t i = 0; i < polygonCount; ++i) center = Add(center, polygon[i]);
                center = Mul(center, 1.0 / polygonCount);
                const D3 axis = std::abs(normal[2]) < 0.85 ? D3{0,0,1} : D3{0,1,0};
                const D3 right = Normalize(Cross(axis, normal));
                const D3 up = Normalize(Cross(normal, right));
                std::array<std::pair<double, D3>, 4> ordered{};
                for (std::size_t i = 0; i < polygonCount; ++i)
                {
                    const D3 delta = Sub(polygon[i], center);
                    ordered[i] = {std::atan2(Dot(delta, up), Dot(delta, right)), polygon[i]};
                }
                std::sort(ordered.begin(), ordered.begin() + polygonCount,
                    [](const auto& a, const auto& b) { return a.first < b.first; });
                for (std::size_t i = 1; i + 1 < polygonCount; ++i)
                {
                    D3 a = ordered[0].second, b = ordered[i].second, c = ordered[i + 1].second;
                    if (Dot(Cross(Sub(b, a), Sub(c, a)), normal) < 0.0) std::swap(b, c);
                    Triangle triangle;
                    for (std::size_t axisIndex = 0; axisIndex < 3; ++axisIndex)
                    {
                        triangle.values[axisIndex] = static_cast<float>(a[axisIndex]);
                        triangle.values[3 + axisIndex] = static_cast<float>(b[axisIndex]);
                        triangle.values[6 + axisIndex] = static_cast<float>(c[axisIndex]);
                        triangle.values[9 + axisIndex] = static_cast<float>(normal[axisIndex]);
                    }
                    triangles.push_back(triangle);
                }
            }
        }

        std::ofstream output(outputPath, std::ios::binary | std::ios::trunc);
        if (!output) throw std::runtime_error("cannot create mesh capture: " + outputPath.string());
        const std::uint32_t count = static_cast<std::uint32_t>(triangles.size());
        output.write("VFMI", 4);
        output.write(reinterpret_cast<const char*>(&count), sizeof(count));
        for (const auto& triangle : triangles)
            output.write(reinterpret_cast<const char*>(triangle.values.data()),
                static_cast<std::streamsize>(sizeof(triangle.values)));
        output.close();
        MeshResult result;
        result.gridSamples = gridCount;
        result.triangles = triangles.size();
        result.bytes = std::filesystem::file_size(outputPath);
        return result;
    }

    void WriteJson(const std::filesystem::path& path, const FloatingIslandManifestV1& manifest,
        const FloatingIslandBoundsAuditResult& bounds, const FloatingIslandRouteAuditResult& route,
        const FloatingIslandRouteAuditResult& raised, const SampleStats& sample,
        const FloatingIslandThicknessAuditResult& thickness,
        const FloatingIslandThicknessAuditResult& thinThickness,
        const OracleAnchor& anchor, const std::array<MeshResult, 4>& meshes,
        const std::array<std::string, 3>& fieldHashes, const ExistingHashes& oldHashes)
    {
        const Point ownerAnchorQ = anchor.surfaceQ;
        const Point coreAnchorQ{ownerAnchorQ[0] / 2, ownerAnchorQ[1] / 2, ownerAnchorQ[2] / 2};
        std::ofstream out(path, std::ios::trunc);
        if (!out) throw std::runtime_error("cannot create audit result JSON");
        out << std::setprecision(12);
        out << "{\n  \"recipe_hash\":\"" << Hex64(manifest.recipe.contentHash)
            << "\",\n  \"plan_hash\":\"" << Hex64(manifest.plan.contentHash)
            << "\",\n  \"manifest_hash\":\"" << Hex64(manifest.contentHash)
            << "\",\n  \"unchanged_existing_hashes\":{\"m1_round1_manifest\":\""
            << oldHashes.round1Manifest << "\",\"surfaceworld_round4_manifest\":\""
            << oldHashes.surfaceWorldRound4Manifest << "\"}"
            << ",\n  \"geometry_version\":" << manifest.recipe.geometryVersion
            << ",\n  \"owner_voxel_size_cm\":12.5,\n  \"core_voxel_size_cm\":25"
            << ",\n  \"owner_to_core_length_scale\":0.5"
            << ",\n  \"plan_islands\":" << manifest.plan.islands.size()
            << ",\n  \"plan_connectors\":" << manifest.plan.clusterEdges.size()
            << ",\n  \"route_points\":" << manifest.plan.route.size()
            << ",\n  \"legacy_anchor_owner_voxels\":[" << ownerAnchorQ[0] / static_cast<double>(Q)
            << ',' << ownerAnchorQ[1] / static_cast<double>(Q) << ','
            << ownerAnchorQ[2] / static_cast<double>(Q) << ']'
            << ",\n  \"comparison_anchor_core_voxels\":[" << coreAnchorQ[0] / static_cast<double>(Q)
            << ',' << coreAnchorQ[1] / static_cast<double>(Q) << ','
            << coreAnchorQ[2] / static_cast<double>(Q) << "],\n  \"legacy_anchor_plateau_samples\":"
            << anchor.supportScore << ",\n  \"recipe_parameters\":{\"seed\":" << manifest.recipe.seed
            << ",\"primary_island_count\":" << manifest.recipe.primaryIslandCount
            << ",\"density_q16\":" << manifest.recipe.islandDensityQ16
            << ",\"spacing_voxels\":" << manifest.recipe.spacingSubvoxels / static_cast<double>(Q)
            << ",\"radius_voxels\":[" << manifest.recipe.minimumRadiusSubvoxels / static_cast<double>(Q)
            << ',' << manifest.recipe.maximumRadiusSubvoxels / static_cast<double>(Q)
            << "],\"thickness_ratio_q16\":" << manifest.recipe.thicknessRatioQ16
            << ",\"vertical_jitter_q16\":" << manifest.recipe.verticalJitterQ16
            << ",\"top_flatten_q16\":" << manifest.recipe.topFlattenQ16
            << ",\"contour_warp_voxels\":" << manifest.recipe.contourWarpSubvoxels / static_cast<double>(Q)
            << ",\"crust_depth_bonus_q16\":" << manifest.recipe.crustDepthBonusQ16
            << ",\"root_depth_ratio_q16\":" << manifest.recipe.rootDepthRatioQ16
            << ",\"roughness_voxels\":" << manifest.recipe.roughnessSubvoxels / static_cast<double>(Q)
            << ",\"blend_voxels\":" << manifest.recipe.blendRadiusSubvoxels / static_cast<double>(Q)
            << ",\"bridge_half_width_voxels\":" << manifest.recipe.bridgeHalfWidthSubvoxels / static_cast<double>(Q)
            << ",\"bridge_thickness_voxels\":" << manifest.recipe.bridgeThicknessSubvoxels / static_cast<double>(Q)
            << "},\n  \"stable_island_ids\":[";
        for (std::size_t i = 0; i < manifest.plan.islands.size(); ++i)
        { if (i) out << ','; out << '\"' << Hex64(manifest.plan.islands[i].stableId) << '\"'; }
        out << "],\n  \"stable_connector_ids\":[";
        for (std::size_t i = 0; i < manifest.plan.clusterEdges.size(); ++i)
        { if (i) out << ','; out << '\"' << Hex64(manifest.plan.clusterEdges[i].stableId) << '\"'; }
        out << "],\n  \"island_profiles\":[";
        for (std::size_t i = 0; i < manifest.plan.islands.size(); ++i)
        {
            if (i) out << ',';
            const auto& island = manifest.plan.islands[i];
            const double ratio = static_cast<double>(island.topZSubvoxels - island.bottomZSubvoxels)
                / (2.0 * static_cast<double>(island.radiusSubvoxels));
            out << "{\"id\":\"" << Hex64(island.stableId) << "\",\"ordinal\":" << island.ordinal
                << ",\"center_voxels\":[" << island.centerSubvoxels[0] / static_cast<double>(Q)
                << ',' << island.centerSubvoxels[1] / static_cast<double>(Q)
                << ',' << island.centerSubvoxels[2] / static_cast<double>(Q)
                << "],\"radius_voxels\":" << island.radiusSubvoxels / static_cast<double>(Q)
                << ",\"top_voxels\":" << island.topZSubvoxels / static_cast<double>(Q)
                << ",\"bottom_voxels\":" << island.bottomZSubvoxels / static_cast<double>(Q)
                << ",\"depth_to_diameter\":" << ratio << '}';
        }
        out << "],\n  \"bounds\":{\"passed\":" << (bounds.passed ? "true" : "false")
            << ",\"negative_control_detected\":" << (bounds.negativeControlDetected ? "true" : "false")
            << ",\"boxes\":" << bounds.boxesTested << ",\"samples\":" << bounds.pointSamples
            << ",\"unknown\":" << bounds.unknownBounds << ",\"missed\":" << bounds.missedBounds
            << ",\"island_boxes\":" << bounds.islandBoundsTested
            << ",\"connector_boxes\":" << bounds.connectorBoundsTested
            << ",\"void_boxes_proved_air\":" << bounds.voidBoxesProvedAir << "},\n  \"thickness\":{\"passed\":"
            << (thickness.passed ? "true" : "false")
            << ",\"islands_tested\":" << thickness.islandsTested
            << ",\"islands_below_minimum\":" << thickness.islandsBelowMinimum
            << ",\"minimum_depth_to_diameter_ratio\":" << thickness.minimumDepthToDiameterRatio
            << ",\"required_depth_to_diameter_ratio\":" << thickness.requiredDepthToDiameterRatio
            << ",\"minimum_island_id\":\"" << Hex64(thickness.minimumIslandStableId) << "\"}"
            << ",\n  \"thin_pad_negative_control\":{\"detected\":"
            << ((!thinThickness.passed && thinThickness.islandsBelowMinimum == thinThickness.islandsTested)
                ? "true" : "false")
            << ",\"failed_islands\":" << thinThickness.islandsBelowMinimum
            << ",\"minimum_depth_to_diameter_ratio\":" << thinThickness.minimumDepthToDiameterRatio << "}"
            << ",\n  \"route\":{\"passed\":"
            << (route.passed ? "true" : "false") << ",\"segments\":" << route.segmentsTested
            << ",\"centers\":" << route.centerSamples << ",\"capsule_samples\":"
            << route.capsuleFieldSamples << ",\"blocked\":" << route.blockedCapsuleSamples
            << ",\"support_samples\":" << route.supportProbeSamples << ",\"unsupported\":"
            << route.unsupportedCenterSamples << ",\"grade_failures\":" << route.gradeFailures
            << ",\"minimum_clearance_subvoxels\":" << route.minimumClearanceSubvoxels
            << "},\n  \"raised_route_negative_control\":{\"detected\":"
            << ((!raised.passed && (raised.unsupportedCenterSamples || raised.blockedCapsuleSamples)) ? "true" : "false")
            << ",\"passed\":" << (raised.passed ? "true" : "false")
            << ",\"unsupported\":" << raised.unsupportedCenterSamples
            << ",\"blocked\":" << raised.blockedCapsuleSamples << "},\n  \"determinism\":{\"W1\":\""
            << fieldHashes[0] << "\",\"W8\":\"" << fieldHashes[1] << "\",\"W16\":\""
            << fieldHashes[2] << "\",\"plan_repeat\":true,\"hash_samples\":32768},\n  \"timing\":{\"same_queries\":32768,\"legacy_ms\":"
            << sample.legacyMilliseconds << ",\"native_ms\":" << sample.nativeMilliseconds
            << ",\"native_over_legacy\":" << (sample.nativeMilliseconds / std::max(sample.legacyMilliseconds, 1e-9))
            << ",\"sign_mismatches\":" << sample.signMismatches << ",\"sign_mismatch_percent\":"
            << (100.0 * sample.signMismatches / 32768.0)
            << ",\"legacy_solid_samples\":" << sample.legacySolidSamples
            << ",\"native_solid_samples\":" << sample.nativeSolidSamples
            << ",\"legacy_only_solid_samples\":" << sample.legacyOnlySolidSamples
            << ",\"native_only_solid_samples\":" << sample.nativeOnlySolidSamples << "},\n  \"mesh\":[";
        constexpr std::array<const char*, 4> names{"cluster_legacy", "cluster_native", "detail_legacy", "detail_native"};
        for (std::size_t i = 0; i < meshes.size(); ++i)
        {
            if (i) out << ',';
            out << "{\"name\":\"" << names[i] << "\",\"grid_samples\":" << meshes[i].gridSamples
                << ",\"triangles\":" << meshes[i].triangles << ",\"bytes\":" << meshes[i].bytes << '}';
        }
        out << "],\n  \"route_points_core_voxels\":[";
        for (std::size_t i = 0; i < manifest.plan.route.size(); ++i)
        {
            if (i) out << ',';
            const auto& point = manifest.plan.route[i];
            out << '[' << point.positionSubvoxels[0] / static_cast<double>(Q) << ','
                << point.positionSubvoxels[1] / static_cast<double>(Q) << ','
                << point.positionSubvoxels[2] / static_cast<double>(Q) << ']';
        }
        out << "]\n}\n";
    }
}

int main(int argc, char** argv)
{
#if defined(_WIN32)
    SetPriorityClass(GetCurrentProcess(), BELOW_NORMAL_PRIORITY_CLASS);
    SetProcessAffinityMask(GetCurrentProcess(), static_cast<DWORD_PTR>(0xFFFF));
#endif
    try
    {
        if (argc != 3)
        { std::cerr << "usage: FloatingIslandsAudit <legacy-probe.json> <output-dir>\n"; return 2; }
        const std::filesystem::path probePath = argv[1];
        const std::filesystem::path outputDir = argv[2];
        std::filesystem::create_directories(outputDir);
        const auto standaloneRoot = probePath.parent_path().parent_path().parent_path()
            .parent_path().parent_path();
        const ExistingHashes oldHashes = VerifyExistingManifests(standaloneRoot);

        FStandaloneOwnerWorld owner;
        std::string error;
        if (!StandaloneOwnerLoader::Load(probePath.string(), owner, error))
            throw std::runtime_error("legacy probe load failed: " + error);
        const OracleAnchor anchor = FindLegacyAnchor(owner);
        if (anchor.supportScore < 9)
            throw std::runtime_error("legacy scan found no broad FloatingIslands top in the frozen strate range");
        const Point coreAnchorQ{anchor.surfaceQ[0] / 2,
            anchor.surfaceQ[1] / 2, anchor.surfaceQ[2] / 2};

        FloatingIslandRecipeV1 recipe;
        recipe.seed = 0;
        recipe.layoutOriginXSubvoxels = coreAnchorQ[0];
        recipe.layoutOriginYSubvoxels = coreAnchorQ[1];
        recipe.layoutCenterTopZSubvoxels = coreAnchorQ[2];
        recipe.geometryVersion = 2;
        recipe.layoutMinimumZSubvoxels = coreAnchorQ[2] - 80 * Q;
        recipe.layoutMaximumZSubvoxels = coreAnchorQ[2] + 48 * Q;
        recipe.primaryIslandCount = 3;
        recipe.islandDensityQ16 = 32768;
        recipe.spacingSubvoxels = 56 * Q;
        recipe.minimumRadiusSubvoxels = 12 * Q;
        recipe.maximumRadiusSubvoxels = 24 * Q;
        recipe.thicknessRatioQ16 = 39322;
        recipe.verticalJitterQ16 = 39322;
        recipe.topFlattenQ16 = 39322;
        recipe.taperEndMinimumQ16 = 29491;
        recipe.taperEndMaximumQ16 = 45875;
        recipe.blendRadiusSubvoxels = 3 * Q;
        recipe.contourWarpSubvoxels = VoxQ(6.3);
        recipe.roughnessSubvoxels = VoxQ(1.5);
        recipe.boundarySealSubvoxels = 2 * Q;
        recipe.bridgeHalfWidthSubvoxels = 6 * Q;
        recipe.bridgeThicknessSubvoxels = 4 * Q;

        const auto beginPlanning = Clock::now();
        auto planned = PlanFloatingIslandCluster(recipe, error);
        const auto endPlanning = Clock::now();
        if (!planned) throw std::runtime_error("native island planner failed: " + error);
        const std::string canonical = planned->ToJson();
        for (int repeat = 0; repeat < 2; ++repeat)
        {
            std::string repeatError;
            const auto repeated = PlanFloatingIslandCluster(recipe, repeatError);
            if (!repeated || repeated->contentHash != planned->contentHash
                || repeated->plan.contentHash != planned->plan.contentHash
                || repeated->ToJson() != canonical)
                throw std::runtime_error("repeated native plan hash or canonical JSON changed");
        }
        const auto parsed = FloatingIslandManifestV1::FromJson(canonical, error);
        if (!parsed || parsed->ToJson() != canonical || !parsed->validate(error))
            throw std::runtime_error("native manifest did not validate after round-trip: " + error);

        FloatingIslandFieldV1 field(planned->recipe, planned->plan);
        if (!field.isValid()) throw std::runtime_error("native field rejected the planned island graph");
        const auto bounds = AuditFloatingIslandBounds(field, 512);
        const auto route = AuditFloatingIslandRoute(field);
        const auto raisedRoute = AuditFloatingIslandRoute(field, 8 * Q);
        const auto thickness = AuditFloatingIslandThickness(planned->plan);
        auto thinPlan = planned->plan;
        for (auto& island : thinPlan.islands)
            island.bottomZSubvoxels = island.topZSubvoxels
                - island.radiusSubvoxels * planned->recipe.thicknessRatioQ16 / 65536;
        const auto thinThickness = AuditFloatingIslandThickness(thinPlan);
        const auto influence = field.influenceBounds();
        const auto points = MakeSamples(influence, 32768, 0x4d33524f554e4431ull);
        const auto w1 = HashWithWorkers(field, owner, points, 1, false);
        const auto w8 = HashWithWorkers(field, owner, points, 8, false);
        const auto w16 = HashWithWorkers(field, owner, points, 16, true);
        const std::array<std::string, 3> hashes{w1.hash, w8.hash, w16.hash};
        if (w1.hash != w8.hash || w1.hash != w16.hash)
            throw std::runtime_error("native field sample hash differs between W1/W8/W16");
        if (!bounds.passed) throw std::runtime_error("native island bounds audit or negative control failed");
        if (!thickness.passed || thinThickness.passed
            || thinThickness.islandsBelowMinimum != thinThickness.islandsTested)
            throw std::runtime_error("native island thickness audit or thin-pad negative control failed");
        std::cout << "route_diagnostics pass=" << route.passed << " segments=" << route.segmentsTested
            << " centers=" << route.centerSamples << " capsule_samples=" << route.capsuleFieldSamples
            << " blocked=" << route.blockedCapsuleSamples << " support_samples=" << route.supportProbeSamples
            << " unsupported=" << route.unsupportedCenterSamples << " grade_failures=" << route.gradeFailures
            << " minimum_clearance_q=" << route.minimumClearanceSubvoxels
            << " first_blocked_segment=" << route.firstBlockedSegment
            << " first_unsupported_segment=" << route.firstUnsupportedSegment
            << " first_grade_segment=" << route.firstGradeFailureSegment
            << " unsupported_center_voxels=" << route.firstUnsupportedCenterSubvoxels[0] / static_cast<double>(Q)
            << ',' << route.firstUnsupportedCenterSubvoxels[1] / static_cast<double>(Q)
            << ',' << route.firstUnsupportedCenterSubvoxels[2] / static_cast<double>(Q) << '\n';
        if (!route.passed)
        {
            const auto& routePoints = planned->plan.route;
            if (route.firstBlockedSegment != ~std::uint64_t{0})
                std::cout << "first_blocked_position_voxels="
                    << route.firstBlockedPositionSubvoxels[0] / static_cast<double>(Q) << ','
                    << route.firstBlockedPositionSubvoxels[1] / static_cast<double>(Q) << ','
                    << route.firstBlockedPositionSubvoxels[2] / static_cast<double>(Q) << '\n';
            const std::size_t debugBegin = route.firstUnsupportedSegment == ~std::uint64_t{0}
                ? 0 : static_cast<std::size_t>(route.firstUnsupportedSegment);
            for (std::size_t i = debugBegin > 1 ? debugBegin - 1 : 0;
                 i < std::min(routePoints.size(), debugBegin + 3); ++i)
            {
                const auto& point = routePoints[i].positionSubvoxels;
                std::cout << "route_point=" << i << ' '
                    << point[0] / static_cast<double>(Q) << ','
                    << point[1] / static_cast<double>(Q) << ','
                    << point[2] / static_cast<double>(Q) << '\n';
            }
            const double maxGrade = static_cast<double>(planned->recipe.routeConstraints.maximumGradeQ16) / 65536.0;
            for (std::size_t i = 0; i + 1 < routePoints.size(); ++i)
            {
                const auto& a = routePoints[i].positionSubvoxels;
                const auto& b = routePoints[i + 1].positionSubvoxels;
                const double horizontal = std::hypot(static_cast<double>(b[0] - a[0]),
                    static_cast<double>(b[1] - a[1]));
                const double grade = horizontal > 0.0 ? std::abs(static_cast<double>(b[2] - a[2])) / horizontal : 1e9;
                if (grade > maxGrade + 1e-7)
                    std::cout << "grade_segment=" << i << " ratio=" << grade << " limit=" << maxGrade
                        << " from=" << a[0] / static_cast<double>(Q) << ',' << a[1] / static_cast<double>(Q)
                        << ',' << a[2] / static_cast<double>(Q) << " to=" << b[0] / static_cast<double>(Q)
                        << ',' << b[1] / static_cast<double>(Q) << ',' << b[2] / static_cast<double>(Q) << '\n';
            }
            for (std::size_t i = 0; i < planned->plan.islands.size(); ++i)
            {
                const auto& island = planned->plan.islands[i];
                std::cout << "island=" << i << " center=" << island.centerSubvoxels[0] / static_cast<double>(Q)
                    << ',' << island.centerSubvoxels[1] / static_cast<double>(Q) << ','
                    << island.topZSubvoxels / static_cast<double>(Q) << " radius="
                    << island.radiusSubvoxels / static_cast<double>(Q) << '\n';
            }
            const auto& constraints = planned->recipe.routeConstraints;
            const auto& unsupported = route.firstUnsupportedCenterSubvoxels;
            const double foot = unsupported[2] / static_cast<double>(Q)
                - constraints.capsuleHeightSubvoxels / (2.0 * Q);
            std::cout << "unsupported_route_probe=" << unsupported[0] / static_cast<double>(Q)
                << ',' << unsupported[1] / static_cast<double>(Q) << ','
                << unsupported[2] / static_cast<double>(Q) << " foot=" << foot << '\n';
            const double probeRadius = constraints.supportProbeRadiusSubvoxels / static_cast<double>(Q);
            const std::array<std::array<double, 2>, 5> supportPoints{{
                {0.0,0.0},{probeRadius,0.0},{-probeRadius,0.0},{0.0,probeRadius},{0.0,-probeRadius}}};
            for (std::size_t probeIndex = 0; probeIndex < supportPoints.size(); ++probeIndex)
            {
                double firstSolidDepth = -1.0;
                for (double depth = (constraints.minimumClearanceSubvoxels + 0.5 * Q) / static_cast<double>(Q);
                    depth <= (constraints.minimumClearanceSubvoxels + constraints.supportSearchDepthSubvoxels)
                        / static_cast<double>(Q) + 0.001; depth += 0.5)
                {
                    const auto density = field.densitySubvoxels(
                        std::llround(unsupported[0] + supportPoints[probeIndex][0] * Q),
                        std::llround(unsupported[1] + supportPoints[probeIndex][1] * Q),
                        std::llround((foot - depth) * Q));
                    if (density < 0) { firstSolidDepth = depth; break; }
                }
                std::cout << "support_probe=" << probeIndex << " first_solid_depth_vox=" << firstSolidDepth << '\n';
                std::cout << "support_profile=" << probeIndex << ':';
                for (double depth = 0.5; depth <= 5.0; depth += 0.5)
                {
                    const auto density = field.densitySubvoxels(
                        std::llround(unsupported[0] + supportPoints[probeIndex][0] * Q),
                        std::llround(unsupported[1] + supportPoints[probeIndex][1] * Q),
                        std::llround((foot - depth) * Q));
                    std::cout << ' ' << depth << '=' << density;
                }
                std::cout << '\n';
            }
            for (std::size_t i = 0; i < planned->plan.clusterEdges.size(); ++i)
            {
                const auto& edge = planned->plan.clusterEdges[i];
                std::cout << "connector=" << i << " start=" << edge.startSubvoxels[0] / static_cast<double>(Q)
                    << ',' << edge.startSubvoxels[1] / static_cast<double>(Q) << ','
                    << edge.startSubvoxels[2] / static_cast<double>(Q) << " end="
                    << edge.endSubvoxels[0] / static_cast<double>(Q) << ','
                    << edge.endSubvoxels[1] / static_cast<double>(Q) << ','
                    << edge.endSubvoxels[2] / static_cast<double>(Q) << '\n';
            }
        }
        if (!route.passed) throw std::runtime_error("native route capsule sweep failed");
        if (raisedRoute.passed || (raisedRoute.unsupportedCenterSamples == 0
            && raisedRoute.blockedCapsuleSamples == 0))
            throw std::runtime_error("raised-route negative control did not detect unsupported traversal");

        std::array<MeshResult, 4> meshes{};
        const auto clusterBounds = field.influenceBounds();
        if (!clusterBounds.finite) throw std::runtime_error("native cluster mesh bounds are not finite");
        const Point clusterMin{clusterBounds.minimumSubvoxels[0] - 8 * Q,
            clusterBounds.minimumSubvoxels[1] - 8 * Q, clusterBounds.minimumSubvoxels[2] - 8 * Q};
        const Point clusterMax{clusterBounds.maximumSubvoxels[0] + 8 * Q,
            clusterBounds.maximumSubvoxels[1] + 8 * Q, clusterBounds.maximumSubvoxels[2] + 8 * Q};
        const Point detailMin{coreAnchorQ[0] - 64 * Q, coreAnchorQ[1] - 64 * Q, coreAnchorQ[2] - 80 * Q};
        const Point detailMax{coreAnchorQ[0] + 64 * Q, coreAnchorQ[1] + 64 * Q, coreAnchorQ[2] + 16 * Q};
        meshes[0] = ExtractMesh(outputDir / "cluster_legacy.vfmesh", clusterMin, clusterMax, 4 * Q,
            [&](const Point& q) { return SampleLegacy(owner, q); });
        meshes[1] = ExtractMesh(outputDir / "cluster_native.vfmesh", clusterMin, clusterMax, 4 * Q,
            [&](const Point& q) { return SampleNative(field, q); });
        meshes[2] = ExtractMesh(outputDir / "detail_legacy.vfmesh", detailMin, detailMax, 2 * Q,
            [&](const Point& q) { return SampleLegacy(owner, q); });
        meshes[3] = ExtractMesh(outputDir / "detail_native.vfmesh", detailMin, detailMax, 2 * Q,
            [&](const Point& q) { return SampleNative(field, q); });
        if (meshes[0].triangles == 0 || meshes[1].triangles == 0
            || meshes[2].triangles == 0 || meshes[3].triangles == 0)
            throw std::runtime_error("matched legacy/native mesh extraction produced an empty field");

        WriteJson(outputDir / "audit.json", *planned, bounds, route, raisedRoute,
            w16, thickness, thinThickness, anchor, meshes, hashes, oldHashes);
        std::cout << std::setprecision(4)
            << "legacy_anchor_owner_voxels=" << anchor.surfaceQ[0] / static_cast<double>(Q) << ','
            << anchor.surfaceQ[1] / static_cast<double>(Q) << ','
            << anchor.surfaceQ[2] / static_cast<double>(Q)
            << " comparison_anchor_core_voxels=" << coreAnchorQ[0] / static_cast<double>(Q) << ','
            << coreAnchorQ[1] / static_cast<double>(Q) << ','
            << coreAnchorQ[2] / static_cast<double>(Q)
            << " plateau_samples=" << anchor.supportScore << '\n'
            << "recipe_hash=" << Hex64(planned->recipe.contentHash)
            << " plan_hash=" << Hex64(planned->plan.contentHash)
            << " manifest_hash=" << Hex64(planned->contentHash)
            << " prior_m1_hash=" << oldHashes.round1Manifest
            << " prior_surfaceworld_hash=" << oldHashes.surfaceWorldRound4Manifest
            << " islands=" << planned->plan.islands.size()
            << " connectors=" << planned->plan.clusterEdges.size()
            << " route_points=" << planned->plan.route.size()
            << " planning_ms=" << std::chrono::duration<double, std::milli>(endPlanning - beginPlanning).count() << '\n'
            << "bounds_pass=" << bounds.passed << " boxes=" << bounds.boxesTested
            << " point_samples=" << bounds.pointSamples << " unknown=" << bounds.unknownBounds
            << " missed=" << bounds.missedBounds << " void_air=" << bounds.voidBoxesProvedAir
            << " negative_control=" << bounds.negativeControlDetected << '\n'
            << "route_pass=" << route.passed << " segments=" << route.segmentsTested
            << " centers=" << route.centerSamples << " capsule_samples=" << route.capsuleFieldSamples
            << " blocked=" << route.blockedCapsuleSamples << " support_probes=" << route.supportProbeSamples
            << " unsupported=" << route.unsupportedCenterSamples << " grade_failures=" << route.gradeFailures
            << " clearance_q=" << route.minimumClearanceSubvoxels << '\n'
            << "raised_route_detected=" << (!raisedRoute.passed && raisedRoute.unsupportedCenterSamples > 0)
            << " unsupported=" << raisedRoute.unsupportedCenterSamples
            << " blocked=" << raisedRoute.blockedCapsuleSamples << '\n'
            << "thickness_pass=" << thickness.passed << " islands=" << thickness.islandsTested
            << " min_depth_diameter=" << thickness.minimumDepthToDiameterRatio
            << " threshold=" << thickness.requiredDepthToDiameterRatio
            << " thin_pad_detected=" << (!thinThickness.passed
                && thinThickness.islandsBelowMinimum == thinThickness.islandsTested)
            << " thin_ratio=" << thinThickness.minimumDepthToDiameterRatio << '\n'
            << "field_hashes W1=" << w1.hash << " W8=" << w8.hash << " W16=" << w16.hash
            << " legacy_ms=" << w16.legacyMilliseconds << " native_ms=" << w16.nativeMilliseconds
            << " mismatch_percent=" << (100.0 * w16.signMismatches / 32768.0)
            << " legacy_solid=" << w16.legacySolidSamples << " native_solid=" << w16.nativeSolidSamples
            << " legacy_only=" << w16.legacyOnlySolidSamples << " native_only=" << w16.nativeOnlySolidSamples << '\n'
            << "mesh_samples=" << meshes[0].gridSamples << ',' << meshes[1].gridSamples
            << ',' << meshes[2].gridSamples << ',' << meshes[3].gridSamples
            << " triangles=" << meshes[0].triangles << ',' << meshes[1].triangles
            << ',' << meshes[2].triangles << ',' << meshes[3].triangles << '\n'
            << "audit_json=" << (outputDir / "audit.json").string() << '\n';
        return 0;
    }
    catch (const std::exception& exception)
    {
        std::cerr << "FLOATING_ISLAND_AUDIT_FAIL: " << exception.what() << '\n';
        return 1;
    }
}