Fr0zkaandClaude Opus 5.5
4661a445be
The round-3 'field contact error' measured hit-to-anchor distance; V2 contacts are 0 cm and on the zero isosurface. Tagged-tile cost 5.7-6.7x -> 1.175-1.270x MC; organic 0.998x. Bound, shelf-thickness and route audits rerun with negative controls; W1/W8/W16 stable; V1 and MarchingCubesV1 unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MmSoF64BbRFmk9C3p1VsBJ
727 lines
34 KiB
C++
727 lines
34 KiB
C++
#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;
|
|
}
|
|
}
|