Fr0zkaandClaude Opus 5.5
6953254349
Standalone only. Hybrid/DC cut sharp residual 16.7 -> 6.25 cm and split creases ~100%, organic wall unchanged under hybrid, +4.4% mesh cost, W1/W8/W16 identical. Gate failed: hybrid open edges on 500/625 tagged tile borders; every method fails LOD0/LOD1 (no transition templates yet). Binaries and OBJ meshes not committed. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MmSoF64BbRFmk9C3p1VsBJ
1175 lines
55 KiB
C++
1175 lines
55 KiB
C++
#include "CoreMinimal.h"
|
|
#include "MarchingCubesTables.h"
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <atomic>
|
|
#include <chrono>
|
|
#include <cmath>
|
|
#include <cstdint>
|
|
#include <cstdio>
|
|
#include <filesystem>
|
|
#include <fstream>
|
|
#include <functional>
|
|
#include <iomanip>
|
|
#include <iostream>
|
|
#include <map>
|
|
#include <numeric>
|
|
#include <string>
|
|
#include <thread>
|
|
#include <tuple>
|
|
#include <unordered_map>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
namespace
|
|
{
|
|
constexpr double VoxelMeters = 0.25;
|
|
constexpr double PositionQuantum = VoxelMeters / 65536.0;
|
|
constexpr double NormalSplitCos = 0.86602540378443864676; // 30 degrees
|
|
constexpr int CornerCount = 8;
|
|
constexpr int EdgeCount = 12;
|
|
constexpr std::array<std::array<int, 3>, CornerCount> CornerOffsets{{
|
|
{{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, 2>, EdgeCount> EdgeCorners{{
|
|
{{0,1}}, {{1,2}}, {{2,3}}, {{3,0}}, {{4,5}}, {{5,6}},
|
|
{{6,7}}, {{7,4}}, {{0,4}}, {{1,5}}, {{2,6}}, {{3,7}}
|
|
}};
|
|
|
|
struct V3
|
|
{
|
|
double x=0, y=0, z=0;
|
|
V3 operator+(V3 b) const { return {x+b.x,y+b.y,z+b.z}; }
|
|
V3 operator-(V3 b) const { return {x-b.x,y-b.y,z-b.z}; }
|
|
V3 operator*(double s) const { return {x*s,y*s,z*s}; }
|
|
V3 operator/(double s) const { return {x/s,y/s,z/s}; }
|
|
V3& operator+=(V3 b) { x+=b.x; y+=b.y; z+=b.z; return *this; }
|
|
};
|
|
double Dot(V3 a,V3 b){return a.x*b.x+a.y*b.y+a.z*b.z;}
|
|
V3 Cross(V3 a,V3 b){return {a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};}
|
|
double Length(V3 a){return std::sqrt(Dot(a,a));}
|
|
V3 Normalize(V3 a)
|
|
{
|
|
const double n=Length(a);
|
|
return n>1.0e-20?a/n:V3{0,0,1};
|
|
}
|
|
V3 Quantize(V3 p)
|
|
{
|
|
const auto q=[](double x){return std::nearbyint(x/PositionQuantum)*PositionQuantum;};
|
|
return {q(p.x),q(p.y),q(p.z)};
|
|
}
|
|
|
|
template<class F> void ParallelFor(std::size_t count,int workers,F&& function)
|
|
{
|
|
workers=std::clamp(workers,1,16);
|
|
if(workers==1||count<512)
|
|
{
|
|
for(std::size_t i=0;i<count;++i)function(i);
|
|
return;
|
|
}
|
|
const int actual=std::min<int>(workers,static_cast<int>(count));
|
|
std::vector<std::thread> threads;
|
|
threads.reserve(actual);
|
|
for(int worker=0;worker<actual;++worker)
|
|
{
|
|
const std::size_t begin=count*static_cast<std::size_t>(worker)/actual;
|
|
const std::size_t end=count*static_cast<std::size_t>(worker+1)/actual;
|
|
threads.emplace_back([begin,end,&function](){for(std::size_t i=begin;i<end;++i)function(i);});
|
|
}
|
|
for(auto& thread:threads)thread.join();
|
|
}
|
|
|
|
struct Plane
|
|
{
|
|
std::uint64_t id=0;
|
|
V3 normal{};
|
|
};
|
|
struct FieldSample
|
|
{
|
|
double density=1.0;
|
|
std::array<Plane,16> planes{};
|
|
std::uint8_t planeCount=0;
|
|
V3 gradient{};
|
|
bool sharp() const { return planeCount!=0; }
|
|
std::uint64_t primaryFeature() const
|
|
{
|
|
if(!planeCount)return 0;
|
|
std::uint64_t id=planes[0].id;
|
|
for(std::uint8_t i=1;i<planeCount;++i)id=std::min(id,planes[i].id);
|
|
return id;
|
|
}
|
|
};
|
|
struct Box
|
|
{
|
|
V3 low{},high{};
|
|
std::uint32_t id=0;
|
|
bool inverted=false;
|
|
};
|
|
struct BoxSample
|
|
{
|
|
double density=1.0;
|
|
std::array<Plane,3> planes{};
|
|
std::uint8_t planeCount=0;
|
|
V3 gradient{};
|
|
};
|
|
|
|
BoxSample EvalBox(const Box& box,V3 p)
|
|
{
|
|
const V3 center=(box.low+box.high)*0.5;
|
|
const V3 half=(box.high-box.low)*0.5;
|
|
const V3 delta=p-center;
|
|
const V3 q{std::abs(delta.x)-half.x,std::abs(delta.y)-half.y,std::abs(delta.z)-half.z};
|
|
const V3 outside{std::max(q.x,0.0),std::max(q.y,0.0),std::max(q.z,0.0)};
|
|
BoxSample out;
|
|
out.density=Length(outside)+std::min(std::max(q.x,std::max(q.y,q.z)),0.0);
|
|
const double maxQ=std::max(q.x,std::max(q.y,q.z));
|
|
V3 gradient{};
|
|
if(Length(outside)>1.0e-14)
|
|
{
|
|
gradient=Normalize({outside.x*(delta.x<0?-1:1),outside.y*(delta.y<0?-1:1),outside.z*(delta.z<0?-1:1)});
|
|
}
|
|
for(int axis=0;axis<3;++axis)
|
|
{
|
|
const double qAxis=axis==0?q.x:(axis==1?q.y:q.z);
|
|
if(std::abs(qAxis-maxQ)>1.0e-10)continue;
|
|
const double dAxis=axis==0?delta.x:(axis==1?delta.y:delta.z);
|
|
V3 normal{};
|
|
if(axis==0)normal.x=dAxis<0?-1:1;
|
|
if(axis==1)normal.y=dAxis<0?-1:1;
|
|
if(axis==2)normal.z=dAxis<0?-1:1;
|
|
const std::uint64_t side=dAxis<0?0ull:1ull;
|
|
const std::uint64_t feature=(static_cast<std::uint64_t>(box.id)<<4)
|
|
|static_cast<std::uint64_t>(axis*2)+side+1ull;
|
|
if(box.inverted)normal=normal*-1.0;
|
|
out.planes[out.planeCount++]={feature,normal};
|
|
gradient+=normal;
|
|
}
|
|
if(box.inverted)
|
|
{
|
|
out.density=-out.density;
|
|
gradient=gradient*-1.0;
|
|
}
|
|
out.gradient=Length(gradient)>1.0e-14?Normalize(gradient):Normalize(outside);
|
|
return out;
|
|
}
|
|
|
|
enum class Piece : int { Corner=0, Doorway=1, Step=2, Ledge=3, Organic=4 };
|
|
const char* PieceName(Piece p)
|
|
{
|
|
switch(p)
|
|
{
|
|
case Piece::Corner:return "corner90";
|
|
case Piece::Doorway:return "doorway";
|
|
case Piece::Step:return "riser1m";
|
|
case Piece::Ledge:return "cliff_overhang";
|
|
default:return "organic_cave_chamber";
|
|
}
|
|
}
|
|
struct Field
|
|
{
|
|
Piece piece=Piece::Corner;
|
|
V3 phaseShift{};
|
|
std::array<Box,4> boxes{};
|
|
std::uint8_t boxCount=0;
|
|
Box doorVoid{};
|
|
void AddBox(V3 low,V3 high,std::uint32_t id)
|
|
{ boxes[boxCount++]={low,high,id,false}; }
|
|
explicit Field(Piece p,V3 phase={}):piece(p),phaseShift(phase)
|
|
{
|
|
switch(p)
|
|
{
|
|
case Piece::Corner:
|
|
AddBox({-2.5,-2.5,-1.5},{0.25,2.5,1.5},1);
|
|
AddBox({-2.5,-0.25,-1.5},{2.5,2.5,1.5},2);
|
|
break;
|
|
case Piece::Doorway:
|
|
AddBox({-0.5,-2.5,-2.0},{0.5,2.5,3.0},3);
|
|
doorVoid={{-0.8,-0.85,-1.75},{0.8,0.85,2.0},4,true};
|
|
break;
|
|
case Piece::Step:
|
|
AddBox({-3.0,-3.0,-2.0},{3.0,0.25,0.0},5);
|
|
AddBox({-3.0,-0.25,-2.0},{3.0,3.0,1.0},6);
|
|
break;
|
|
case Piece::Ledge:
|
|
AddBox({-3.0,-3.0,-2.0},{0.25,3.0,3.0},7);
|
|
AddBox({-0.25,-3.0,1.5},{2.2,3.0,2.0},8);
|
|
break;
|
|
case Piece::Organic:break;
|
|
}
|
|
}
|
|
|
|
FieldSample Eval(V3 p) const
|
|
{
|
|
p=p+phaseShift;
|
|
if(piece==Piece::Organic)
|
|
{
|
|
const V3 d=p-V3{0,0,0};
|
|
FieldSample result;
|
|
// Exact spherical chamber cavity: negative solid outside, positive air inside.
|
|
result.density=2.0-Length(d);
|
|
result.gradient=Normalize(d)*-1.0;
|
|
return result;
|
|
}
|
|
FieldSample result;
|
|
result.density=std::numeric_limits<double>::infinity();
|
|
for(std::uint8_t i=0;i<boxCount;++i)
|
|
{
|
|
const BoxSample candidate=EvalBox(boxes[i],p);
|
|
if(candidate.density<result.density-1.0e-10)
|
|
{
|
|
result.density=candidate.density;
|
|
result.planeCount=candidate.planeCount;
|
|
for(std::uint8_t k=0;k<candidate.planeCount;++k)result.planes[k]=candidate.planes[k];
|
|
}
|
|
else if(std::abs(candidate.density-result.density)<=1.0e-10)
|
|
{
|
|
for(std::uint8_t k=0;k<candidate.planeCount&&result.planeCount<result.planes.size();++k)
|
|
result.planes[result.planeCount++]=candidate.planes[k];
|
|
}
|
|
}
|
|
if(piece==Piece::Doorway)
|
|
{
|
|
const BoxSample cut=EvalBox(doorVoid,p);
|
|
const double cutDensity=cut.density;
|
|
if(cutDensity>result.density+1.0e-10)
|
|
{
|
|
result.density=cutDensity;
|
|
result.planeCount=cut.planeCount;
|
|
for(std::uint8_t k=0;k<cut.planeCount;++k)result.planes[k]=cut.planes[k];
|
|
}
|
|
else if(std::abs(cutDensity-result.density)<=1.0e-10)
|
|
{
|
|
for(std::uint8_t k=0;k<cut.planeCount&&result.planeCount<result.planes.size();++k)
|
|
result.planes[result.planeCount++]=cut.planes[k];
|
|
}
|
|
}
|
|
std::sort(result.planes.begin(),result.planes.begin()+result.planeCount,
|
|
[](const Plane& a,const Plane& b){return a.id<b.id;});
|
|
std::uint8_t unique=0;
|
|
for(std::uint8_t i=0;i<result.planeCount;++i)
|
|
if(unique==0||result.planes[i].id!=result.planes[unique-1].id)
|
|
result.planes[unique++]=result.planes[i];
|
|
result.planeCount=unique;
|
|
V3 sum{};
|
|
for(std::uint8_t i=0;i<result.planeCount;++i)sum+=result.planes[i].normal;
|
|
result.gradient=Length(sum)>1.0e-14?Normalize(sum):V3{0,0,1};
|
|
return result;
|
|
}
|
|
};
|
|
|
|
struct Grid
|
|
{
|
|
int nx=0,ny=0,nz=0,step=1;
|
|
V3 origin{};
|
|
std::vector<double> density;
|
|
std::size_t NodeIndex(int x,int y,int z) const
|
|
{ return (static_cast<std::size_t>(z)*(ny+1)+y)*(nx+1)+x; }
|
|
V3 Position(int x,int y,int z) const
|
|
{ return origin+V3{x*step*VoxelMeters,y*step*VoxelMeters,z*step*VoxelMeters}; }
|
|
std::size_t CellIndex(int x,int y,int z) const
|
|
{ return (static_cast<std::size_t>(z)*ny+y)*nx+x; }
|
|
std::size_t CellCount() const { return static_cast<std::size_t>(nx)*ny*nz; }
|
|
};
|
|
Grid SampleGrid(const Field& field,int nx,int ny,int nz,int step,V3 origin,int workers)
|
|
{
|
|
Grid grid;grid.nx=nx;grid.ny=ny;grid.nz=nz;grid.step=step;grid.origin=origin;
|
|
const std::size_t count=static_cast<std::size_t>(nx+1)*(ny+1)*(nz+1);
|
|
grid.density.resize(count);
|
|
ParallelFor(count,workers,[&](std::size_t index)
|
|
{
|
|
const int x=static_cast<int>(index%(nx+1));
|
|
const int y=static_cast<int>((index/(nx+1))%(ny+1));
|
|
const int z=static_cast<int>(index/((nx+1)*(ny+1)));
|
|
const double density=field.Eval(grid.Position(x,y,z)).density;
|
|
// Stable iso-level tie rule: zero remains air, matching the Core sign convention, but
|
|
// receives a sub-cell positive bias so edge roots do not collapse onto one lattice
|
|
// vertex from several unrelated edges. This is only 0.01 mm and is below the output
|
|
// surface-error precision; it is shared by all three candidates.
|
|
grid.density[index]=density==0.0?1.0e-5:density;
|
|
});
|
|
return grid;
|
|
}
|
|
|
|
struct Root
|
|
{
|
|
bool active=false;
|
|
V3 position{},gradient{};
|
|
std::array<Plane,16> planes{};
|
|
std::uint8_t planeCount=0;
|
|
std::uint64_t feature=0;
|
|
};
|
|
struct EdgeArrays
|
|
{
|
|
std::array<std::vector<Root>,3> roots;
|
|
std::array<std::vector<int>,3> vertexIds;
|
|
std::array<std::size_t,3> counts{};
|
|
std::size_t Total() const{return counts[0]+counts[1]+counts[2];}
|
|
std::size_t Index(int axis,int x,int y,int z,const Grid& g) const
|
|
{
|
|
if(axis==0)return (static_cast<std::size_t>(z)*(g.ny+1)+y)*g.nx+x;
|
|
if(axis==1)return (static_cast<std::size_t>(z)*g.ny+y)*(g.nx+1)+x;
|
|
return (static_cast<std::size_t>(z)*(g.ny+1)+y)*(g.nx+1)+x;
|
|
}
|
|
Root& Get(int axis,int x,int y,int z,const Grid& g){return roots[axis][Index(axis,x,y,z,g)];}
|
|
const Root& Get(int axis,int x,int y,int z,const Grid& g) const{return roots[axis][Index(axis,x,y,z,g)];}
|
|
int& Vertex(int axis,int x,int y,int z,const Grid& g){return vertexIds[axis][Index(axis,x,y,z,g)];}
|
|
int Vertex(int axis,int x,int y,int z,const Grid& g) const{return vertexIds[axis][Index(axis,x,y,z,g)];}
|
|
};
|
|
|
|
EdgeArrays BuildRoots(const Grid& grid,const Field& field,int workers)
|
|
{
|
|
EdgeArrays edges;
|
|
edges.counts={static_cast<std::size_t>(grid.nx)*(grid.ny+1)*(grid.nz+1),
|
|
static_cast<std::size_t>(grid.nx+1)*grid.ny*(grid.nz+1),
|
|
static_cast<std::size_t>(grid.nx+1)*(grid.ny+1)*grid.nz};
|
|
for(int axis=0;axis<3;++axis)
|
|
{
|
|
edges.roots[axis].resize(edges.counts[axis]);
|
|
edges.vertexIds[axis].assign(edges.counts[axis],-1);
|
|
}
|
|
const std::size_t total=edges.Total();
|
|
ParallelFor(total,workers,[&](std::size_t flat)
|
|
{
|
|
int axis=0;std::size_t index=flat;
|
|
while(axis<2&&index>=edges.counts[axis]){index-=edges.counts[axis];++axis;}
|
|
int x=0,y=0,z=0;
|
|
if(axis==0)
|
|
{
|
|
x=static_cast<int>(index%grid.nx);y=static_cast<int>((index/grid.nx)%(grid.ny+1));
|
|
z=static_cast<int>(index/(grid.nx*(grid.ny+1)));
|
|
}
|
|
else if(axis==1)
|
|
{
|
|
x=static_cast<int>(index%(grid.nx+1));y=static_cast<int>((index/(grid.nx+1))%grid.ny);
|
|
z=static_cast<int>(index/((grid.nx+1)*grid.ny));
|
|
}
|
|
else
|
|
{
|
|
x=static_cast<int>(index%(grid.nx+1));y=static_cast<int>((index/(grid.nx+1))%(grid.ny+1));
|
|
z=static_cast<int>(index/((grid.nx+1)*(grid.ny+1)));
|
|
}
|
|
int x2=x,y2=y,z2=z;
|
|
if(axis==0)++x2;else if(axis==1)++y2;else ++z2;
|
|
const double d0=grid.density[grid.NodeIndex(x,y,z)];
|
|
const double d1=grid.density[grid.NodeIndex(x2,y2,z2)];
|
|
if((d0<0.0)==(d1<0.0))return;
|
|
double t=d0/(d0-d1);t=std::clamp(t,0.0,1.0);
|
|
Root root;root.active=true;
|
|
root.position=Quantize(grid.Position(x,y,z)+(grid.Position(x2,y2,z2)-grid.Position(x,y,z))*t);
|
|
const FieldSample sample=field.Eval(root.position);
|
|
const double h=grid.step*VoxelMeters;
|
|
const auto gradientAt=[&](V3 point)
|
|
{
|
|
const double dx=field.Eval(point+V3{h,0,0}).density-field.Eval(point-V3{h,0,0}).density;
|
|
const double dy=field.Eval(point+V3{0,h,0}).density-field.Eval(point-V3{0,h,0}).density;
|
|
const double dz=field.Eval(point+V3{0,0,h}).density-field.Eval(point-V3{0,0,h}).density;
|
|
return Normalize({dx,dy,dz});
|
|
};
|
|
const V3 g0=gradientAt(grid.Position(x,y,z));
|
|
const V3 g1=gradientAt(grid.Position(x2,y2,z2));
|
|
root.gradient=Normalize(g0*(1.0-t)+g1*t);
|
|
root.planeCount=sample.planeCount;
|
|
for(std::uint8_t i=0;i<root.planeCount;++i)root.planes[i]=sample.planes[i];
|
|
root.feature=sample.primaryFeature();
|
|
edges.roots[axis][index]=root;
|
|
});
|
|
return edges;
|
|
}
|
|
|
|
struct Vertex
|
|
{
|
|
V3 position{},normal{};
|
|
};
|
|
struct Triangle
|
|
{
|
|
std::array<std::uint32_t,3> index{};
|
|
std::uint64_t feature=0;
|
|
};
|
|
struct Mesh
|
|
{
|
|
std::vector<Vertex> vertices;
|
|
std::vector<Triangle> triangles;
|
|
std::uint64_t hash=0;
|
|
std::uint64_t skippedDualEdges=0;
|
|
};
|
|
enum class Method : int { MC=0, Hybrid=1, DC=2 };
|
|
const char* MethodName(Method m){return m==Method::MC?"MC":m==Method::Hybrid?"Hybrid": "DC";}
|
|
|
|
struct Cell
|
|
{
|
|
std::uint8_t cubeCase=0;
|
|
bool active=false;
|
|
std::uint8_t componentCount=0;
|
|
std::array<std::uint8_t,EdgeCount> edgeComponent{};
|
|
std::array<bool,EdgeCount> componentSharp{},componentUseDC{};
|
|
std::array<V3,EdgeCount> dual{};
|
|
};
|
|
|
|
int FindLocal(std::array<std::uint8_t,EdgeCount>& parent,int value)
|
|
{
|
|
while(parent[value]!=value)
|
|
{
|
|
parent[value]=parent[parent[value]];
|
|
value=parent[value];
|
|
}
|
|
return value;
|
|
}
|
|
void UniteLocal(std::array<std::uint8_t,EdgeCount>& parent,int a,int b)
|
|
{
|
|
a=FindLocal(parent,a);b=FindLocal(parent,b);if(a==b)return;
|
|
if(a>b)std::swap(a,b);
|
|
parent[b]=static_cast<std::uint8_t>(a);
|
|
}
|
|
|
|
struct ComponentSets
|
|
{
|
|
std::vector<std::uint32_t> parent;
|
|
std::vector<std::uint8_t> rank,sharp;
|
|
explicit ComponentSets(std::size_t size):parent(size),rank(size,0),sharp(size,0)
|
|
{std::iota(parent.begin(),parent.end(),0u);}
|
|
std::uint32_t Find(std::uint32_t value)
|
|
{
|
|
while(parent[value]!=value){parent[value]=parent[parent[value]];value=parent[value];}
|
|
return value;
|
|
}
|
|
void Unite(std::uint32_t a,std::uint32_t b)
|
|
{
|
|
a=Find(a);b=Find(b);if(a==b)return;
|
|
if(rank[a]<rank[b])std::swap(a,b);
|
|
parent[b]=a;sharp[a]=static_cast<std::uint8_t>(sharp[a]|sharp[b]);
|
|
if(rank[a]==rank[b])++rank[a];
|
|
}
|
|
};
|
|
|
|
std::array<int,4> LocalEdge(int x,int y,int z,int edge)
|
|
{
|
|
const auto c=EdgeCorners[edge];const auto a=CornerOffsets[c[0]],b=CornerOffsets[c[1]];
|
|
const int axis=a[0]!=b[0]?0:a[1]!=b[1]?1:2;
|
|
return {axis,x+std::min(a[0],b[0]),y+std::min(a[1],b[1]),z+std::min(a[2],b[2])};
|
|
}
|
|
int FindLocalEdge(int cellX,int cellY,int cellZ,int axis,int edgeX,int edgeY,int edgeZ)
|
|
{
|
|
for(int edge=0;edge<EdgeCount;++edge)
|
|
{
|
|
const auto local=LocalEdge(cellX,cellY,cellZ,edge);
|
|
if(local[0]==axis&&local[1]==edgeX&&local[2]==edgeY&&local[3]==edgeZ)return edge;
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
V3 SolveQEF(const std::vector<std::pair<V3,V3>>& constraints,V3 fallback,V3 low,V3 high)
|
|
{
|
|
if(constraints.empty())return fallback;
|
|
double a[3][4]{};
|
|
// Fixed constraint order; a small center penalty resolves rank-deficient planes without
|
|
// moving a two/three-plane intersection at measurable scale.
|
|
constexpr double lambda=1.0e-10;
|
|
const double c[3]={fallback.x,fallback.y,fallback.z};
|
|
for(int i=0;i<3;++i){a[i][i]=lambda;a[i][3]=lambda*c[i];}
|
|
for(const auto& constraint:constraints)
|
|
{
|
|
const V3 n=Normalize(constraint.second);const double b=Dot(n,constraint.first);
|
|
const double nv[3]={n.x,n.y,n.z};
|
|
for(int r=0;r<3;++r){for(int k=0;k<3;++k)a[r][k]+=nv[r]*nv[k];a[r][3]+=nv[r]*b;}
|
|
}
|
|
bool singular=false;
|
|
for(int column=0;column<3;++column)
|
|
{
|
|
int pivot=column;
|
|
for(int row=column+1;row<3;++row)if(std::abs(a[row][column])>std::abs(a[pivot][column]))pivot=row;
|
|
if(std::abs(a[pivot][column])<1.0e-14){singular=true;break;}
|
|
if(pivot!=column)for(int k=column;k<4;++k)std::swap(a[pivot][k],a[column][k]);
|
|
const double divisor=a[column][column];
|
|
for(int k=column;k<4;++k)a[column][k]/=divisor;
|
|
for(int row=0;row<3;++row)if(row!=column)
|
|
{
|
|
const double factor=a[row][column];
|
|
for(int k=column;k<4;++k)a[row][k]-=factor*a[column][k];
|
|
}
|
|
}
|
|
if(singular)return fallback;
|
|
return Quantize({std::clamp(a[0][3],low.x,high.x),
|
|
std::clamp(a[1][3],low.y,high.y),std::clamp(a[2][3],low.z,high.z)});
|
|
}
|
|
|
|
std::uint64_t SelectFeature(const Root& root,int edgeAxis)
|
|
{
|
|
if(!root.planeCount)return 0;
|
|
std::uint8_t selected=0;
|
|
for(std::uint8_t i=1;i<root.planeCount;++i)
|
|
{
|
|
const double a=edgeAxis==0?std::abs(root.planes[i].normal.x):
|
|
edgeAxis==1?std::abs(root.planes[i].normal.y):std::abs(root.planes[i].normal.z);
|
|
const double b=edgeAxis==0?std::abs(root.planes[selected].normal.x):
|
|
edgeAxis==1?std::abs(root.planes[selected].normal.y):std::abs(root.planes[selected].normal.z);
|
|
if(a>b+1.0e-12||(std::abs(a-b)<=1.0e-12&&root.planes[i].id<root.planes[selected].id))selected=i;
|
|
}
|
|
return root.planes[selected].id;
|
|
}
|
|
V3 NormalForFeature(const Root& root,std::uint64_t id)
|
|
{
|
|
for(std::uint8_t i=0;i<root.planeCount;++i)if(root.planes[i].id==id)return root.planes[i].normal;
|
|
return root.gradient;
|
|
}
|
|
|
|
void AddTriangle(Mesh& mesh,std::uint32_t a,std::uint32_t b,std::uint32_t c,
|
|
std::uint64_t feature,V3 outward)
|
|
{
|
|
if(a==b||b==c||c==a)return;
|
|
const V3 pa=mesh.vertices[a].position,pb=mesh.vertices[b].position,pc=mesh.vertices[c].position;
|
|
if(Length(Cross(pb-pa,pc-pa))<1.0e-14)return;
|
|
if(Dot(Cross(pb-pa,pc-pa),outward)<0.0)std::swap(b,c);
|
|
mesh.triangles.push_back({{a,b,c},feature});
|
|
}
|
|
|
|
std::uint64_t HashMesh(const Mesh& mesh)
|
|
{
|
|
std::uint64_t h=14695981039346656037ull;
|
|
const auto add=[&](std::uint64_t word,int bytes)
|
|
{
|
|
for(int i=0;i<bytes;++i){h^=static_cast<std::uint8_t>(word>>(i*8));h*=1099511628211ull;}
|
|
};
|
|
add(mesh.vertices.size(),8);
|
|
for(const Vertex& v:mesh.vertices)
|
|
{
|
|
for(double x:{v.position.x,v.position.y,v.position.z,v.normal.x,v.normal.y,v.normal.z})
|
|
add(std::bit_cast<std::uint64_t>(x),8);
|
|
}
|
|
add(mesh.triangles.size(),8);
|
|
for(const Triangle& t:mesh.triangles)
|
|
{add(t.index[0],4);add(t.index[1],4);add(t.index[2],4);add(t.feature,8);}
|
|
return h;
|
|
}
|
|
|
|
Mesh BuildMesh(const Grid& grid,const Field& field,Method method,int workers)
|
|
{
|
|
EdgeArrays edges=BuildRoots(grid,field,workers);
|
|
std::vector<Cell> cells(grid.CellCount());
|
|
ParallelFor(cells.size(),workers,[&](std::size_t ci)
|
|
{
|
|
const int x=static_cast<int>(ci%grid.nx);
|
|
const int y=static_cast<int>((ci/grid.nx)%grid.ny);
|
|
const int z=static_cast<int>(ci/(grid.nx*grid.ny));
|
|
Cell& cell=cells[ci];cell.edgeComponent.fill(0xff);
|
|
double values[CornerCount];
|
|
for(int corner=0;corner<CornerCount;++corner)
|
|
{
|
|
const auto o=CornerOffsets[corner];
|
|
values[corner]=grid.density[grid.NodeIndex(x+o[0],y+o[1],z+o[2])];
|
|
if(values[corner]<0.0)cell.cubeCase|=static_cast<std::uint8_t>(1u<<corner);
|
|
}
|
|
if(cell.cubeCase==0||cell.cubeCase==255)return;
|
|
cell.active=true;
|
|
if(method==Method::MC)
|
|
{
|
|
cell.componentCount=1;
|
|
for(int edge=0;edge<EdgeCount;++edge)
|
|
{
|
|
const auto ec=EdgeCorners[edge];
|
|
if((values[ec[0]]<0.0)!=(values[ec[1]]<0.0))cell.edgeComponent[edge]=0;
|
|
}
|
|
return;
|
|
}
|
|
std::array<std::uint8_t,EdgeCount> localParent{};
|
|
for(int edge=0;edge<EdgeCount;++edge)localParent[edge]=static_cast<std::uint8_t>(edge);
|
|
for(int k=0;k<16&&MCTriTable[cell.cubeCase][k]>=0;k+=3)
|
|
{
|
|
const int a=MCTriTable[cell.cubeCase][k],b=MCTriTable[cell.cubeCase][k+1],c=MCTriTable[cell.cubeCase][k+2];
|
|
UniteLocal(localParent,a,b);UniteLocal(localParent,b,c);
|
|
}
|
|
std::array<std::uint8_t,EdgeCount> rootToComponent{};rootToComponent.fill(0xff);
|
|
for(int edge=0;edge<EdgeCount;++edge)
|
|
{
|
|
const auto ec=EdgeCorners[edge];
|
|
if((values[ec[0]]<0.0)==(values[ec[1]]<0.0))continue;
|
|
const auto local=LocalEdge(x,y,z,edge);
|
|
const Root& root=edges.Get(local[0],local[1],local[2],local[3],grid);
|
|
if(!root.active)continue;
|
|
const int localRoot=FindLocal(localParent,edge);
|
|
if(rootToComponent[localRoot]==0xff)
|
|
rootToComponent[localRoot]=cell.componentCount++;
|
|
const std::uint8_t component=rootToComponent[localRoot];
|
|
cell.edgeComponent[edge]=component;
|
|
cell.componentSharp[component]|=root.planeCount!=0;
|
|
}
|
|
});
|
|
|
|
// Link local surface patches only through the same sign-changing primal edge. A cube may
|
|
// contain multiple disconnected patches; keeping their QEFs separate avoids DC pinches.
|
|
std::unique_ptr<ComponentSets> patches;
|
|
if(method==Method::Hybrid)
|
|
{
|
|
patches=std::make_unique<ComponentSets>(cells.size()*EdgeCount);
|
|
for(std::size_t cellIndex=0;cellIndex<cells.size();++cellIndex)
|
|
for(std::uint8_t component=0;component<cells[cellIndex].componentCount;++component)
|
|
patches->sharp[cellIndex*EdgeCount+component]=cells[cellIndex].componentSharp[component]?1:0;
|
|
for(int axis=0;axis<3;++axis)
|
|
{
|
|
for(std::size_t i=0;i<edges.roots[axis].size();++i)
|
|
{
|
|
if(!edges.roots[axis][i].active)continue;
|
|
int x=0,y=0,z=0;
|
|
if(axis==0){x=static_cast<int>(i%grid.nx);y=static_cast<int>((i/grid.nx)%(grid.ny+1));z=static_cast<int>(i/(grid.nx*(grid.ny+1)));}
|
|
else if(axis==1){x=static_cast<int>(i%(grid.nx+1));y=static_cast<int>((i/(grid.nx+1))%grid.ny);z=static_cast<int>(i/((grid.nx+1)*grid.ny));}
|
|
else{x=static_cast<int>(i%(grid.nx+1));y=static_cast<int>((i/(grid.nx+1))%(grid.ny+1));z=static_cast<int>(i/((grid.nx+1)*(grid.ny+1)));}
|
|
std::array<std::array<int,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<std::uint32_t,4> patchIds{};bool valid=true;
|
|
for(int corner=0;corner<4;++corner)
|
|
{
|
|
const auto& c=around[corner];
|
|
if(c[0]<0||c[1]<0||c[2]<0||c[0]>=grid.nx||c[1]>=grid.ny||c[2]>=grid.nz){valid=false;break;}
|
|
const std::size_t cellIndex=grid.CellIndex(c[0],c[1],c[2]);
|
|
const int localEdge=FindLocalEdge(c[0],c[1],c[2],axis,x,y,z);
|
|
if(localEdge<0||cells[cellIndex].edgeComponent[localEdge]==0xff){valid=false;break;}
|
|
patchIds[corner]=static_cast<std::uint32_t>(cellIndex*EdgeCount+cells[cellIndex].edgeComponent[localEdge]);
|
|
}
|
|
if(valid)for(int corner=1;corner<4;++corner)patches->Unite(patchIds[0],patchIds[corner]);
|
|
}
|
|
}
|
|
}
|
|
for(std::size_t cellIndex=0;cellIndex<cells.size();++cellIndex)
|
|
for(std::uint8_t component=0;component<cells[cellIndex].componentCount;++component)
|
|
{
|
|
const std::uint32_t patch=static_cast<std::uint32_t>(cellIndex*EdgeCount+component);
|
|
cells[cellIndex].componentUseDC[component]=method==Method::DC||
|
|
(method==Method::Hybrid&&patches->sharp[patches->Find(patch)]!=0);
|
|
}
|
|
|
|
ParallelFor(cells.size(),workers,[&](std::size_t ci)
|
|
{
|
|
Cell& cell=cells[ci];if(cell.componentCount==0)return;
|
|
const int x=static_cast<int>(ci%grid.nx);
|
|
const int y=static_cast<int>((ci/grid.nx)%grid.ny);
|
|
const int z=static_cast<int>(ci/(grid.nx*grid.ny));
|
|
for(std::uint8_t component=0;component<cell.componentCount;++component)
|
|
{
|
|
if(!cell.componentUseDC[component])continue;
|
|
std::vector<std::pair<V3,V3>> constraints;
|
|
V3 rootSum{};int rootCount=0;
|
|
for(int edge=0;edge<EdgeCount;++edge)
|
|
{
|
|
if(cell.edgeComponent[edge]!=component)continue;
|
|
const auto local=LocalEdge(x,y,z,edge);
|
|
const Root& root=edges.Get(local[0],local[1],local[2],local[3],grid);
|
|
if(!root.active)continue;
|
|
rootSum+=root.position;++rootCount;
|
|
if(root.planeCount)
|
|
{
|
|
for(std::uint8_t i=0;i<root.planeCount;++i)
|
|
constraints.emplace_back(root.position,root.planes[i].normal);
|
|
}
|
|
else constraints.emplace_back(root.position,root.gradient);
|
|
}
|
|
const V3 low=grid.Position(x,y,z),high=grid.Position(x+1,y+1,z+1);
|
|
const V3 center=rootCount?rootSum/static_cast<double>(rootCount):(low+high)*0.5;
|
|
cell.dual[component]=SolveQEF(constraints,center,low,high);
|
|
}
|
|
});
|
|
|
|
Mesh mesh;
|
|
std::unordered_map<std::uint64_t,std::uint32_t> dcVertexMap;
|
|
if(method==Method::MC||method==Method::Hybrid)
|
|
{
|
|
std::map<std::tuple<std::int64_t,std::int64_t,std::int64_t>,std::uint32_t> positionMap;
|
|
for(int axis=0;axis<3;++axis)for(std::size_t index=0;index<edges.roots[axis].size();++index)
|
|
{
|
|
const Root& root=edges.roots[axis][index];if(!root.active)continue;
|
|
const auto q=[](double x){return static_cast<std::int64_t>(std::llround(x/PositionQuantum));};
|
|
const V3 p=root.position;
|
|
const auto key=std::tuple{q(p.x),q(p.y),q(p.z)};
|
|
auto found=positionMap.find(key);
|
|
std::uint32_t vertex=0;
|
|
if(found==positionMap.end())
|
|
{
|
|
vertex=static_cast<std::uint32_t>(mesh.vertices.size());
|
|
mesh.vertices.push_back({p,root.gradient});positionMap.emplace(key,vertex);
|
|
}
|
|
else vertex=found->second;
|
|
edges.vertexIds[axis][index]=static_cast<int>(vertex);
|
|
}
|
|
// Coincident zero-at-node roots use the deterministic first edge normal. The shared
|
|
// vertex position, matching the owner's 0.01 cm-style MC welding contract.
|
|
}
|
|
const auto getDCVertex=[&](int cellId,std::uint8_t component,const Root& root,int axis)->std::uint32_t
|
|
{
|
|
const std::uint64_t feature=SelectFeature(root,axis);
|
|
const std::uint64_t patch=static_cast<std::uint64_t>(cellId)*EdgeCount+component;
|
|
const std::uint64_t key=(patch<<32)^(feature&0xffffffffull);
|
|
const auto existing=dcVertexMap.find(key);
|
|
if(existing!=dcVertexMap.end())return existing->second;
|
|
const V3 normal=feature?NormalForFeature(root,feature):root.gradient;
|
|
const std::uint32_t index=static_cast<std::uint32_t>(mesh.vertices.size());
|
|
mesh.vertices.push_back({cells[cellId].dual[component],Normalize(normal)});
|
|
dcVertexMap.emplace(key,index);
|
|
return index;
|
|
};
|
|
for(int z=0;z<grid.nz;++z)for(int y=0;y<grid.ny;++y)for(int x=0;x<grid.nx;++x)
|
|
{
|
|
const std::size_t ci=grid.CellIndex(x,y,z);const Cell& cell=cells[ci];
|
|
if(!cell.active)continue;
|
|
std::array<int,EdgeCount> localVertices{};localVertices.fill(-1);
|
|
for(int edge=0;edge<EdgeCount;++edge)
|
|
{
|
|
const auto ec=EdgeCorners[edge];
|
|
const auto oa=CornerOffsets[ec[0]],ob=CornerOffsets[ec[1]];
|
|
const double da=grid.density[grid.NodeIndex(x+oa[0],y+oa[1],z+oa[2])];
|
|
const double db=grid.density[grid.NodeIndex(x+ob[0],y+ob[1],z+ob[2])];
|
|
if((da<0.0)==(db<0.0))continue;
|
|
const auto e=LocalEdge(x,y,z,edge);
|
|
localVertices[edge]=edges.Vertex(e[0],e[1],e[2],e[3],grid);
|
|
}
|
|
for(int k=0;k<16&&MCTriTable[cell.cubeCase][k]>=0;k+=3)
|
|
{
|
|
const int e0=MCTriTable[cell.cubeCase][k],e1=MCTriTable[cell.cubeCase][k+1],e2=MCTriTable[cell.cubeCase][k+2];
|
|
if(localVertices[e0]<0||localVertices[e1]<0||localVertices[e2]<0)continue;
|
|
const std::uint8_t component=cell.edgeComponent[e0];
|
|
if(component==0xff||cell.componentUseDC[component])continue;
|
|
std::array<std::uint64_t,3> labels{};
|
|
const std::array<int,3> triEdges{{e0,e1,e2}};
|
|
for(int edgeIndex=0;edgeIndex<3;++edgeIndex)
|
|
{
|
|
const auto le=LocalEdge(x,y,z,triEdges[edgeIndex]);
|
|
const Root& root=edges.Get(le[0],le[1],le[2],le[3],grid);
|
|
labels[edgeIndex]=root.feature;
|
|
}
|
|
std::uint64_t label=0;int bestVotes=0;
|
|
for(std::uint64_t candidate:labels)if(candidate)
|
|
{
|
|
const int votes=static_cast<int>(std::count(labels.begin(),labels.end(),candidate));
|
|
if(votes>bestVotes||(votes==bestVotes&&(!label||candidate<label)))
|
|
{label=candidate;bestVotes=votes;}
|
|
}
|
|
const V3 center=(mesh.vertices[localVertices[e0]].position+
|
|
mesh.vertices[localVertices[e1]].position+mesh.vertices[localVertices[e2]].position)/3.0;
|
|
AddTriangle(mesh,static_cast<std::uint32_t>(localVertices[e0]),
|
|
static_cast<std::uint32_t>(localVertices[e1]),static_cast<std::uint32_t>(localVertices[e2]),
|
|
label,field.Eval(center).gradient);
|
|
}
|
|
}
|
|
for(int axis=0;axis<3;++axis)
|
|
{
|
|
for(std::size_t index=0;index<edges.roots[axis].size();++index)
|
|
{
|
|
const Root& root=edges.roots[axis][index];if(!root.active)continue;
|
|
int x=0,y=0,z=0;
|
|
if(axis==0){x=static_cast<int>(index%grid.nx);y=static_cast<int>((index/grid.nx)%(grid.ny+1));z=static_cast<int>(index/(grid.nx*(grid.ny+1)));}
|
|
else if(axis==1){x=static_cast<int>(index%(grid.nx+1));y=static_cast<int>((index/(grid.nx+1))%grid.ny);z=static_cast<int>(index/((grid.nx+1)*grid.ny));}
|
|
else{x=static_cast<int>(index%(grid.nx+1));y=static_cast<int>((index/(grid.nx+1))%(grid.ny+1));z=static_cast<int>(index/((grid.nx+1)*(grid.ny+1)));}
|
|
std::array<std::array<int,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<int,4> incident{};bool valid=true;bool anyDC=false;
|
|
std::array<std::uint8_t,4> incidentComponents{};
|
|
for(int i=0;i<4;++i)
|
|
{
|
|
const auto& c=around[i];
|
|
if(c[0]<0||c[1]<0||c[2]<0||c[0]>=grid.nx||c[1]>=grid.ny||c[2]>=grid.nz)
|
|
{valid=false;break;}
|
|
incident[i]=static_cast<int>(grid.CellIndex(c[0],c[1],c[2]));
|
|
const int localEdge=FindLocalEdge(c[0],c[1],c[2],axis,x,y,z);
|
|
if(localEdge<0){valid=false;break;}
|
|
incidentComponents[i]=cells[incident[i]].edgeComponent[localEdge];
|
|
if(incidentComponents[i]==0xff){valid=false;break;}
|
|
anyDC|=cells[incident[i]].componentUseDC[incidentComponents[i]];
|
|
}
|
|
if(!valid||!anyDC||!cells[incident[0]].componentUseDC[incidentComponents[0]])continue;
|
|
std::array<std::uint32_t,4> v{};
|
|
for(int i=0;i<4;++i)
|
|
{
|
|
if(!cells[incident[i]].componentUseDC[incidentComponents[i]]){valid=false;break;}
|
|
v[i]=getDCVertex(incident[i],incidentComponents[i],root,axis);
|
|
}
|
|
if(!valid){++mesh.skippedDualEdges;continue;}
|
|
if(Dot(Cross(mesh.vertices[v[1]].position-mesh.vertices[v[0]].position,
|
|
mesh.vertices[v[2]].position-mesh.vertices[v[0]].position),root.gradient)<0.0)
|
|
std::swap(v[1],v[3]);
|
|
const std::uint64_t label=SelectFeature(root,axis);
|
|
const std::uint32_t edgeTie=static_cast<std::uint32_t>(x+y+z+axis)&1u;
|
|
if(edgeTie==0)
|
|
{
|
|
AddTriangle(mesh,v[0],v[1],v[2],label,root.gradient);
|
|
AddTriangle(mesh,v[0],v[2],v[3],label,root.gradient);
|
|
}
|
|
else
|
|
{
|
|
AddTriangle(mesh,v[0],v[1],v[3],label,root.gradient);
|
|
AddTriangle(mesh,v[1],v[2],v[3],label,root.gradient);
|
|
}
|
|
}
|
|
}
|
|
mesh.hash=HashMesh(mesh);
|
|
return mesh;
|
|
}
|
|
|
|
struct ErrorMetrics{double mean=0,maximum=0;std::uint64_t samples=0;};
|
|
ErrorMetrics SurfaceError(const Mesh& mesh,const Field& field)
|
|
{
|
|
ErrorMetrics out;double total=0;
|
|
const auto sample=[&](V3 p)
|
|
{
|
|
const double error=std::abs(field.Eval(p).density);
|
|
total+=error;out.maximum=std::max(out.maximum,error);++out.samples;
|
|
};
|
|
for(const Triangle& t:mesh.triangles)
|
|
{
|
|
const V3 a=mesh.vertices[t.index[0]].position,b=mesh.vertices[t.index[1]].position,c=mesh.vertices[t.index[2]].position;
|
|
sample(a);sample(b);sample(c);sample((a+b)*0.5);sample((b+c)*0.5);sample((c+a)*0.5);sample((a+b+c)/3.0);
|
|
}
|
|
if(out.samples)out.mean=total/out.samples;
|
|
return out;
|
|
}
|
|
|
|
struct GeomKey
|
|
{
|
|
std::int64_t x=0,y=0,z=0;
|
|
friend bool operator<(const GeomKey& a,const GeomKey& b)
|
|
{return std::tie(a.x,a.y,a.z)<std::tie(b.x,b.y,b.z);}
|
|
friend bool operator==(const GeomKey&,const GeomKey&)=default;
|
|
};
|
|
GeomKey Key(V3 p)
|
|
{return {static_cast<std::int64_t>(std::llround(p.x/PositionQuantum)),
|
|
static_cast<std::int64_t>(std::llround(p.y/PositionQuantum)),
|
|
static_cast<std::int64_t>(std::llround(p.z/PositionQuantum))};}
|
|
struct EdgeInfo{int count=0;std::vector<std::pair<int,V3>> faces;};
|
|
std::map<std::pair<GeomKey,GeomKey>,EdgeInfo> GeomEdges(const Mesh& mesh)
|
|
{
|
|
std::map<std::pair<GeomKey,GeomKey>,EdgeInfo> result;
|
|
for(std::size_t ti=0;ti<mesh.triangles.size();++ti)
|
|
{
|
|
const Triangle& t=mesh.triangles[ti];
|
|
const V3 a=mesh.vertices[t.index[0]].position,b=mesh.vertices[t.index[1]].position,c=mesh.vertices[t.index[2]].position;
|
|
const V3 normal=Normalize(Cross(b-a,c-a));
|
|
const std::array<std::pair<std::uint32_t,std::uint32_t>,3> edges{{{t.index[0],t.index[1]},{t.index[1],t.index[2]},{t.index[2],t.index[0]}}};
|
|
for(const auto& e:edges)
|
|
{
|
|
GeomKey ka=Key(mesh.vertices[e.first].position),kb=Key(mesh.vertices[e.second].position);
|
|
if(kb<ka)std::swap(ka,kb);
|
|
auto& info=result[{ka,kb}];++info.count;info.faces.emplace_back(static_cast<int>(ti),normal);
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
struct TopologyMetrics
|
|
{
|
|
std::uint64_t boundaryEdges=0,nonManifoldEdges=0,creaseEdgePairs=0,splitPairs=0;
|
|
};
|
|
TopologyMetrics AnalyzeTopology(const Mesh& mesh)
|
|
{
|
|
TopologyMetrics out;
|
|
const auto edges=GeomEdges(mesh);
|
|
for(const auto& [key,info]:edges)
|
|
{
|
|
if(info.count==1)++out.boundaryEdges;
|
|
else if(info.count!=2)++out.nonManifoldEdges;
|
|
if(info.faces.size()==2)
|
|
{
|
|
const Triangle& a=mesh.triangles[info.faces[0].first];
|
|
const Triangle& b=mesh.triangles[info.faces[1].first];
|
|
if(a.feature&&b.feature&&a.feature!=b.feature&&
|
|
Dot(Normalize(info.faces[0].second),Normalize(info.faces[1].second))<NormalSplitCos)
|
|
{
|
|
++out.creaseEdgePairs;
|
|
const auto avgNormal=[&](const Triangle& t)
|
|
{return Normalize(mesh.vertices[t.index[0]].normal+mesh.vertices[t.index[1]].normal+mesh.vertices[t.index[2]].normal);};
|
|
if(Dot(avgNormal(a),avgNormal(b))<NormalSplitCos)++out.splitPairs;
|
|
}
|
|
}
|
|
}
|
|
return out;
|
|
}
|
|
|
|
std::uint64_t FnvBytes(const std::uint8_t* p,std::size_t n)
|
|
{
|
|
std::uint64_t h=14695981039346656037ull;
|
|
for(std::size_t i=0;i<n;++i){h^=p[i];h*=1099511628211ull;}
|
|
return h;
|
|
}
|
|
|
|
void WriteObj(const std::filesystem::path& path,const Mesh& mesh,double sideFilter=0.0,int side=0)
|
|
{
|
|
std::ofstream out(path,std::ios::binary);
|
|
out<<std::setprecision(9);
|
|
for(const Vertex& v:mesh.vertices)out<<"v "<<v.position.x<<' '<<v.position.y<<' '<<v.position.z<<"\n";
|
|
for(const Vertex& v:mesh.vertices)out<<"vn "<<v.normal.x<<' '<<v.normal.y<<' '<<v.normal.z<<"\n";
|
|
for(const Triangle& t:mesh.triangles)
|
|
{
|
|
const V3 center=(mesh.vertices[t.index[0]].position+mesh.vertices[t.index[1]].position+mesh.vertices[t.index[2]].position)/3.0;
|
|
if(side<0&¢er.x>sideFilter)continue;
|
|
if(side>0&¢er.x<sideFilter)continue;
|
|
out<<"f "<<t.index[0]+1<<"//"<<t.index[0]+1<<' '
|
|
<<t.index[1]+1<<"//"<<t.index[1]+1<<' '<<t.index[2]+1<<"//"<<t.index[2]+1<<"\n";
|
|
}
|
|
}
|
|
|
|
struct Segment{V3 a{},b{};};
|
|
std::vector<Segment> SeamSegments(const Mesh& mesh,double xPlane,int side)
|
|
{
|
|
std::vector<Segment> segments;
|
|
for(const Triangle& triangle:mesh.triangles)
|
|
{
|
|
std::array<V3,3> input{{mesh.vertices[triangle.index[0]].position,
|
|
mesh.vertices[triangle.index[1]].position,mesh.vertices[triangle.index[2]].position}};
|
|
std::vector<V3> output;
|
|
for(int i=0;i<3;++i)
|
|
{
|
|
const V3 a=input[i],b=input[(i+1)%3];
|
|
const bool ina=side<0?a.x<=xPlane+PositionQuantum:a.x>=xPlane-PositionQuantum;
|
|
const bool inb=side<0?b.x<=xPlane+PositionQuantum:b.x>=xPlane-PositionQuantum;
|
|
if(ina)output.push_back(a);
|
|
if(ina!=inb)
|
|
{
|
|
const double t=(xPlane-a.x)/(b.x-a.x);
|
|
output.push_back(Quantize(a+(b-a)*t));
|
|
}
|
|
}
|
|
std::vector<V3> onPlane;
|
|
for(V3 p:output)if(std::abs(p.x-xPlane)<=2.0*PositionQuantum)
|
|
{
|
|
bool found=false;for(V3 q:onPlane)if(Length(p-q)<PositionQuantum){found=true;break;}
|
|
if(!found)onPlane.push_back(p);
|
|
}
|
|
if(onPlane.size()==2&&Length(onPlane[0]-onPlane[1])>PositionQuantum)
|
|
segments.push_back({onPlane[0],onPlane[1]});
|
|
}
|
|
return segments;
|
|
}
|
|
double PointSegmentDistance(V3 p,const Segment& s)
|
|
{
|
|
const V3 d=s.b-s.a;const double den=Dot(d,d);
|
|
const double t=den>0?std::clamp(Dot(p-s.a,d)/den,0.0,1.0):0.0;
|
|
return Length(p-(s.a+d*t));
|
|
}
|
|
struct SeamMetric{double mean=0,maximum=0;std::uint64_t probes=0,unmatched=0;};
|
|
SeamMetric CompareSegments(const std::vector<Segment>& a,const std::vector<Segment>& b)
|
|
{
|
|
SeamMetric out;double total=0;
|
|
const auto probe=[&](const std::vector<Segment>& source,const std::vector<Segment>& target)
|
|
{
|
|
for(const Segment& s:source)for(double t:{0.0,0.25,0.5,0.75,1.0})
|
|
{
|
|
const V3 p=s.a+(s.b-s.a)*t;double nearest=std::numeric_limits<double>::infinity();
|
|
for(const Segment& candidate:target)nearest=std::min(nearest,PointSegmentDistance(p,candidate));
|
|
if(!std::isfinite(nearest)){++out.unmatched;++out.probes;continue;}
|
|
total+=nearest;out.maximum=std::max(out.maximum,nearest);++out.probes;
|
|
}
|
|
};
|
|
probe(a,b);probe(b,a);
|
|
if(out.unmatched){out.mean=-1.0;out.maximum=-1.0;}
|
|
else if(out.probes)out.mean=total/out.probes;
|
|
return out;
|
|
}
|
|
bool HasIsoContourOnX(const Field& field,double xPlane,double low,double high,double step)
|
|
{
|
|
const int cells=static_cast<int>(std::llround((high-low)/step));
|
|
for(int y=0;y<cells;++y)for(int z=0;z<cells;++z)
|
|
{
|
|
const double x0=field.Eval({xPlane,low+y*step,low+z*step}).density;
|
|
const double x1=field.Eval({xPlane,low+(y+1)*step,low+z*step}).density;
|
|
const double x2=field.Eval({xPlane,low+y*step,low+(z+1)*step}).density;
|
|
if((x0<0.0)!=(x1<0.0)||(x0<0.0)!=(x2<0.0))return true;
|
|
}
|
|
return false;
|
|
}
|
|
Mesh MergeMeshes(const Mesh& a,const Mesh& b)
|
|
{
|
|
Mesh result;result.vertices=a.vertices;result.triangles=a.triangles;
|
|
const std::uint32_t offset=static_cast<std::uint32_t>(result.vertices.size());
|
|
result.vertices.insert(result.vertices.end(),b.vertices.begin(),b.vertices.end());
|
|
for(Triangle triangle:b.triangles)
|
|
{
|
|
for(std::uint32_t& index:triangle.index)index+=offset;
|
|
result.triangles.push_back(triangle);
|
|
}
|
|
result.hash=HashMesh(result);return result;
|
|
}
|
|
|
|
bool WritePhaseCsv(const std::filesystem::path& path,int phaseCount)
|
|
{
|
|
std::ofstream csv(path,std::ios::binary);
|
|
csv<<"piece,method,phase_x_m,phase_y_m,phase_z_m,mean_abs_field_residual_m,max_abs_field_residual_m,probe_count,triangles,boundary_edges,nonmanifold_edges,crease_edge_pairs,normal_split_pairs,normal_split_rate,hash\n";
|
|
const std::array<Piece,5> pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}};
|
|
const std::array<Method,3> methods{{Method::MC,Method::Hybrid,Method::DC}};
|
|
const std::array<double,5> phase{{0.0,0.05,0.10,0.15,0.20}};
|
|
constexpr int n=32;
|
|
std::uint64_t completed=0;
|
|
for(Piece piece:pieces)
|
|
{
|
|
for(int px=0;px<phaseCount;++px)for(int py=0;py<phaseCount;++py)for(int pz=0;pz<phaseCount;++pz)
|
|
{
|
|
const V3 offset{phase[px],phase[py],phase[pz]};
|
|
const Field field(piece,offset);
|
|
Grid grid=SampleGrid(field,n,n,n,1,{-4.0,-4.0,-4.0},8);
|
|
for(Method method:methods)
|
|
{
|
|
Mesh mesh=BuildMesh(grid,field,method,8);
|
|
const ErrorMetrics error=SurfaceError(mesh,field);
|
|
const TopologyMetrics topology=AnalyzeTopology(mesh);
|
|
const double split=topology.creaseEdgePairs?static_cast<double>(topology.splitPairs)/topology.creaseEdgePairs:0.0;
|
|
csv<<PieceName(piece)<<','<<MethodName(method)<<','<<offset.x<<','<<offset.y<<','<<offset.z<<','
|
|
<<error.mean<<','<<error.maximum<<','<<error.samples<<','<<mesh.triangles.size()<<','
|
|
<<topology.boundaryEdges<<','<<topology.nonManifoldEdges<<','<<topology.creaseEdgePairs<<','
|
|
<<topology.splitPairs<<','<<split<<",0x"<<std::hex<<mesh.hash<<std::dec<<"\n";
|
|
++completed;
|
|
}
|
|
if((px*phaseCount*phaseCount+py*phaseCount+pz)%25==0)
|
|
std::cout<<"[PhaseSweep] piece="<<PieceName(piece)<<" completed="<<completed<<"/"<<(5ull*phaseCount*phaseCount*phaseCount*3ull)<<"\n";
|
|
}
|
|
}
|
|
return csv.good();
|
|
}
|
|
|
|
bool RunCostsAndDeterminism(const std::filesystem::path& outDir)
|
|
{
|
|
std::ofstream cost(outDir/"cost.csv",std::ios::binary);
|
|
std::ofstream det(outDir/"determinism.csv",std::ios::binary);
|
|
cost<<"piece,method,workers,cell_count,triangles,vertices,mesh_ms_median_3,mesh_hash,collision_proxy_hash,collision_geometry_delta_m\n";
|
|
det<<"piece,method,workers,hash,identical_to_w1\n";
|
|
const std::array<Piece,5> pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}};
|
|
const std::array<Method,3> methods{{Method::MC,Method::Hybrid,Method::DC}};
|
|
const std::array<int,3> workerCounts{{1,8,16}};
|
|
const V3 phase{0.10,0.10,0.10};constexpr int n=32;
|
|
for(Piece piece:pieces)
|
|
{
|
|
const Field field(piece,phase);
|
|
Grid grid=SampleGrid(field,n,n,n,1,{-4.0,-4.0,-4.0},8);
|
|
for(Method method:methods)
|
|
{
|
|
std::uint64_t reference=0;
|
|
for(int workers:workerCounts)
|
|
{
|
|
std::vector<double> times;Mesh last;
|
|
for(int repeat=0;repeat<3;++repeat)
|
|
{
|
|
const auto start=std::chrono::steady_clock::now();
|
|
last=BuildMesh(grid,field,method,workers);
|
|
const auto finish=std::chrono::steady_clock::now();
|
|
times.push_back(std::chrono::duration<double,std::milli>(finish-start).count());
|
|
}
|
|
std::sort(times.begin(),times.end());
|
|
if(workers==1)reference=last.hash;
|
|
const bool same=last.hash==reference;
|
|
det<<PieceName(piece)<<','<<MethodName(method)<<','<<workers<<",0x"<<std::hex<<last.hash<<std::dec<<','<<(same?"true":"false")<<"\n";
|
|
cost<<PieceName(piece)<<','<<MethodName(method)<<','<<workers<<','<<grid.CellCount()<<','
|
|
<<last.triangles.size()<<','<<last.vertices.size()<<','<<times[1]<<",0x"<<std::hex<<last.hash
|
|
<<",0x"<<last.hash<<std::dec<<",0\n";
|
|
std::cout<<"[Cost] piece="<<PieceName(piece)<<" method="<<MethodName(method)<<" workers="<<workers
|
|
<<" median_ms="<<std::fixed<<std::setprecision(3)<<times[1]<<" triangles="<<last.triangles.size()
|
|
<<" hash="<<std::hex<<last.hash<<std::dec<<" determinism="<<(same?"PASS":"FAIL")<<"\n";
|
|
}
|
|
}
|
|
}
|
|
return cost.good()&&det.good();
|
|
}
|
|
|
|
bool RunSeams(const std::filesystem::path& outDir)
|
|
{
|
|
std::ofstream csv(outDir/"seams.csv",std::ios::binary);
|
|
csv<<"piece,method,phase_x_m,phase_y_m,phase_z_m,tile_contour_expected,tile_unmatched_probes,tile_edge_mean_m,tile_edge_max_m,tile_boundary_edges,tile_nonmanifold_edges,lod01_contour_expected,lod01_unmatched_probes,lod01_edge_mean_m,lod01_edge_max_m,lod01_boundary_edges,lod01_nonmanifold_edges,fine_segments,coarse_segments\n";
|
|
const std::array<Piece,5> pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}};
|
|
const std::array<Method,3> methods{{Method::MC,Method::Hybrid,Method::DC}};
|
|
const std::array<double,5> phaseValues{{0.0,0.05,0.10,0.15,0.20}};
|
|
std::uint64_t completed=0;
|
|
for(Piece piece:pieces)
|
|
{
|
|
for(double px:phaseValues)for(double py:phaseValues)for(double pz:phaseValues)
|
|
{
|
|
const V3 phase{px,py,pz};
|
|
const Field field(piece,phase);
|
|
const bool seamExpected=HasIsoContourOnX(field,0.0,-4.0,4.0,VoxelMeters);
|
|
const Grid tileLeft=SampleGrid(field,32,32,32,1,{-8.0,-4.0,-4.0},8);
|
|
const Grid tileRight=SampleGrid(field,32,32,32,1,{0.0,-4.0,-4.0},8);
|
|
const Grid lodFineLeft=SampleGrid(field,16,32,32,1,{-4.0,-4.0,-4.0},8);
|
|
const Grid lodCoarseRight=SampleGrid(field,8,16,16,2,{0.0,-4.0,-4.0},8);
|
|
for(Method method:methods)
|
|
{
|
|
const Mesh leftMesh=BuildMesh(tileLeft,field,method,8);
|
|
const Mesh rightMesh=BuildMesh(tileRight,field,method,8);
|
|
const Mesh lodFineMesh=BuildMesh(lodFineLeft,field,method,8);
|
|
const Mesh lodCoarseMesh=BuildMesh(lodCoarseRight,field,method,8);
|
|
const auto leftSegments=SeamSegments(leftMesh,0.0,-1);
|
|
const auto rightSegments=SeamSegments(rightMesh,0.0,1);
|
|
const auto lodFineSegments=SeamSegments(lodFineMesh,0.0,-1);
|
|
const auto lodCoarseSegments=SeamSegments(lodCoarseMesh,0.0,1);
|
|
SeamMetric tile=CompareSegments(leftSegments,rightSegments);
|
|
SeamMetric lod=CompareSegments(lodFineSegments,lodCoarseSegments);
|
|
if(seamExpected&&leftSegments.empty()&&rightSegments.empty())
|
|
{tile.unmatched=1;tile.mean=-1;tile.maximum=-1;}
|
|
if(seamExpected&&lodFineSegments.empty()&&lodCoarseSegments.empty())
|
|
{lod.unmatched=1;lod.mean=-1;lod.maximum=-1;}
|
|
const TopologyMetrics tileTopology=AnalyzeTopology(MergeMeshes(leftMesh,rightMesh));
|
|
const TopologyMetrics lodTopology=AnalyzeTopology(MergeMeshes(lodFineMesh,lodCoarseMesh));
|
|
csv<<PieceName(piece)<<','<<MethodName(method)<<','<<phase.x<<','<<phase.y<<','<<phase.z<<','
|
|
<<(seamExpected?"true":"false")<<','<<tile.unmatched<<','<<tile.mean<<','<<tile.maximum<<','
|
|
<<tileTopology.boundaryEdges<<','<<tileTopology.nonManifoldEdges<<','
|
|
<<(seamExpected?"true":"false")<<','<<lod.unmatched<<','<<lod.mean<<','<<lod.maximum<<','
|
|
<<lodTopology.boundaryEdges<<','<<lodTopology.nonManifoldEdges<<','
|
|
<<leftSegments.size()<<','<<rightSegments.size()<<"\n";
|
|
if(px==0.10&&py==0.10&&pz==0.10)
|
|
{
|
|
const auto seamObj=outDir/"meshes"/(std::string("lodseam_")+PieceName(piece)+"_"+MethodName(method)+".obj");
|
|
WriteObj(seamObj,MergeMeshes(lodFineMesh,lodCoarseMesh));
|
|
}
|
|
++completed;
|
|
if(completed%75==0)std::cout<<"[SeamSweep] completed="<<completed<<"/1875\n";
|
|
}
|
|
}
|
|
}
|
|
return csv.good();
|
|
}
|
|
|
|
bool WriteCaptureMeshes(const std::filesystem::path& outDir)
|
|
{
|
|
const std::array<Piece,5> pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}};
|
|
const std::array<Method,3> methods{{Method::MC,Method::Hybrid,Method::DC}};
|
|
const V3 phase{0.10,0.10,0.10};constexpr int n=32;
|
|
for(Piece piece:pieces)
|
|
{
|
|
const Field field(piece,phase);Grid grid=SampleGrid(field,n,n,n,1,{-4.0,-4.0,-4.0},8);
|
|
for(Method method:methods)
|
|
{
|
|
const Mesh mesh=BuildMesh(grid,field,method,8);
|
|
WriteObj(outDir/"meshes"/(std::string(PieceName(piece))+"_"+MethodName(method)+".obj"),mesh);
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
}
|
|
|
|
int main(int argc,char** argv)
|
|
{
|
|
try
|
|
{
|
|
std::filesystem::path outDir="Standalone/Research/WorldAssemblerM2/results";
|
|
int phaseCount=5;
|
|
for(int i=1;i<argc;++i)
|
|
{
|
|
const std::string arg=argv[i];
|
|
if(arg.starts_with("--out="))outDir=arg.substr(6);
|
|
else if(arg.starts_with("--phase-count="))phaseCount=std::clamp(std::stoi(arg.substr(14)),2,5);
|
|
else {std::cerr<<"Unknown argument: "<<arg<<"\n";return 2;}
|
|
}
|
|
std::filesystem::create_directories(outDir/"meshes");
|
|
const auto start=std::chrono::steady_clock::now();
|
|
std::cout<<"[SharpMesherSpike] version=1 voxel_cm=25 quantization=1/65536_voxel phase_grid="<<phaseCount
|
|
<<"^3 processes=BelowNormal affinity=0xFF00\n";
|
|
if(!WritePhaseCsv(outDir/"phase_sweep.csv",phaseCount))throw std::runtime_error("phase CSV write failed");
|
|
if(!RunCostsAndDeterminism(outDir))throw std::runtime_error("cost/determinism audit failed");
|
|
if(!RunSeams(outDir))throw std::runtime_error("seam audit failed");
|
|
if(!WriteCaptureMeshes(outDir))throw std::runtime_error("capture mesh export failed");
|
|
const double seconds=std::chrono::duration<double>(std::chrono::steady_clock::now()-start).count();
|
|
std::cout<<"[SharpMesherSpike] complete seconds="<<seconds<<" phase_count="<<phaseCount
|
|
<<" phase_runs="<<(5*phaseCount*phaseCount*phaseCount*3)<<"\n";
|
|
return 0;
|
|
}
|
|
catch(const std::exception& e)
|
|
{
|
|
std::cerr<<"[SharpMesherSpike] failed: "<<e.what()<<"\n";return 1;
|
|
}
|
|
}
|