#pragma once

#include "MarchingCubesTables.h"
#include "VoxelMarchingCubesMesher.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <functional>
#include <limits>
#include <numeric>
#include <tuple>
#include <utility>
#include <vector>

// Feature-aware dual-contouring support for the opt-in Core recipe.
// Support du dual contouring sensible aux fonctions pour la recette Core optionnelle.
namespace VoxelHybridSharpMesher
{
    constexpr int32 EdgeCount = 12;
    constexpr int32 CornerCount = 8;
    constexpr double SubvoxelsPerVoxel = 256.0;

    constexpr int32 CornerOffsets[CornerCount][3] = {
        {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 int32 EdgeCorners[EdgeCount][2] = {
        {0,1}, {1,2}, {2,3}, {3,0}, {4,5}, {5,6},
        {6,7}, {7,4}, {0,4}, {1,5}, {2,6}, {3,7}
    };

    inline double AxisValue(const FVector& V, int32 Axis)
    {
        return Axis == 0 ? V.X : (Axis == 1 ? V.Y : V.Z);
    }

    inline void SetAxisValue(FVector& V, int32 Axis, double Value)
    {
        if (Axis == 0) V.X = Value;
        else if (Axis == 1) V.Y = Value;
        else V.Z = Value;
    }

    inline int32 IntAxisValue(const FIntVector& V, int32 Axis)
    {
        return Axis == 0 ? V.X : (Axis == 1 ? V.Y : V.Z);
    }

    struct FPlane
    {
        uint64 Id = 0;
        FVector Normal = FVector::UpVector;
    };

    struct FRoot
    {
        bool bActive = false;
        FVector Position = FVector::ZeroVector;
        FVector Gradient = FVector::UpVector;
        std::vector<FPlane> Planes;
        uint64 FeatureId = 0;
    };

    struct FTriangle
    {
        FVector Positions[3]{};
        FVector Normals[3]{};
        uint64 SplitKeys[3]{};
    };

    struct FResult
    {
        std::vector<uint8> CoreCellDCEdges;
        std::vector<FTriangle> Triangles;
        uint32 SharpCells = 0;
        uint32 DualCells = 0;
        uint32 ActiveEdgeRoots = 0;
        uint32 TaggedEdgeRoots = 0;
        double MinimumFeatureRootAbsSdf = std::numeric_limits<double>::infinity();
        uint32 FeatureRootSamples = 0;
        uint32 NearFeatureRootSamples = 0;
        uint32 BorderFallbackPatches = 0;

        bool UsesDC(int32 X, int32 Y, int32 Z, int32 Cells) const
        {
            if (X < 0 || Y < 0 || Z < 0 || X >= Cells || Y >= Cells || Z >= Cells)
                return false;
            const std::size_t Cell = (static_cast<std::size_t>(Z) * Cells + Y) * Cells + X;
            const std::size_t Offset = Cell * EdgeCount;
            for (int32 Edge = 0; Edge < EdgeCount && Offset + Edge < CoreCellDCEdges.size(); ++Edge)
                if (CoreCellDCEdges[Offset + Edge] != 0) return true;
            return false;
        }

        bool UsesDCEdge(int32 X, int32 Y, int32 Z, int32 Edge, int32 Cells) const
        {
            if (X < 0 || Y < 0 || Z < 0 || X >= Cells || Y >= Cells || Z >= Cells)
                return false;
            const std::size_t Cell = (static_cast<std::size_t>(Z) * Cells + Y) * Cells + X;
            const std::size_t Offset = Cell * EdgeCount + static_cast<std::size_t>(Edge);
            return Offset < CoreCellDCEdges.size() && (CoreCellDCEdges[Offset] & 1u) != 0;
        }
    };

    struct FCell
    {
        uint8 Case = 0;
        uint8 ComponentCount = 0;
        std::array<uint8, EdgeCount> EdgeComponent{};
        std::array<uint8, EdgeCount> ComponentSharp{};
        std::array<uint8, EdgeCount> ComponentUseDC{};
        std::array<FVector, EdgeCount> Dual{};

        FCell() { EdgeComponent.fill(0xff); }
    };

    inline int32 FindLocal(std::array<uint8, EdgeCount>& Parent, int32 Value)
    {
        while (Parent[Value] != Value)
        {
            Parent[Value] = Parent[Parent[Value]];
            Value = Parent[Value];
        }
        return Value;
    }

    inline void UniteLocal(std::array<uint8, EdgeCount>& Parent, int32 A, int32 B)
    {
        A = FindLocal(Parent, A);
        B = FindLocal(Parent, B);
        if (A == B) return;
        if (A > B) std::swap(A, B);
        Parent[B] = static_cast<uint8>(A);
    }

    inline int32 FindLocalEdge(int32 CellX, int32 CellY, int32 CellZ,
                               int32 Axis, int32 EdgeX, int32 EdgeY, int32 EdgeZ)
    {
        for (int32 Edge = 0; Edge < EdgeCount; ++Edge)
        {
            const int32 A = EdgeCorners[Edge][0];
            const int32 B = EdgeCorners[Edge][1];
            const int32 AX = CornerOffsets[A][0], AY = CornerOffsets[A][1], AZ = CornerOffsets[A][2];
            const int32 BX = CornerOffsets[B][0], BY = CornerOffsets[B][1], BZ = CornerOffsets[B][2];
            const int32 LocalAxis = AX != BX ? 0 : AY != BY ? 1 : 2;
            if (LocalAxis == Axis
                && CellX + FMath::Min(AX, BX) == EdgeX
                && CellY + FMath::Min(AY, BY) == EdgeY
                && CellZ + FMath::Min(AZ, BZ) == EdgeZ)
                return Edge;
        }
        return INDEX_NONE;
    }

    inline FVector SolveQEF(const std::vector<std::pair<FVector, FVector>>& Constraints,
                            const FVector& Fallback, const FVector& Low, const FVector& High)
    {
        if (Constraints.empty()) return Fallback;
        double A[3][4]{};
        constexpr double Lambda = 1.0e-10;
        const double Center[3] = {Fallback.X, Fallback.Y, Fallback.Z};
        for (int32 I = 0; I < 3; ++I) { A[I][I] = Lambda; A[I][3] = Lambda * Center[I]; }
        for (const auto& Constraint : Constraints)
        {
            FVector N = Constraint.second.GetSafeNormal();
            const double B = FVector::DotProduct(N, Constraint.first);
            const double NV[3] = {N.X, N.Y, N.Z};
            for (int32 R = 0; R < 3; ++R)
            {
                for (int32 C = 0; C < 3; ++C) A[R][C] += NV[R] * NV[C];
                A[R][3] += NV[R] * B;
            }
        }
        for (int32 Column = 0; Column < 3; ++Column)
        {
            int32 Pivot = Column;
            for (int32 Row = Column + 1; Row < 3; ++Row)
                if (FMath::Abs(A[Row][Column]) > FMath::Abs(A[Pivot][Column])) Pivot = Row;
            if (FMath::Abs(A[Pivot][Column]) < 1.0e-14) return Fallback;
            if (Pivot != Column)
                for (int32 K = Column; K < 4; ++K) std::swap(A[Pivot][K], A[Column][K]);
            const double Divisor = A[Column][Column];
            for (int32 K = Column; K < 4; ++K) A[Column][K] /= Divisor;
            for (int32 Row = 0; Row < 3; ++Row)
            {
                if (Row == Column) continue;
                const double Factor = A[Row][Column];
                for (int32 K = Column; K < 4; ++K) A[Row][K] -= Factor * A[Column][K];
            }
        }
        return FVector(FMath::Clamp(A[0][3], Low.X, High.X),
                       FMath::Clamp(A[1][3], Low.Y, High.Y),
                       FMath::Clamp(A[2][3], Low.Z, High.Z));
    }

    inline double BoxSDF(const FVoxelHybridSharpFeatureBox& Box, const FVector& P,
                         std::vector<FPlane>& OutPlanes)
    {
        FVector Low, High;
        for (int32 Axis = 0; Axis < 3; ++Axis)
        {
            SetAxisValue(Low, Axis,
                static_cast<double>(Box.MinimumSubvoxels[Axis]) / SubvoxelsPerVoxel);
            SetAxisValue(High, Axis,
                static_cast<double>(Box.MaximumSubvoxels[Axis]) / SubvoxelsPerVoxel);
        }
        const FVector Center = (Low + High) * 0.5;
        const FVector Half = (High - Low) * 0.5;
        const FVector Delta = P - Center;
        const FVector Q(FMath::Abs(Delta.X) - Half.X,
                        FMath::Abs(Delta.Y) - Half.Y,
                        FMath::Abs(Delta.Z) - Half.Z);
        const double MaxQ = FMath::Max(Q.X, FMath::Max(Q.Y, Q.Z));
        const FVector Outside(FMath::Max(Q.X, 0.0), FMath::Max(Q.Y, 0.0), FMath::Max(Q.Z, 0.0));
        const double OutsideLength = Outside.Size();
        const double SDF = OutsideLength + FMath::Min(MaxQ, 0.0);
        const double Tolerance = 1.0e-4;
        if (FMath::Abs(MaxQ) <= Tolerance)
        {
            for (int32 Axis = 0; Axis < 3; ++Axis)
            {
                if (FMath::Abs(AxisValue(Q, Axis) - MaxQ) > Tolerance) continue;
                FVector Normal = FVector::ZeroVector;
                const bool bNegativeFace = AxisValue(Delta, Axis) < 0.0;
                SetAxisValue(Normal, Axis, bNegativeFace ? -1.0 : 1.0);
                uint64 FaceId = Box.StableId ^ (0x9e3779b97f4a7c15ull
                    + static_cast<uint64>(Axis * 2 + (bNegativeFace ? 0 : 1) + 1));
                if (FaceId == 0) FaceId = static_cast<uint64>(Axis + 1);
                OutPlanes.push_back({FaceId, Normal});
            }
        }
        return SDF;
    }

    struct FPreparedFeatureBox
    {
        uint64 StableId = 0;
        double Minimum[3]{};
        double Maximum[3]{};
    };

    inline std::vector<FPreparedFeatureBox> PrepareFeatureBoxes(
        const TArray<FVoxelHybridSharpFeatureBox>& Features)
    {
        std::vector<FPreparedFeatureBox> Prepared;
        Prepared.reserve(static_cast<std::size_t>(Features.Num()));
        for (const FVoxelHybridSharpFeatureBox& Feature : Features)
        {
            FPreparedFeatureBox Box;
            Box.StableId = Feature.StableId;
            for (int32 Axis = 0; Axis < 3; ++Axis)
            {
                Box.Minimum[Axis] = static_cast<double>(Feature.MinimumSubvoxels[Axis])
                    / SubvoxelsPerVoxel;
                Box.Maximum[Axis] = static_cast<double>(Feature.MaximumSubvoxels[Axis])
                    / SubvoxelsPerVoxel;
            }
            Prepared.push_back(Box);
        }
        return Prepared;
    }

    // MC interpolates the sampled field. A finite hard box is only included by the field at
    // samples inside it, so an edge can interpolate away from the box plane when the box face
    // lies between lattice points. For V2, keep the MC connectivity and move only roots whose
    // lattice edge crosses a tagged analytic face. This is deterministic and leaves tile edges
    // to the unchanged MC/skirt path.
    inline bool TrySnapFeatureEdge(const std::vector<FPreparedFeatureBox>& Features,
        const FIntVector& OriginVoxels, const FIntVector& LocalStart, int32 Axis,
        const FVector& ChunkWorldPos, float GeometryVoxelSize,
        FVector& InOutWorldPosition, FVector& OutNormal, uint64& OutSplitKey)
    {
        const double EdgeStart = static_cast<double>(IntAxisValue(OriginVoxels, Axis)
            + IntAxisValue(LocalStart, Axis));
        const double EdgeEnd = EdgeStart + 1.0;
        const double EdgeLow = FMath::Min(EdgeStart, EdgeEnd);
        const double EdgeHigh = FMath::Max(EdgeStart, EdgeEnd);
        const double RootVoxel = static_cast<double>(IntAxisValue(OriginVoxels, Axis))
            + (AxisValue(InOutWorldPosition, Axis) - AxisValue(ChunkWorldPos, Axis))
                / GeometryVoxelSize;
        constexpr double Tolerance = 1.0e-9;
        double BestDistance = std::numeric_limits<double>::infinity();
        uint64 BestStableId = std::numeric_limits<uint64>::max();
        uint64 BestFaceId = 0;
        double BestPlane = 0.0;
        int32 BestSide = 0;

        for (const FPreparedFeatureBox& Feature : Features)
        {
            bool bWithinFaceSpan = true;
            for (int32 OtherAxis = 0; OtherAxis < 3; ++OtherAxis)
            {
                if (OtherAxis == Axis) continue;
                const double Coordinate = static_cast<double>(IntAxisValue(OriginVoxels, OtherAxis)
                    + IntAxisValue(LocalStart, OtherAxis));
                if (Coordinate < Feature.Minimum[OtherAxis] - Tolerance
                    || Coordinate > Feature.Maximum[OtherAxis] + Tolerance)
                {
                    bWithinFaceSpan = false;
                    break;
                }
            }
            if (!bWithinFaceSpan) continue;

            for (int32 Side = 0; Side < 2; ++Side)
            {
                const double Plane = Side == 0
                    ? Feature.Minimum[Axis] : Feature.Maximum[Axis];
                if (Plane < EdgeLow - Tolerance || Plane > EdgeHigh + Tolerance) continue;
                const double Distance = FMath::Abs(Plane - RootVoxel);
                const uint64 FaceId = Feature.StableId ^ (0x9e3779b97f4a7c15ull
                    + static_cast<uint64>(Axis * 2 + Side + 1));
                if (Distance > BestDistance + Tolerance
                    || (FMath::Abs(Distance - BestDistance) <= Tolerance
                        && std::tie(Feature.StableId, FaceId)
                            >= std::tie(BestStableId, BestFaceId)))
                {
                    continue;
                }
                BestDistance = Distance;
                BestStableId = Feature.StableId;
                BestFaceId = FaceId != 0 ? FaceId : static_cast<uint64>(Axis + 1);
                BestPlane = Plane;
                BestSide = Side;
            }
        }

        if (BestFaceId == 0) return false;
        SetAxisValue(InOutWorldPosition, Axis, AxisValue(ChunkWorldPos, Axis)
            + (BestPlane - static_cast<double>(IntAxisValue(OriginVoxels, Axis)))
                * GeometryVoxelSize);
        OutNormal = FVector::ZeroVector;
        SetAxisValue(OutNormal, Axis, BestSide == 0 ? -1.0 : 1.0);
        OutSplitKey = BestFaceId;
        return true;
    }

    inline uint64 SelectFeature(const FRoot& Root, int32 Axis)
    {
        if (Root.Planes.empty()) return 0;
        std::size_t Selected = 0;
        const auto Alignment = [Axis](const FVector& Normal)
        { return FMath::Abs(AxisValue(Normal, Axis)); };
        for (std::size_t I = 1; I < Root.Planes.size(); ++I)
        {
            const double Candidate = Alignment(Root.Planes[I].Normal);
            const double Current = Alignment(Root.Planes[Selected].Normal);
            if (Candidate > Current + 1.0e-12
                || (FMath::Abs(Candidate - Current) <= 1.0e-12
                    && Root.Planes[I].Id < Root.Planes[Selected].Id)) Selected = I;
        }
        return Root.Planes[Selected].Id;
    }

    inline FVector NormalForFeature(const FRoot& Root, uint64 Feature)
    {
        for (const FPlane& Plane : Root.Planes) if (Plane.Id == Feature) return Plane.Normal;
        return Root.Gradient;
    }

    inline FResult Generate(int32 Cells, int32 Step, float IsoLevel, float GeometryVoxelSize,
                            const FIntVector& OriginVoxels, const FVector& ChunkWorldPos,
                            const TArray<FVoxelHybridSharpFeatureBox>& Features,
                            const std::function<float(int32,int32,int32)>& Sample,
                            const std::function<FVector(int32,int32,int32)>& Gradient)
    {
        FResult Result;
        if (Cells <= 0 || Features.IsEmpty()) return Result;
        const int32 CellDim = Cells + 2;
        const int32 NodeDim = Cells + 3;
        const std::size_t CellCount = static_cast<std::size_t>(CellDim) * CellDim * CellDim;
        std::vector<FCell> CellData(CellCount);
        const auto CellIndex = [CellDim](int32 X, int32 Y, int32 Z)
        { return (static_cast<std::size_t>(Z + 1) * CellDim + (Y + 1)) * CellDim + (X + 1); };
        const auto NodePosition = [&](int32 X, int32 Y, int32 Z)
        { return FVector(X * Step * GeometryVoxelSize, Y * Step * GeometryVoxelSize,
                         Z * Step * GeometryVoxelSize) + ChunkWorldPos; };

        std::array<std::vector<FRoot>, 3> Roots;
        Roots[0].resize(static_cast<std::size_t>(CellDim) * NodeDim * NodeDim);
        Roots[1].resize(static_cast<std::size_t>(NodeDim) * CellDim * NodeDim);
        Roots[2].resize(static_cast<std::size_t>(NodeDim) * NodeDim * CellDim);
        const auto RootIndex = [CellDim, NodeDim](int32 Axis, int32 X, int32 Y, int32 Z)
        {
            const std::size_t NX = static_cast<std::size_t>(X + 1);
            const std::size_t NY = static_cast<std::size_t>(Y + 1);
            const std::size_t NZ = static_cast<std::size_t>(Z + 1);
            if (Axis == 0) return (NZ * NodeDim + NY) * CellDim + NX;
            if (Axis == 1) return (NZ * CellDim + NY) * NodeDim + NX;
            return (NZ * NodeDim + NY) * NodeDim + NX;
        };

        for (int32 Z = -1; Z <= Cells + 1; ++Z)
        for (int32 Y = -1; Y <= Cells + 1; ++Y)
        for (int32 X = -1; X <= Cells; ++X)
        {
            const float D0 = Sample(X, Y, Z), D1 = Sample(X + 1, Y, Z);
            if ((D0 < IsoLevel) == (D1 < IsoLevel)) continue;
            const double T = FMath::Clamp(static_cast<double>(IsoLevel - D0) / (D1 - D0), 0.0, 1.0);
            FRoot& Root = Roots[0][RootIndex(0, X, Y, Z)];
            Root.bActive = true;
            ++Result.ActiveEdgeRoots;
            Root.Position = FMath::Lerp(NodePosition(X, Y, Z), NodePosition(X + 1, Y, Z), T);
            Root.Gradient = FMath::Lerp(Gradient(X, Y, Z), Gradient(X + 1, Y, Z), T).GetSafeNormal();
        }
        for (int32 Z = -1; Z <= Cells + 1; ++Z)
        for (int32 Y = -1; Y <= Cells; ++Y)
        for (int32 X = -1; X <= Cells + 1; ++X)
        {
            const float D0 = Sample(X, Y, Z), D1 = Sample(X, Y + 1, Z);
            if ((D0 < IsoLevel) == (D1 < IsoLevel)) continue;
            const double T = FMath::Clamp(static_cast<double>(IsoLevel - D0) / (D1 - D0), 0.0, 1.0);
            FRoot& Root = Roots[1][RootIndex(1, X, Y, Z)];
            Root.bActive = true;
            ++Result.ActiveEdgeRoots;
            Root.Position = FMath::Lerp(NodePosition(X, Y, Z), NodePosition(X, Y + 1, Z), T);
            Root.Gradient = FMath::Lerp(Gradient(X, Y, Z), Gradient(X, Y + 1, Z), T).GetSafeNormal();
        }
        for (int32 Z = -1; Z <= Cells; ++Z)
        for (int32 Y = -1; Y <= Cells + 1; ++Y)
        for (int32 X = -1; X <= Cells + 1; ++X)
        {
            const float D0 = Sample(X, Y, Z), D1 = Sample(X, Y, Z + 1);
            if ((D0 < IsoLevel) == (D1 < IsoLevel)) continue;
            const double T = FMath::Clamp(static_cast<double>(IsoLevel - D0) / (D1 - D0), 0.0, 1.0);
            FRoot& Root = Roots[2][RootIndex(2, X, Y, Z)];
            Root.bActive = true;
            ++Result.ActiveEdgeRoots;
            Root.Position = FMath::Lerp(NodePosition(X, Y, Z), NodePosition(X, Y, Z + 1), T);
            Root.Gradient = FMath::Lerp(Gradient(X, Y, Z), Gradient(X, Y, Z + 1), T).GetSafeNormal();
        }

        const auto GetRoot = [&](int32 Axis, int32 X, int32 Y, int32 Z) -> FRoot&
        { return Roots[Axis][RootIndex(Axis, X, Y, Z)]; };
        const auto GetCellEdgeRoot = [&](int32 X, int32 Y, int32 Z, int32 Edge) -> FRoot&
        {
            const int32 A = EdgeCorners[Edge][0], B = EdgeCorners[Edge][1];
            const int32 AX = CornerOffsets[A][0], AY = CornerOffsets[A][1], AZ = CornerOffsets[A][2];
            const int32 BX = CornerOffsets[B][0], BY = CornerOffsets[B][1], BZ = CornerOffsets[B][2];
            const int32 Axis = AX != BX ? 0 : AY != BY ? 1 : 2;
            return GetRoot(Axis, X + FMath::Min(AX, BX), Y + FMath::Min(AY, BY), Z + FMath::Min(AZ, BZ));
        };

        for (int32 Z = -1; Z <= Cells; ++Z)
        for (int32 Y = -1; Y <= Cells; ++Y)
        for (int32 X = -1; X <= Cells; ++X)
        {
            FCell& Cell = CellData[CellIndex(X, Y, Z)];
            float Densities[CornerCount]{};
            for (int32 Corner = 0; Corner < CornerCount; ++Corner)
            {
                const int32 GX = X + CornerOffsets[Corner][0];
                const int32 GY = Y + CornerOffsets[Corner][1];
                const int32 GZ = Z + CornerOffsets[Corner][2];
                Densities[Corner] = Sample(GX, GY, GZ);
                if (Densities[Corner] >= IsoLevel) Cell.Case |= static_cast<uint8>(1u << Corner);
            }
            if (Cell.Case == 0 || Cell.Case == 255 || MCEdgeTable[Cell.Case] == 0) continue;
            std::array<uint8, EdgeCount> Parent{};
            for (int32 Edge = 0; Edge < EdgeCount; ++Edge) Parent[Edge] = static_cast<uint8>(Edge);
            for (int32 I = 0; I < 16 && MCTriTable[Cell.Case][I] >= 0; I += 3)
            {
                const int32 A = MCTriTable[Cell.Case][I];
                const int32 B = MCTriTable[Cell.Case][I + 1];
                const int32 C = MCTriTable[Cell.Case][I + 2];
                UniteLocal(Parent, A, B); UniteLocal(Parent, B, C);
            }
            std::array<uint8, EdgeCount> RootToComponent{};
            RootToComponent.fill(0xff);
            for (int32 Edge = 0; Edge < EdgeCount; ++Edge)
            {
                const int32 A = EdgeCorners[Edge][0], B = EdgeCorners[Edge][1];
                if ((Densities[A] < IsoLevel) == (Densities[B] < IsoLevel)) continue;
                FRoot& Root = GetCellEdgeRoot(X, Y, Z, Edge);
                if (!Root.bActive) continue;
                const int32 LocalRoot = FindLocal(Parent, Edge);
                if (RootToComponent[LocalRoot] == 0xff)
                    RootToComponent[LocalRoot] = Cell.ComponentCount++;
                const uint8 Component = RootToComponent[LocalRoot];
                Cell.EdgeComponent[Edge] = Component;

                const int32 OA = EdgeCorners[Edge][0];
                const int32 OB = EdgeCorners[Edge][1];
                const int32 AX = X + CornerOffsets[OA][0], AY = Y + CornerOffsets[OA][1], AZ = Z + CornerOffsets[OA][2];
                const int32 BX = X + CornerOffsets[OB][0], BY = Y + CornerOffsets[OB][1], BZ = Z + CornerOffsets[OB][2];
                const FVector RootVoxel = FVector(OriginVoxels.X, OriginVoxels.Y, OriginVoxels.Z)
                    + FVector(AX, AY, AZ) * Step
                    + FVector(BX - AX, BY - AY, BZ - AZ) * (Step *
                        FMath::Clamp(static_cast<double>(IsoLevel - Densities[OA])
                            / (Densities[OB] - Densities[OA]), 0.0, 1.0));
                for (const FVoxelHybridSharpFeatureBox& Feature : Features)
                {
                    std::vector<FPlane> CandidatePlanes;
                    const double SDF = BoxSDF(Feature, RootVoxel, CandidatePlanes);
                    const double AbsSDF = FMath::Abs(SDF);
                    Result.MinimumFeatureRootAbsSdf = FMath::Min(
                        Result.MinimumFeatureRootAbsSdf, AbsSDF);
                    ++Result.FeatureRootSamples;
                    if (AbsSDF <= 1.0e-2) ++Result.NearFeatureRootSamples;
                    if (AbsSDF > 1.0e-4 || CandidatePlanes.empty()) continue;
                    for (const FPlane& Plane : CandidatePlanes) Root.Planes.push_back(Plane);
                }
                if (!Root.Planes.empty())
                {
                    ++Result.TaggedEdgeRoots;
                    std::sort(Root.Planes.begin(), Root.Planes.end(),
                        [](const FPlane& A, const FPlane& B) { return A.Id < B.Id; });
                    Root.Planes.erase(std::unique(Root.Planes.begin(), Root.Planes.end(),
                        [](const FPlane& A, const FPlane& B) { return A.Id == B.Id; }), Root.Planes.end());
                    Cell.ComponentSharp[Component] = 1;
                    if (X >= 0 && Y >= 0 && Z >= 0 && X < Cells && Y < Cells && Z < Cells)
                        ++Result.SharpCells;
                }
            }
        }

        const std::size_t PatchCount = CellCount * EdgeCount;
        std::vector<uint32> PatchParent(PatchCount);
        std::vector<uint8> PatchRank(PatchCount, 0), PatchSharp(PatchCount, 0), PatchTouchesBorder(PatchCount, 0);
        std::iota(PatchParent.begin(), PatchParent.end(), 0u);
        const auto FindPatch = [&](uint32 Value, auto&& Self) -> uint32
        {
            if (PatchParent[Value] != Value) PatchParent[Value] = Self(PatchParent[Value], Self);
            return PatchParent[Value];
        };
        const auto UnitePatch = [&](uint32 A, uint32 B)
        {
            uint32 RA = FindPatch(A, FindPatch), RB = FindPatch(B, FindPatch);
            if (RA == RB) return;
            if (PatchRank[RA] < PatchRank[RB]) std::swap(RA, RB);
            PatchParent[RB] = RA;
            PatchSharp[RA] = static_cast<uint8>(PatchSharp[RA] | PatchSharp[RB]);
            PatchTouchesBorder[RA] = static_cast<uint8>(PatchTouchesBorder[RA] | PatchTouchesBorder[RB]);
            if (PatchRank[RA] == PatchRank[RB]) ++PatchRank[RA];
        };
        for (int32 Z = -1; Z <= Cells; ++Z)
        for (int32 Y = -1; Y <= Cells; ++Y)
        for (int32 X = -1; X <= Cells; ++X)
        {
            const FCell& Cell = CellData[CellIndex(X, Y, Z)];
            for (uint8 Component = 0; Component < Cell.ComponentCount; ++Component)
            {
                const uint32 Id = static_cast<uint32>(CellIndex(X, Y, Z) * EdgeCount + Component);
                PatchSharp[Id] = Cell.ComponentSharp[Component];
                const bool bCoreBorder = X >= 0 && Y >= 0 && Z >= 0 && X < Cells && Y < Cells && Z < Cells
                    && (X == 0 || Y == 0 || Z == 0 || X == Cells - 1 || Y == Cells - 1 || Z == Cells - 1);
                PatchTouchesBorder[Id] = bCoreBorder ? 1 : 0;
            }
        }

        for (int32 Axis = 0; Axis < 3; ++Axis)
        for (int32 Z = -1; Z <= Cells + (Axis == 2 ? 0 : 1); ++Z)
        for (int32 Y = -1; Y <= Cells + (Axis == 1 ? 0 : 1); ++Y)
        for (int32 X = -1; X <= Cells + (Axis == 0 ? 0 : 1); ++X)
        {
            if (!GetRoot(Axis, X, Y, Z).bActive) continue;
            std::array<std::array<int32, 3>, 4> Around{};
            if (Axis == 0) Around = {{{X,Y-1,Z-1},{X,Y,Z-1},{X,Y,Z},{X,Y-1,Z}}};
            else if (Axis == 1) Around = {{{X-1,Y,Z-1},{X-1,Y,Z},{X,Y,Z},{X,Y,Z-1}}};
            else Around = {{{X-1,Y-1,Z},{X,Y-1,Z},{X,Y,Z},{X-1,Y,Z}}};
            std::array<uint32, 4> Ids{};
            bool bValid = true;
            for (int32 I = 0; I < 4; ++I)
            {
                const int32 CX = Around[I][0], CY = Around[I][1], CZ = Around[I][2];
                if (CX < -1 || CY < -1 || CZ < -1 || CX > Cells || CY > Cells || CZ > Cells)
                { bValid = false; break; }
                const int32 LocalEdge = FindLocalEdge(CX, CY, CZ, Axis, X, Y, Z);
                const FCell& Cell = CellData[CellIndex(CX, CY, CZ)];
                if (LocalEdge == INDEX_NONE || Cell.EdgeComponent[LocalEdge] == 0xff)
                { bValid = false; break; }
                Ids[I] = static_cast<uint32>(CellIndex(CX, CY, CZ) * EdgeCount + Cell.EdgeComponent[LocalEdge]);
            }
            if (bValid) for (int32 I = 1; I < 4; ++I) UnitePatch(Ids[0], Ids[I]);
        }

        bool bSawBorderPatch = false;
        std::vector<uint8> PatchUseDC(PatchCount, 0);
        for (int32 Z = -1; Z <= Cells; ++Z)
        for (int32 Y = -1; Y <= Cells; ++Y)
        for (int32 X = -1; X <= Cells; ++X)
        {
            const FCell& Cell = CellData[CellIndex(X, Y, Z)];
            for (uint8 Component = 0; Component < Cell.ComponentCount; ++Component)
            {
                const uint32 Id = static_cast<uint32>(CellIndex(X, Y, Z) * EdgeCount + Component);
                const uint32 Root = FindPatch(Id, FindPatch);
                PatchUseDC[Root] = static_cast<uint8>(PatchSharp[Root] && !PatchTouchesBorder[Root]);
                if (PatchSharp[Root] && PatchTouchesBorder[Root]) bSawBorderPatch = true;
            }
        }
        Result.BorderFallbackPatches = bSawBorderPatch ? 1u : 0u;

        for (int32 Z = -1; Z <= Cells; ++Z)
        for (int32 Y = -1; Y <= Cells; ++Y)
        for (int32 X = -1; X <= Cells; ++X)
        {
            FCell& Cell = CellData[CellIndex(X, Y, Z)];
            if (Cell.ComponentCount == 0) continue;
            for (uint8 Component = 0; Component < Cell.ComponentCount; ++Component)
            {
                const uint32 Id = static_cast<uint32>(CellIndex(X, Y, Z) * EdgeCount + Component);
                Cell.ComponentUseDC[Component] = PatchUseDC[FindPatch(Id, FindPatch)];
                if (!Cell.ComponentUseDC[Component]) continue;
                std::vector<std::pair<FVector, FVector>> Constraints;
                FVector RootSum = FVector::ZeroVector;
                int32 RootCount = 0;
                for (int32 Edge = 0; Edge < EdgeCount; ++Edge)
                {
                    if (Cell.EdgeComponent[Edge] != Component) continue;
                    FRoot& Root = GetCellEdgeRoot(X, Y, Z, Edge);
                    if (!Root.bActive) continue;
                    RootSum += Root.Position;
                    ++RootCount;
                    if (!Root.Planes.empty())
                        for (const FPlane& Plane : Root.Planes) Constraints.emplace_back(Root.Position, Plane.Normal);
                    else Constraints.emplace_back(Root.Position, Root.Gradient);
                }
                const FVector Low = NodePosition(X, Y, Z);
                const FVector High = NodePosition(X + 1, Y + 1, Z + 1);
                const FVector Center = RootCount > 0 ? RootSum / RootCount : (Low + High) * 0.5;
                Cell.Dual[Component] = SolveQEF(Constraints, Center, Low, High);
                if (X >= 0 && Y >= 0 && Z >= 0 && X < Cells && Y < Cells && Z < Cells)
                    ++Result.DualCells;
            }
        }

        Result.CoreCellDCEdges.assign(static_cast<std::size_t>(Cells) * Cells * Cells * EdgeCount, 0);
        const auto AddTriangle = [&](const FVector& A, const FVector& B, const FVector& C,
                                     const FVector& NA, const FVector& NB, const FVector& NC,
                                     uint64 FA, uint64 FB, uint64 FC)
        {
            if (FVector::CrossProduct(B - A, C - A).SizeSquared() < 1.0e-12) return;
            FTriangle T;
            T.Positions[0] = A; T.Positions[1] = B; T.Positions[2] = C;
            T.Normals[0] = NA.GetSafeNormal(); T.Normals[1] = NB.GetSafeNormal(); T.Normals[2] = NC.GetSafeNormal();
            T.SplitKeys[0] = FA; T.SplitKeys[1] = FB; T.SplitKeys[2] = FC;
            Result.Triangles.push_back(T);
        };

        for (int32 Axis = 0; Axis < 3; ++Axis)
        for (int32 Z = -1; Z <= Cells + (Axis == 2 ? 0 : 1); ++Z)
        for (int32 Y = -1; Y <= Cells + (Axis == 1 ? 0 : 1); ++Y)
        for (int32 X = -1; X <= Cells + (Axis == 0 ? 0 : 1); ++X)
        {
            FRoot& Root = GetRoot(Axis, X, Y, Z);
            if (!Root.bActive) continue;
            std::array<std::array<int32, 3>, 4> Around{};
            if (Axis == 0) Around = {{{X,Y-1,Z-1},{X,Y,Z-1},{X,Y,Z},{X,Y-1,Z}}};
            else if (Axis == 1) Around = {{{X-1,Y,Z-1},{X-1,Y,Z},{X,Y,Z},{X,Y,Z-1}}};
            else Around = {{{X-1,Y-1,Z},{X,Y-1,Z},{X,Y,Z},{X-1,Y,Z}}};
            std::array<const FCell*, 4> Incident{};
            std::array<uint8, 4> Components{};
            bool bAllDC = true;
            bool bValid = true;
            for (int32 I = 0; I < 4; ++I)
            {
                const int32 CX = Around[I][0], CY = Around[I][1], CZ = Around[I][2];
                if (CX < -1 || CY < -1 || CZ < -1 || CX > Cells || CY > Cells || CZ > Cells)
                { bValid = false; break; }
                const int32 LocalEdge = FindLocalEdge(CX, CY, CZ, Axis, X, Y, Z);
                const FCell& Cell = CellData[CellIndex(CX, CY, CZ)];
                if (LocalEdge == INDEX_NONE || Cell.EdgeComponent[LocalEdge] == 0xff)
                { bValid = false; break; }
                const uint8 Component = Cell.EdgeComponent[LocalEdge];
                Incident[I] = &Cell;
                Components[I] = Component;
                bAllDC &= Cell.ComponentUseDC[Component] != 0;
            }
            if (!bValid || !bAllDC) continue;
            const int32 OwnerX = FMath::Min(FMath::Min(Around[0][0], Around[1][0]), FMath::Min(Around[2][0], Around[3][0]));
            const int32 OwnerY = FMath::Min(FMath::Min(Around[0][1], Around[1][1]), FMath::Min(Around[2][1], Around[3][1]));
            const int32 OwnerZ = FMath::Min(FMath::Min(Around[0][2], Around[1][2]), FMath::Min(Around[2][2], Around[3][2]));
            if (OwnerX < 0 || OwnerY < 0 || OwnerZ < 0 || OwnerX >= Cells || OwnerY >= Cells || OwnerZ >= Cells)
                continue;

            std::array<FVector, 4> Positions{};
            std::array<FVector, 4> Normals{};
            std::array<uint64, 4> SplitKeys{};
            for (int32 I = 0; I < 4; ++I)
            {
                const int32 CX = Around[I][0], CY = Around[I][1], CZ = Around[I][2];
                const uint8 Component = Components[I];
                Positions[I] = Incident[I]->Dual[Component];
                SplitKeys[I] = SelectFeature(Root, Axis);
                Normals[I] = SplitKeys[I] != 0 ? NormalForFeature(Root, SplitKeys[I]) : Root.Gradient;
            }
            if (FVector::DotProduct(FVector::CrossProduct(Positions[1] - Positions[0],
                    Positions[2] - Positions[0]), Root.Gradient) < 0.0)
            {
                std::swap(Positions[1], Positions[3]);
                std::swap(Normals[1], Normals[3]);
                std::swap(SplitKeys[1], SplitKeys[3]);
            }
            const int32 Tie = (X + Y + Z + Axis) & 1;
            if (Tie == 0)
            {
                AddTriangle(Positions[0], Positions[1], Positions[2], Normals[0], Normals[1], Normals[2], SplitKeys[0], SplitKeys[1], SplitKeys[2]);
                AddTriangle(Positions[0], Positions[2], Positions[3], Normals[0], Normals[2], Normals[3], SplitKeys[0], SplitKeys[2], SplitKeys[3]);
            }
            else
            {
                AddTriangle(Positions[0], Positions[1], Positions[3], Normals[0], Normals[1], Normals[3], SplitKeys[0], SplitKeys[1], SplitKeys[3]);
                AddTriangle(Positions[1], Positions[2], Positions[3], Normals[1], Normals[2], Normals[3], SplitKeys[1], SplitKeys[2], SplitKeys[3]);
            }
        }

        for (int32 Z = 0; Z < Cells; ++Z)
        for (int32 Y = 0; Y < Cells; ++Y)
        for (int32 X = 0; X < Cells; ++X)
        {
            const FCell& Cell = CellData[CellIndex(X, Y, Z)];
            if (!Cell.ComponentCount) continue;
            const std::size_t CoreCell = (static_cast<std::size_t>(Z) * Cells + Y) * Cells + X;
            for (int32 Edge = 0; Edge < EdgeCount; ++Edge)
            {
                const uint8 Component = Cell.EdgeComponent[Edge];
                if (Component != 0xff && Cell.ComponentUseDC[Component])
                    Result.CoreCellDCEdges[CoreCell * EdgeCount + Edge] = 1;
            }
        }
        return Result;
    }
}