Fr0zkaandClaude Opus 5.5
a3b7d2015e
Shared halo + deterministic border-quad ownership + T-junction conformance: 625/625 same-LOD for all three. LOD 2:1 prototype (collar zipper) fails MC 5 / hybrid 38 / DC 40 cases (corners, overhangs). Hybrid is bit-identical to MC on a sampled native SurfaceWorld cliff; DC rougher. Standalone only. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MmSoF64BbRFmk9C3p1VsBJ
2227 lines
119 KiB
C++
2227 lines
119 KiB
C++
#include "CoreMinimal.h"
|
|
#include "MarchingCubesTables.h"
|
|
#include "VoxelForgeCore/VoxelForgeCore.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 <set>
|
|
#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 LodContourErrorBudgetMeters = 0.50; // One LOD1 cell (50 cm) physical mismatch.
|
|
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, NativeSurface=5 };
|
|
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";
|
|
case Piece::Organic:return "organic_cave_chamber";
|
|
default:return "native_surface_cliff";
|
|
}
|
|
}
|
|
struct Field
|
|
{
|
|
Piece piece=Piece::Corner;
|
|
V3 phaseShift{};
|
|
const VoxelForgeCore::Heightfield2D* nativeSurface=nullptr;
|
|
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:case Piece::NativeSurface:break;
|
|
}
|
|
}
|
|
explicit Field(const VoxelForgeCore::Heightfield2D& surface,V3 phase={})
|
|
:piece(Piece::NativeSurface),phaseShift(phase),nativeSurface(&surface){}
|
|
|
|
FieldSample Eval(V3 p) const
|
|
{
|
|
p=p+phaseShift;
|
|
if(nativeSurface)
|
|
{
|
|
constexpr double SubvoxelsPerMeter=1024.0;
|
|
constexpr std::int64_t GradientDeltaSubvoxels=64;
|
|
const auto x=static_cast<std::int64_t>(std::llround(p.x*SubvoxelsPerMeter));
|
|
const auto y=static_cast<std::int64_t>(std::llround(p.y*SubvoxelsPerMeter));
|
|
const auto z=static_cast<std::int64_t>(std::llround(p.z*SubvoxelsPerMeter));
|
|
const auto center=nativeSurface->point(x,y,z);
|
|
FieldSample result;result.density=static_cast<double>(center.densitySubvoxels)/SubvoxelsPerMeter;
|
|
const auto xp=nativeSurface->point(x+GradientDeltaSubvoxels,y,z).densitySubvoxels;
|
|
const auto xm=nativeSurface->point(x-GradientDeltaSubvoxels,y,z).densitySubvoxels;
|
|
const auto yp=nativeSurface->point(x,y+GradientDeltaSubvoxels,z).densitySubvoxels;
|
|
const auto ym=nativeSurface->point(x,y-GradientDeltaSubvoxels,z).densitySubvoxels;
|
|
const auto zp=nativeSurface->point(x,y,z+GradientDeltaSubvoxels).densitySubvoxels;
|
|
const auto zm=nativeSurface->point(x,y,z-GradientDeltaSubvoxels).densitySubvoxels;
|
|
result.gradient=Normalize({static_cast<double>(xp-xm),static_cast<double>(yp-ym),static_cast<double>(zp-zm)});
|
|
return result;
|
|
}
|
|
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;
|
|
std::array<int,3> coreMin{{0,0,0}},coreMax{{0,0,0}};
|
|
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; }
|
|
bool IsCoreCell(int x,int y,int z) const
|
|
{ return x>=coreMin[0]&&y>=coreMin[1]&&z>=coreMin[2]&&
|
|
x<coreMax[0]&&y<coreMax[1]&&z<coreMax[2]; }
|
|
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,int halo=0)
|
|
{
|
|
Grid grid;grid.step=step;
|
|
grid.coreMin={halo,halo,halo};grid.coreMax={halo+nx,halo+ny,halo+nz};
|
|
grid.nx=nx+2*halo;grid.ny=ny+2*halo;grid.nz=nz+2*halo;
|
|
grid.origin=origin-V3{halo*step*VoxelMeters,halo*step*VoxelMeters,halo*step*VoxelMeters};
|
|
const std::size_t count=static_cast<std::size_t>(grid.nx+1)*(grid.ny+1)*(grid.nz+1);
|
|
grid.density.resize(count);
|
|
ParallelFor(count,workers,[&](std::size_t index)
|
|
{
|
|
const int x=static_cast<int>(index%(grid.nx+1));
|
|
const int y=static_cast<int>((index/(grid.nx+1))%(grid.ny+1));
|
|
const int z=static_cast<int>(index/((grid.nx+1)*(grid.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,bool enforceQuadOwnership=true)
|
|
{
|
|
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)
|
|
{
|
|
if(!grid.IsCoreCell(x,y,z))continue;
|
|
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]];
|
|
}
|
|
// Own each dual quad by the lexicographically lowest incident cell. A one-cell
|
|
// sampled halo supplies the neighbour's QEF vertex, so independently built tiles
|
|
// emit the same complete border quad exactly once.
|
|
const int ownerX=std::min({around[0][0],around[1][0],around[2][0],around[3][0]});
|
|
const int ownerY=std::min({around[0][1],around[1][1],around[2][1],around[3][1]});
|
|
const int ownerZ=std::min({around[0][2],around[1][2],around[2][2],around[3][2]});
|
|
if(enforceQuadOwnership&&!grid.IsCoreCell(ownerX,ownerY,ownerZ))continue;
|
|
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 BoundaryLoop
|
|
{
|
|
std::vector<GeomKey> keys;
|
|
std::vector<V3> points;
|
|
V3 center{};
|
|
double perimeter=0.0;
|
|
};
|
|
struct BoundaryLoopSet
|
|
{
|
|
std::vector<BoundaryLoop> loops;
|
|
std::uint64_t boundaryEdges=0,badVertices=0;
|
|
};
|
|
std::pair<GeomKey,GeomKey> OrderedEdge(GeomKey a,GeomKey b)
|
|
{if(b<a)std::swap(a,b);return {a,b};}
|
|
using EdgeKey=std::pair<GeomKey,GeomKey>;
|
|
BoundaryLoopSet ExtractBoundaryLoops(const Mesh& mesh)
|
|
{
|
|
BoundaryLoopSet result;
|
|
const auto edges=GeomEdges(mesh);
|
|
std::map<GeomKey,std::vector<GeomKey>> adjacency;
|
|
for(const auto& [key,info]:edges)if(info.count==1)
|
|
{
|
|
++result.boundaryEdges;
|
|
adjacency[key.first].push_back(key.second);
|
|
adjacency[key.second].push_back(key.first);
|
|
}
|
|
for(const auto& [key,neighbors]:adjacency)if(neighbors.size()!=2)++result.badVertices;
|
|
std::set<std::pair<GeomKey,GeomKey>> visited;
|
|
for(const auto& [start,startNeighbors]:adjacency)
|
|
{
|
|
for(const GeomKey first:startNeighbors)
|
|
{
|
|
if(visited.contains(OrderedEdge(start,first)))continue;
|
|
BoundaryLoop loop;
|
|
GeomKey previous=start,current=first;
|
|
visited.insert(OrderedEdge(previous,current));
|
|
loop.keys.push_back(start);
|
|
bool closed=false;
|
|
for(std::size_t guard=0;guard<=result.boundaryEdges;++guard)
|
|
{
|
|
loop.keys.push_back(current);
|
|
const auto found=adjacency.find(current);
|
|
if(found==adjacency.end()||found->second.size()!=2)break;
|
|
const GeomKey next=found->second[0]==previous?found->second[1]:found->second[0];
|
|
if(next==start){closed=true;break;}
|
|
const auto edge=OrderedEdge(current,next);
|
|
if(visited.contains(edge))break;
|
|
visited.insert(edge);previous=current;current=next;
|
|
}
|
|
if(!closed||loop.keys.size()<3)continue;
|
|
BoundaryLoop clean;
|
|
clean.keys=std::move(loop.keys);
|
|
for(const GeomKey& key:clean.keys)
|
|
{
|
|
const V3 p{key.x*PositionQuantum,key.y*PositionQuantum,key.z*PositionQuantum};
|
|
clean.points.push_back(p);clean.center+=p;
|
|
}
|
|
clean.center=clean.center/static_cast<double>(clean.points.size());
|
|
for(std::size_t i=0;i<clean.points.size();++i)
|
|
clean.perimeter+=Length(clean.points[(i+1)%clean.points.size()]-clean.points[i]);
|
|
result.loops.push_back(std::move(clean));
|
|
}
|
|
}
|
|
return result;
|
|
}
|
|
|
|
struct TransitionAudit
|
|
{
|
|
std::uint64_t pairedLoops=0,unmatchedLoops=0,cappedLoops=0,transitionTriangles=0;
|
|
double contourMean=0.0,contourMaximum=0.0;
|
|
};
|
|
V3 PointAlongLoop(const BoundaryLoop& loop,std::size_t start,int direction,double fraction)
|
|
{
|
|
const std::size_t count=loop.points.size();
|
|
if(!count)return {};
|
|
double remaining=std::clamp(fraction,0.0,1.0)*loop.perimeter;
|
|
for(std::size_t walked=0;walked<count;++walked)
|
|
{
|
|
const std::size_t current=direction>0?(start+walked)%count:(start+count-(walked%count))%count;
|
|
const std::size_t next=direction>0?(current+1)%count:(current+count-1)%count;
|
|
const V3 a=loop.points[current],b=loop.points[next];
|
|
const double edgeLength=Length(b-a);
|
|
if(remaining<=edgeLength||walked+1==count)
|
|
return edgeLength>0.0?a+(b-a)*(remaining/edgeLength):a;
|
|
remaining-=edgeLength;
|
|
}
|
|
return loop.points[start];
|
|
}
|
|
int BestLoopDirection(const BoundaryLoop& a,const BoundaryLoop& b,std::size_t start)
|
|
{
|
|
double scores[2]{};
|
|
constexpr int samples=24;
|
|
for(int dirIndex=0;dirIndex<2;++dirIndex)
|
|
{
|
|
const int direction=dirIndex==0?1:-1;
|
|
for(int i=0;i<samples;++i)
|
|
{
|
|
const double f=static_cast<double>(i)/samples;
|
|
scores[dirIndex]+=Dot(PointAlongLoop(a,0,1,f)-PointAlongLoop(b,start,direction,f),
|
|
PointAlongLoop(a,0,1,f)-PointAlongLoop(b,start,direction,f));
|
|
}
|
|
}
|
|
return scores[0]<=scores[1]?1:-1;
|
|
}
|
|
TransitionAudit StitchBoundaryLoops(Mesh& mesh,const BoundaryLoopSet& fine,
|
|
const BoundaryLoopSet& coarse,const Field& field)
|
|
{
|
|
TransitionAudit result;
|
|
const std::size_t patchStart=mesh.triangles.size();
|
|
result.unmatchedLoops=fine.badVertices+coarse.badVertices;
|
|
std::map<GeomKey,std::uint32_t> vertexByKey;
|
|
for(std::uint32_t i=0;i<mesh.vertices.size();++i)vertexByKey.try_emplace(Key(mesh.vertices[i].position),i);
|
|
const auto loopEdges=[](const BoundaryLoopSet& loops)
|
|
{
|
|
std::set<EdgeKey> edges;
|
|
for(const BoundaryLoop& loop:loops.loops)
|
|
for(std::size_t i=0;i<loop.keys.size();++i)
|
|
edges.insert(OrderedEdge(loop.keys[i],loop.keys[(i+1)%loop.keys.size()]));
|
|
return edges;
|
|
};
|
|
const std::set<EdgeKey> fineEdges=loopEdges(fine),coarseEdges=loopEdges(coarse);
|
|
std::set<EdgeKey> sharedEdges;
|
|
std::set_intersection(fineEdges.begin(),fineEdges.end(),coarseEdges.begin(),coarseEdges.end(),
|
|
std::inserter(sharedEdges,sharedEdges.end()));
|
|
constexpr double collarOffset=0.001;
|
|
if(!sharedEdges.empty()&&fine.loops.size()==coarse.loops.size())
|
|
{
|
|
struct Chain{std::vector<GeomKey> keys;};
|
|
const auto buildChains=[&](const BoundaryLoopSet& source)
|
|
{
|
|
std::vector<Chain> chains;
|
|
for(const BoundaryLoop& loop:source.loops)
|
|
{
|
|
const std::size_t n=loop.keys.size();
|
|
std::vector<bool> isShared(n,false);
|
|
for(std::size_t i=0;i<n;++i)
|
|
isShared[i]=sharedEdges.contains(OrderedEdge(loop.keys[i],loop.keys[(i+1)%n]));
|
|
if(std::none_of(isShared.begin(),isShared.end(),[](bool value){return value;}))
|
|
{chains.push_back({loop.keys});continue;}
|
|
for(std::size_t edge=0;edge<n;++edge)if(!isShared[edge]&&isShared[(edge+n-1)%n])
|
|
{
|
|
Chain chain;std::size_t current=edge;chain.keys.push_back(loop.keys[current]);
|
|
for(std::size_t guard=0;guard<n&&!isShared[current];++guard)
|
|
{
|
|
current=(current+1)%n;
|
|
chain.keys.push_back(loop.keys[current]);
|
|
}
|
|
if(chain.keys.size()>1)chains.push_back(std::move(chain));
|
|
}
|
|
}
|
|
return chains;
|
|
};
|
|
const std::vector<Chain> fineChains=buildChains(fine),coarseChains=buildChains(coarse);
|
|
std::vector<bool> coarseUsed(coarseChains.size(),false);
|
|
const std::size_t transitionStart=mesh.triangles.size();
|
|
const auto offsetChain=[&](const Chain& chain,double offset)
|
|
{
|
|
std::vector<std::uint32_t> original,shifted;
|
|
for(const GeomKey& key:chain.keys)
|
|
{
|
|
const auto found=vertexByKey.find(key);
|
|
if(found==vertexByKey.end())return std::pair{original,shifted};
|
|
original.push_back(found->second);
|
|
Vertex vertex=mesh.vertices[found->second];vertex.position.x+=offset;
|
|
const GeomKey shiftedKey=Key(vertex.position);
|
|
const auto existing=vertexByKey.find(shiftedKey);
|
|
if(existing!=vertexByKey.end())shifted.push_back(existing->second);
|
|
else
|
|
{
|
|
const auto index=static_cast<std::uint32_t>(mesh.vertices.size());
|
|
mesh.vertices.push_back(vertex);vertexByKey.emplace(shiftedKey,index);shifted.push_back(index);
|
|
}
|
|
}
|
|
for(std::size_t i=0;i+1<original.size();++i)
|
|
{
|
|
const auto a=original[i],b=original[i+1],sa=shifted[i],sb=shifted[i+1];
|
|
const V3 center=(mesh.vertices[a].position+mesh.vertices[b].position+
|
|
mesh.vertices[sb].position+mesh.vertices[sa].position)/4.0;
|
|
const FieldSample sample=field.Eval(center);
|
|
AddTriangle(mesh,a,b,sb,sample.primaryFeature(),sample.gradient);
|
|
AddTriangle(mesh,a,sb,sa,sample.primaryFeature(),sample.gradient);
|
|
}
|
|
return std::pair{original,shifted};
|
|
};
|
|
std::uint64_t matchedChains=0;
|
|
const auto endpointDistance=[](const GeomKey& a,const GeomKey& b)
|
|
{
|
|
const V3 pa{a.x*PositionQuantum,a.y*PositionQuantum,a.z*PositionQuantum};
|
|
const V3 pb{b.x*PositionQuantum,b.y*PositionQuantum,b.z*PositionQuantum};
|
|
return Length(pa-pb);
|
|
};
|
|
for(const Chain& fineChain:fineChains)
|
|
{
|
|
std::size_t match=coarseChains.size();bool reverse=false;
|
|
double bestEndpointScore=std::numeric_limits<double>::infinity();
|
|
for(std::size_t j=0;j<coarseChains.size();++j)if(!coarseUsed[j])
|
|
{
|
|
const Chain& candidate=coarseChains[j];
|
|
const double directA=endpointDistance(fineChain.keys.front(),candidate.keys.front());
|
|
const double directB=endpointDistance(fineChain.keys.back(),candidate.keys.back());
|
|
const double reverseA=endpointDistance(fineChain.keys.front(),candidate.keys.back());
|
|
const double reverseB=endpointDistance(fineChain.keys.back(),candidate.keys.front());
|
|
const double directScore=std::max(directA,directB)+1.0e-6*(directA+directB);
|
|
const double reverseScore=std::max(reverseA,reverseB)+1.0e-6*(reverseA+reverseB);
|
|
const bool candidateReverse=reverseScore<directScore;
|
|
const double candidateScore=std::min(directScore,reverseScore);
|
|
const double maxDistance=candidateReverse?std::max(reverseA,reverseB):std::max(directA,directB);
|
|
if(maxDistance<=LodContourErrorBudgetMeters&&candidateScore<bestEndpointScore)
|
|
{match=j;reverse=candidateReverse;bestEndpointScore=candidateScore;}
|
|
}
|
|
if(match==coarseChains.size()){++result.unmatchedLoops;continue;}
|
|
coarseUsed[match]=true;++matchedChains;
|
|
Chain alignedCoarse=coarseChains[match];
|
|
if(reverse)std::reverse(alignedCoarse.keys.begin(),alignedCoarse.keys.end());
|
|
const auto [fineOriginal,fineOffset]=offsetChain(fineChain,-collarOffset);
|
|
const auto [coarseOriginal,coarseOffset]=offsetChain(alignedCoarse,collarOffset);
|
|
if(fineOffset.size()<2||coarseOffset.size()<2||fineOriginal.size()!=fineOffset.size()||
|
|
coarseOriginal.size()!=coarseOffset.size())
|
|
{++result.unmatchedLoops;continue;}
|
|
for(const auto endpoint:{std::size_t{0},fineOffset.size()-1})
|
|
{
|
|
const auto fa=fineOffset[endpoint];
|
|
const auto ca=coarseOffset[endpoint==0?0:coarseOffset.size()-1];
|
|
const auto originalFine=fineOriginal[endpoint];
|
|
const auto originalCoarse=coarseOriginal[endpoint==0?0:coarseOriginal.size()-1];
|
|
std::vector<std::uint32_t> capRing{originalFine,fa,ca,originalCoarse};
|
|
capRing.erase(std::unique(capRing.begin(),capRing.end(),[&](std::uint32_t a,std::uint32_t b)
|
|
{return Key(mesh.vertices[a].position)==Key(mesh.vertices[b].position);}),capRing.end());
|
|
if(capRing.size()>2&&Key(mesh.vertices[capRing.front()].position)==Key(mesh.vertices[capRing.back()].position))capRing.pop_back();
|
|
V3 capCenter{};for(const std::uint32_t index:capRing)capCenter+=mesh.vertices[index].position;
|
|
capCenter=capCenter/static_cast<double>(capRing.size())+V3{0.0,0.001,0.001};
|
|
Vertex capVertex;capVertex.position=capCenter;
|
|
capVertex.normal=Normalize(mesh.vertices[fa].normal+mesh.vertices[ca].normal);
|
|
const GeomKey capKey=Key(capVertex.position);
|
|
auto capFound=vertexByKey.find(capKey);
|
|
std::uint32_t cap=0;
|
|
if(capFound!=vertexByKey.end())cap=capFound->second;
|
|
else
|
|
{
|
|
cap=static_cast<std::uint32_t>(mesh.vertices.size());
|
|
mesh.vertices.push_back(capVertex);vertexByKey.emplace(capKey,cap);
|
|
}
|
|
for(std::size_t i=0;i<capRing.size();++i)
|
|
{
|
|
const std::uint32_t a=capRing[i],b=capRing[(i+1)%capRing.size()];
|
|
const V3 center=(mesh.vertices[cap].position+mesh.vertices[a].position+
|
|
mesh.vertices[b].position)/3.0;
|
|
const FieldSample sample=field.Eval(center);
|
|
AddTriangle(mesh,cap,a,b,sample.primaryFeature(),sample.gradient);
|
|
}
|
|
}
|
|
const std::size_t na=fineOffset.size(),nb=coarseOffset.size();
|
|
std::vector<double> cumulativeA(na,0.0),cumulativeB(nb,0.0);
|
|
for(std::size_t i=1;i<na;++i)cumulativeA[i]=cumulativeA[i-1]+Length(mesh.vertices[fineOffset[i]].position-mesh.vertices[fineOffset[i-1]].position);
|
|
for(std::size_t i=1;i<nb;++i)cumulativeB[i]=cumulativeB[i-1]+Length(mesh.vertices[coarseOffset[i]].position-mesh.vertices[coarseOffset[i-1]].position);
|
|
const double perimA=cumulativeA.back(),perimB=cumulativeB.back();
|
|
if(perimA<=0.0||perimB<=0.0){++result.unmatchedLoops;continue;}
|
|
std::size_t i=0,j=0;
|
|
while(i+1<na||j+1<nb)
|
|
{
|
|
const double nextA=i+1<na?cumulativeA[i+1]/perimA:std::numeric_limits<double>::infinity();
|
|
const double nextB=j+1<nb?cumulativeB[j+1]/perimB:std::numeric_limits<double>::infinity();
|
|
if(nextA<=nextB)
|
|
{
|
|
const V3 center=(mesh.vertices[fineOffset[i]].position+mesh.vertices[fineOffset[i+1]].position+
|
|
mesh.vertices[coarseOffset[j]].position)/3.0;
|
|
const FieldSample sample=field.Eval(center);
|
|
AddTriangle(mesh,fineOffset[i],fineOffset[i+1],coarseOffset[j],sample.primaryFeature(),sample.gradient);++i;
|
|
}
|
|
else
|
|
{
|
|
const V3 center=(mesh.vertices[fineOffset[i]].position+mesh.vertices[coarseOffset[j+1]].position+
|
|
mesh.vertices[coarseOffset[j]].position)/3.0;
|
|
const FieldSample sample=field.Eval(center);
|
|
AddTriangle(mesh,fineOffset[i],coarseOffset[j+1],coarseOffset[j],sample.primaryFeature(),sample.gradient);++j;
|
|
}
|
|
}
|
|
}
|
|
for(bool used:coarseUsed)if(!used)++result.unmatchedLoops;
|
|
result.pairedLoops=matchedChains;
|
|
result.transitionTriangles=mesh.triangles.size()-transitionStart;
|
|
return result;
|
|
}
|
|
// Give both sides a one-millimetre transition skin before zipping. This avoids degenerate
|
|
// triangles when otherwise independent contours share collinear runs.
|
|
const auto offsetLoops=[&](const BoundaryLoopSet& source,double offset)
|
|
{
|
|
std::vector<BoundaryLoop> output;output.reserve(source.loops.size());
|
|
for(const BoundaryLoop& input:source.loops)
|
|
{
|
|
BoundaryLoop loop;bool valid=true;
|
|
for(const GeomKey& key:input.keys)
|
|
{
|
|
const auto found=vertexByKey.find(key);
|
|
if(found==vertexByKey.end()){valid=false;break;}
|
|
Vertex vertex=mesh.vertices[found->second];vertex.position.x+=offset;
|
|
const std::uint32_t index=static_cast<std::uint32_t>(mesh.vertices.size());
|
|
mesh.vertices.push_back(vertex);vertexByKey.try_emplace(Key(vertex.position),index);
|
|
loop.keys.push_back(Key(vertex.position));loop.points.push_back(vertex.position);loop.center+=vertex.position;
|
|
}
|
|
if(!valid||loop.points.size()<3){++result.unmatchedLoops;continue;}
|
|
loop.center=loop.center/static_cast<double>(loop.points.size());
|
|
for(std::size_t i=0;i<loop.points.size();++i)
|
|
{
|
|
const std::size_t next=(i+1)%loop.points.size();
|
|
loop.perimeter+=Length(loop.points[next]-loop.points[i]);
|
|
const auto originalA=vertexByKey.find(input.keys[i]);
|
|
const auto originalB=vertexByKey.find(input.keys[next]);
|
|
const auto shiftedA=vertexByKey.find(loop.keys[i]);
|
|
const auto shiftedB=vertexByKey.find(loop.keys[next]);
|
|
if(originalA==vertexByKey.end()||originalB==vertexByKey.end()||
|
|
shiftedA==vertexByKey.end()||shiftedB==vertexByKey.end())
|
|
{++result.unmatchedLoops;continue;}
|
|
const auto a=originalA->second,b=originalB->second,sa=shiftedA->second,sb=shiftedB->second;
|
|
const V3 center=(mesh.vertices[a].position+mesh.vertices[b].position+
|
|
mesh.vertices[sb].position+mesh.vertices[sa].position)/4.0;
|
|
const FieldSample sample=field.Eval(center);
|
|
AddTriangle(mesh,a,b,sb,sample.primaryFeature(),sample.gradient);
|
|
AddTriangle(mesh,a,sb,sa,sample.primaryFeature(),sample.gradient);
|
|
}
|
|
output.push_back(std::move(loop));
|
|
}
|
|
return output;
|
|
};
|
|
const std::vector<BoundaryLoop> fineOffset=offsetLoops(fine,-collarOffset);
|
|
const std::vector<BoundaryLoop> coarseOffset=offsetLoops(coarse,collarOffset);
|
|
const auto capLoop=[&](const BoundaryLoop& loop)
|
|
{
|
|
if(loop.keys.size()<3)return false;
|
|
V3 center{};for(const V3 point:loop.points)center+=point;
|
|
center=center/static_cast<double>(loop.points.size())+V3{0.001,0.001,0.001};
|
|
Vertex centerVertex;centerVertex.position=center;centerVertex.normal=Normalize(field.Eval(center).gradient);
|
|
const GeomKey centerKey=Key(center);
|
|
auto centerIt=vertexByKey.find(centerKey);std::uint32_t centerIndex=0;
|
|
if(centerIt!=vertexByKey.end())centerIndex=centerIt->second;
|
|
else
|
|
{
|
|
centerIndex=static_cast<std::uint32_t>(mesh.vertices.size());
|
|
mesh.vertices.push_back(centerVertex);vertexByKey.emplace(centerKey,centerIndex);
|
|
}
|
|
const std::size_t before=mesh.triangles.size();
|
|
for(std::size_t i=0;i<loop.keys.size();++i)
|
|
{
|
|
const auto a=vertexByKey.find(loop.keys[i]),b=vertexByKey.find(loop.keys[(i+1)%loop.keys.size()]);
|
|
if(a==vertexByKey.end()||b==vertexByKey.end())continue;
|
|
const V3 at=(mesh.vertices[centerIndex].position+mesh.vertices[a->second].position+
|
|
mesh.vertices[b->second].position)/3.0;
|
|
const FieldSample sample=field.Eval(at);
|
|
AddTriangle(mesh,centerIndex,a->second,b->second,sample.primaryFeature(),sample.gradient);
|
|
}
|
|
if(mesh.triangles.size()==before)return false;
|
|
++result.cappedLoops;return true;
|
|
};
|
|
std::vector<bool> used(coarseOffset.size(),false);
|
|
double totalResidual=0.0;std::uint64_t residualSamples=0;
|
|
for(const BoundaryLoop& a:fineOffset)
|
|
{
|
|
std::size_t match=coarseOffset.size();double best=std::numeric_limits<double>::infinity();
|
|
for(std::size_t j=0;j<coarseOffset.size();++j)if(!used[j])
|
|
{
|
|
const double d=Length(a.center-coarseOffset[j].center);
|
|
if(d<best){best=d;match=j;}
|
|
}
|
|
if(match==coarseOffset.size()){if(!capLoop(a))++result.unmatchedLoops;continue;}
|
|
used[match]=true;++result.pairedLoops;
|
|
const BoundaryLoop& b=coarseOffset[match];
|
|
if(a.keys.size()<3||b.keys.size()<3){++result.unmatchedLoops;continue;}
|
|
std::size_t bStart=0;double nearest=std::numeric_limits<double>::infinity();
|
|
for(std::size_t j=0;j<b.points.size();++j)
|
|
{
|
|
const double d=Length(a.points[0]-b.points[j]);
|
|
if(d<nearest){nearest=d;bStart=j;}
|
|
}
|
|
const int bDirection=BestLoopDirection(a,b,bStart);
|
|
std::vector<std::uint32_t> aIndices,bIndices;
|
|
for(const GeomKey& key:a.keys)
|
|
{
|
|
const auto it=vertexByKey.find(key);
|
|
if(it==vertexByKey.end()){++result.unmatchedLoops;break;}
|
|
aIndices.push_back(it->second);
|
|
}
|
|
for(std::size_t j=0;j<b.keys.size();++j)
|
|
{
|
|
const std::size_t bi=bDirection>0?(bStart+j)%b.keys.size():(bStart+b.keys.size()-(j%b.keys.size()))%b.keys.size();
|
|
const auto it=vertexByKey.find(b.keys[bi]);
|
|
if(it==vertexByKey.end()){++result.unmatchedLoops;break;}
|
|
bIndices.push_back(it->second);
|
|
}
|
|
if(aIndices.size()!=a.keys.size()||bIndices.size()!=b.keys.size())continue;
|
|
const std::size_t na=aIndices.size(),nb=bIndices.size();
|
|
std::vector<double> cumulativeA(na+1,0.0),cumulativeB(nb+1,0.0);
|
|
for(std::size_t i=0;i<na;++i)cumulativeA[i+1]=cumulativeA[i]+Length(mesh.vertices[aIndices[(i+1)%na]].position-mesh.vertices[aIndices[i]].position);
|
|
for(std::size_t i=0;i<nb;++i)cumulativeB[i+1]=cumulativeB[i]+Length(mesh.vertices[bIndices[(i+1)%nb]].position-mesh.vertices[bIndices[i]].position);
|
|
const double perimA=cumulativeA.back(),perimB=cumulativeB.back();
|
|
if(perimA<=0.0||perimB<=0.0){++result.unmatchedLoops;continue;}
|
|
std::size_t i=0,j=0;
|
|
while(i<na||j<nb)
|
|
{
|
|
const double nextA=i<na?cumulativeA[i+1]/perimA:std::numeric_limits<double>::infinity();
|
|
const double nextB=j<nb?cumulativeB[j+1]/perimB:std::numeric_limits<double>::infinity();
|
|
if(nextA<=nextB)
|
|
{
|
|
const auto a0=aIndices[i%na],a1=aIndices[(i+1)%na],b0=bIndices[j%nb];
|
|
const V3 center=(mesh.vertices[a0].position+mesh.vertices[a1].position+mesh.vertices[b0].position)/3.0;
|
|
const FieldSample sample=field.Eval(center);
|
|
AddTriangle(mesh,a0,a1,b0,sample.primaryFeature(),sample.gradient);++i;
|
|
}
|
|
else
|
|
{
|
|
const auto a0=aIndices[i%na],b1=bIndices[(j+1)%nb],b0=bIndices[j%nb];
|
|
const V3 center=(mesh.vertices[a0].position+mesh.vertices[b1].position+mesh.vertices[b0].position)/3.0;
|
|
const FieldSample sample=field.Eval(center);
|
|
AddTriangle(mesh,a0,b1,b0,sample.primaryFeature(),sample.gradient);++j;
|
|
}
|
|
}
|
|
const std::size_t sampleCount=std::max(na,nb);
|
|
for(std::size_t s=0;s<sampleCount;++s)
|
|
{
|
|
const double f=static_cast<double>(s)/sampleCount;
|
|
const V3 pa=PointAlongLoop(a,0,1,f),pb=PointAlongLoop(b,bStart,bDirection,f);
|
|
const double d=Length(pa-pb);totalResidual+=d;result.contourMaximum=std::max(result.contourMaximum,d);++residualSamples;
|
|
}
|
|
}
|
|
for(std::size_t j=0;j<used.size();++j)if(!used[j]&&!capLoop(coarseOffset[j]))++result.unmatchedLoops;
|
|
if(residualSamples)result.contourMean=totalResidual/residualSamples;
|
|
result.transitionTriangles=mesh.triangles.size()-patchStart;
|
|
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;
|
|
}
|
|
TopologyMetrics AnalyzeTopologyInXBand(const Mesh& mesh,double plane,double halfWidth)
|
|
{
|
|
TopologyMetrics out;
|
|
for(const auto& [edge,info]:GeomEdges(mesh))
|
|
{
|
|
const double ax=edge.first.x*PositionQuantum,bx=edge.second.x*PositionQuantum;
|
|
if(std::abs(ax-plane)>halfWidth||std::abs(bx-plane)>halfWidth)continue;
|
|
if(info.count==1)++out.boundaryEdges;
|
|
else if(info.count!=2)++out.nonManifoldEdges;
|
|
}
|
|
return out;
|
|
}
|
|
struct SeamCapAudit
|
|
{
|
|
std::uint64_t cappedLoops=0,unmatchedLoops=0,triangles=0;
|
|
};
|
|
SeamCapAudit CapClosedBoundaryLoopsInXBand(Mesh& mesh,const Field& field,double plane,double halfWidth)
|
|
{
|
|
SeamCapAudit result;
|
|
std::map<GeomKey,std::vector<GeomKey>> adjacency;
|
|
std::set<EdgeKey> boundaryEdges;
|
|
std::map<GeomKey,std::uint32_t> vertexByKey;
|
|
for(std::uint32_t i=0;i<mesh.vertices.size();++i)vertexByKey.try_emplace(Key(mesh.vertices[i].position),i);
|
|
for(const auto& [edge,info]:GeomEdges(mesh))if(info.count==1&&
|
|
std::abs(edge.first.x*PositionQuantum-plane)<=halfWidth&&
|
|
std::abs(edge.second.x*PositionQuantum-plane)<=halfWidth)
|
|
{
|
|
boundaryEdges.insert(edge);
|
|
adjacency[edge.first].push_back(edge.second);
|
|
adjacency[edge.second].push_back(edge.first);
|
|
}
|
|
std::set<EdgeKey> visited;
|
|
for(const EdgeKey& firstEdge:boundaryEdges)
|
|
{
|
|
if(visited.contains(firstEdge))continue;
|
|
std::vector<GeomKey> loop{firstEdge.first};
|
|
GeomKey previous=firstEdge.first,current=firstEdge.second;
|
|
visited.insert(firstEdge);
|
|
bool closed=false,valid=true;
|
|
for(std::size_t guard=0;guard<=boundaryEdges.size();++guard)
|
|
{
|
|
if(current==loop.front()){closed=true;break;}
|
|
loop.push_back(current);
|
|
const auto found=adjacency.find(current);
|
|
if(found==adjacency.end()||found->second.size()!=2){valid=false;break;}
|
|
const GeomKey next=found->second[0]==previous?found->second[1]:found->second[0];
|
|
const EdgeKey edge=OrderedEdge(current,next);
|
|
if(visited.contains(edge)){valid=false;break;}
|
|
visited.insert(edge);previous=current;current=next;
|
|
}
|
|
if(!closed||!valid||loop.size()<3)
|
|
{++result.unmatchedLoops;continue;}
|
|
V3 center{};bool haveVertices=true;
|
|
for(const GeomKey& key:loop)
|
|
{
|
|
const auto found=vertexByKey.find(key);
|
|
if(found==vertexByKey.end()){haveVertices=false;break;}
|
|
center+=mesh.vertices[found->second].position;
|
|
}
|
|
if(!haveVertices){++result.unmatchedLoops;continue;}
|
|
center=center/static_cast<double>(loop.size());
|
|
center+=Normalize(field.Eval(center).gradient)*0.0005;
|
|
Vertex centerVertex;centerVertex.position=center;centerVertex.normal=Normalize(field.Eval(center).gradient);
|
|
const GeomKey centerKey=Key(center);
|
|
const std::size_t verticesBefore=mesh.vertices.size();
|
|
auto centerIt=vertexByKey.find(centerKey);std::uint32_t centerIndex=0;
|
|
if(centerIt!=vertexByKey.end())centerIndex=centerIt->second;
|
|
else
|
|
{
|
|
centerIndex=static_cast<std::uint32_t>(mesh.vertices.size());
|
|
mesh.vertices.push_back(centerVertex);vertexByKey.emplace(centerKey,centerIndex);
|
|
}
|
|
const std::size_t before=mesh.triangles.size();
|
|
for(std::size_t i=0;i<loop.size();++i)
|
|
{
|
|
const auto a=vertexByKey.find(loop[i]),b=vertexByKey.find(loop[(i+1)%loop.size()]);
|
|
if(a==vertexByKey.end()||b==vertexByKey.end())continue;
|
|
const V3 samplePoint=(mesh.vertices[centerIndex].position+
|
|
mesh.vertices[a->second].position+mesh.vertices[b->second].position)/3.0;
|
|
const FieldSample sample=field.Eval(samplePoint);
|
|
AddTriangle(mesh,centerIndex,a->second,b->second,sample.primaryFeature(),sample.gradient);
|
|
}
|
|
if(mesh.triangles.size()==before)
|
|
{
|
|
if(mesh.vertices.size()>verticesBefore)
|
|
{mesh.vertices.resize(verticesBefore);vertexByKey.erase(centerKey);}
|
|
++result.unmatchedLoops;continue;
|
|
}
|
|
Mesh capOnly;capOnly.vertices=mesh.vertices;
|
|
capOnly.triangles.assign(mesh.triangles.begin()+static_cast<std::ptrdiff_t>(before),mesh.triangles.end());
|
|
if(SurfaceError(capOnly,field).maximum>LodContourErrorBudgetMeters)
|
|
{
|
|
mesh.triangles.resize(before);
|
|
if(mesh.vertices.size()>verticesBefore)
|
|
{mesh.vertices.resize(verticesBefore);vertexByKey.erase(centerKey);}
|
|
++result.unmatchedLoops;continue;
|
|
}
|
|
++result.cappedLoops;result.triangles+=mesh.triangles.size()-before;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
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;
|
|
}
|
|
Mesh ClipMeshAtX(const Mesh& source,const Field& field,bool keepLeft)
|
|
{
|
|
Mesh result;
|
|
std::map<GeomKey,std::uint32_t> vertexByKey;
|
|
const auto addVertex=[&](Vertex vertex)
|
|
{
|
|
if(std::abs(vertex.position.x)<PositionQuantum)vertex.position.x=0.0;
|
|
const GeomKey key=Key(vertex.position);
|
|
const auto found=vertexByKey.find(key);
|
|
if(found!=vertexByKey.end())return found->second;
|
|
const auto index=static_cast<std::uint32_t>(result.vertices.size());
|
|
result.vertices.push_back(vertex);vertexByKey.emplace(key,index);return index;
|
|
};
|
|
for(const Triangle& triangle:source.triangles)
|
|
{
|
|
std::vector<Vertex> polygon;polygon.reserve(4);
|
|
for(int edge=0;edge<3;++edge)
|
|
{
|
|
const Vertex a=source.vertices[triangle.index[edge]];
|
|
const Vertex b=source.vertices[triangle.index[(edge+1)%3]];
|
|
const bool insideA=keepLeft?a.position.x<=0.0:a.position.x>=0.0;
|
|
const bool insideB=keepLeft?b.position.x<=0.0:b.position.x>=0.0;
|
|
if(insideA)polygon.push_back(a);
|
|
if(insideA!=insideB)
|
|
{
|
|
const double t=std::clamp(-a.position.x/(b.position.x-a.position.x),0.0,1.0);
|
|
Vertex cut;
|
|
cut.position=a.position+(b.position-a.position)*t;cut.position.x=0.0;
|
|
cut.normal=Normalize(a.normal+(b.normal-a.normal)*t);
|
|
polygon.push_back(cut);
|
|
}
|
|
}
|
|
std::vector<std::uint32_t> indices;indices.reserve(polygon.size());
|
|
for(const Vertex& vertex:polygon)
|
|
{
|
|
const auto index=addVertex(vertex);
|
|
if(indices.empty()||indices.back()!=index)indices.push_back(index);
|
|
}
|
|
if(indices.size()>2&&indices.front()==indices.back())indices.pop_back();
|
|
if(indices.size()<3)continue;
|
|
for(std::size_t i=1;i+1<indices.size();++i)
|
|
{
|
|
const V3 center=(result.vertices[indices[0]].position+
|
|
result.vertices[indices[i]].position+result.vertices[indices[i+1]].position)/3.0;
|
|
AddTriangle(result,indices[0],indices[i],indices[i+1],triangle.feature,field.Eval(center).gradient);
|
|
}
|
|
}
|
|
result.hash=HashMesh(result);return result;
|
|
}
|
|
std::vector<Segment> BoundarySegmentsAtX(const Mesh& mesh,double xPlane)
|
|
{
|
|
std::vector<Segment> segments;
|
|
for(const auto& [key,info]:GeomEdges(mesh))
|
|
{
|
|
if(info.count!=1)continue;
|
|
const V3 a{key.first.x*PositionQuantum,key.first.y*PositionQuantum,key.first.z*PositionQuantum};
|
|
const V3 b{key.second.x*PositionQuantum,key.second.y*PositionQuantum,key.second.z*PositionQuantum};
|
|
if(std::abs(a.x-xPlane)<=PositionQuantum&&std::abs(b.x-xPlane)<=PositionQuantum)
|
|
segments.push_back({a,b});
|
|
}
|
|
return segments;
|
|
}
|
|
struct BoundaryRepairAudit
|
|
{
|
|
std::uint64_t triangles=0,unmatchedEdges=0;
|
|
double maximumResidual=0.0;
|
|
};
|
|
std::uint64_t ConformSeamBoundaryEdges(Mesh& a,Mesh& b,double xPlane,const Field& field)
|
|
{
|
|
std::map<GeomKey,V3> uniquePoints;
|
|
for(const Mesh* mesh:{&a,&b})for(const Vertex& vertex:mesh->vertices)
|
|
if(std::abs(vertex.position.x-xPlane)<=PositionQuantum)
|
|
uniquePoints.try_emplace(Key(vertex.position),vertex.position);
|
|
std::vector<V3> points;points.reserve(uniquePoints.size());
|
|
for(const auto& [key,p]:uniquePoints)points.push_back(p);
|
|
const auto conform=[&](Mesh& source)
|
|
{
|
|
Mesh output;output.vertices=source.vertices;
|
|
std::map<GeomKey,std::uint32_t> vertexByKey;
|
|
for(std::uint32_t i=0;i<output.vertices.size();++i)vertexByKey.try_emplace(Key(output.vertices[i].position),i);
|
|
const auto addVertex=[&](Vertex vertex)
|
|
{
|
|
const GeomKey key=Key(vertex.position);
|
|
const auto found=vertexByKey.find(key);
|
|
if(found!=vertexByKey.end())return found->second;
|
|
const auto index=static_cast<std::uint32_t>(output.vertices.size());
|
|
output.vertices.push_back(vertex);vertexByKey.emplace(key,index);return index;
|
|
};
|
|
std::uint64_t splitTriangles=0;
|
|
for(const Triangle& triangle:source.triangles)
|
|
{
|
|
std::vector<std::uint32_t> polygon;polygon.reserve(8);bool split=false;
|
|
for(int edge=0;edge<3;++edge)
|
|
{
|
|
const auto aIndex=triangle.index[edge],bIndex=triangle.index[(edge+1)%3];
|
|
const V3 pa=source.vertices[aIndex].position,pb=source.vertices[bIndex].position;
|
|
polygon.push_back(addVertex(source.vertices[aIndex]));
|
|
if(std::abs(pa.x-xPlane)>PositionQuantum||std::abs(pb.x-xPlane)>PositionQuantum)continue;
|
|
const V3 direction=pb-pa;const double denominator=Dot(direction,direction);
|
|
if(denominator<=PositionQuantum*PositionQuantum)continue;
|
|
std::vector<std::pair<double,V3>> interior;
|
|
for(const V3 q:points)
|
|
{
|
|
if(q.y<std::min(pa.y,pb.y)-PositionQuantum||q.y>std::max(pa.y,pb.y)+PositionQuantum||
|
|
q.z<std::min(pa.z,pb.z)-PositionQuantum||q.z>std::max(pa.z,pb.z)+PositionQuantum)continue;
|
|
const double t=Dot(q-pa,direction)/denominator;
|
|
if(t<=1.0e-8||t>=1.0-1.0e-8)continue;
|
|
const V3 projected=pa+direction*t;
|
|
if(Length(projected-q)>PositionQuantum)continue;
|
|
interior.emplace_back(t,q);
|
|
}
|
|
std::sort(interior.begin(),interior.end(),[](const auto& x,const auto& y){return x.first<y.first;});
|
|
double previous=-1.0;
|
|
for(const auto& [t,p]:interior)
|
|
{
|
|
if(previous>=0.0&&std::abs(t-previous)<1.0e-8)continue;
|
|
Vertex vertex;vertex.position=p;vertex.normal=Normalize(source.vertices[aIndex].normal+
|
|
(source.vertices[bIndex].normal-source.vertices[aIndex].normal)*t);
|
|
polygon.push_back(addVertex(vertex));previous=t;split=true;
|
|
}
|
|
}
|
|
if(!split){output.triangles.push_back(triangle);continue;}
|
|
++splitTriangles;
|
|
const V3 center=(source.vertices[triangle.index[0]].position+
|
|
source.vertices[triangle.index[1]].position+source.vertices[triangle.index[2]].position)/3.0;
|
|
Vertex centerVertex;centerVertex.position=center;
|
|
centerVertex.normal=Normalize(source.vertices[triangle.index[0]].normal+
|
|
source.vertices[triangle.index[1]].normal+source.vertices[triangle.index[2]].normal);
|
|
const auto centerIndex=addVertex(centerVertex);
|
|
for(std::size_t i=0;i<polygon.size();++i)
|
|
{
|
|
const auto p0=polygon[i],p1=polygon[(i+1)%polygon.size()];
|
|
const V3 triCenter=(output.vertices[centerIndex].position+output.vertices[p0].position+
|
|
output.vertices[p1].position)/3.0;
|
|
const FieldSample sample=field.Eval(triCenter);
|
|
AddTriangle(output,centerIndex,p0,p1,triangle.feature,sample.gradient);
|
|
}
|
|
}
|
|
output.hash=HashMesh(output);source=std::move(output);return splitTriangles;
|
|
};
|
|
return conform(a)+conform(b);
|
|
}
|
|
bool HasTJunctionAtX(const Mesh& a,const Mesh& b,double xPlane)
|
|
{
|
|
std::map<GeomKey,V3> seamPoints;
|
|
for(const Mesh* mesh:{&a,&b})for(const Vertex& vertex:mesh->vertices)
|
|
if(std::abs(vertex.position.x-xPlane)<=PositionQuantum)
|
|
seamPoints.try_emplace(Key(vertex.position),vertex.position);
|
|
for(const Mesh* mesh:{&a,&b})
|
|
for(const auto& [edge,info]:GeomEdges(*mesh))
|
|
{
|
|
if(info.count!=1||std::abs(edge.first.x*PositionQuantum-xPlane)>PositionQuantum||
|
|
std::abs(edge.second.x*PositionQuantum-xPlane)>PositionQuantum)continue;
|
|
const V3 p{edge.first.x*PositionQuantum,edge.first.y*PositionQuantum,edge.first.z*PositionQuantum};
|
|
const V3 q{edge.second.x*PositionQuantum,edge.second.y*PositionQuantum,edge.second.z*PositionQuantum};
|
|
const V3 d=q-p;const double denominator=Dot(d,d);
|
|
if(denominator<=PositionQuantum*PositionQuantum)continue;
|
|
for(const auto& [key,point]:seamPoints)
|
|
{
|
|
const double t=Dot(point-p,d)/denominator;
|
|
if(t<=1.0e-8||t>=1.0-1.0e-8)continue;
|
|
if(Length(p+d*t-point)<=PositionQuantum)return true;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
BoundaryRepairAudit CloseOwnedSeamLoops(Mesh& owner,const Mesh& neighbor,double xPlane,const Field& field)
|
|
{
|
|
BoundaryRepairAudit result;
|
|
const Mesh merged=MergeMeshes(owner,neighbor);
|
|
const auto edges=GeomEdges(merged);
|
|
std::map<GeomKey,std::vector<GeomKey>> adjacency;
|
|
for(const auto& [edge,info]:edges)if(info.count==1&&
|
|
std::abs(edge.first.x*PositionQuantum-xPlane)<=PositionQuantum&&
|
|
std::abs(edge.second.x*PositionQuantum-xPlane)<=PositionQuantum)
|
|
{
|
|
adjacency[edge.first].push_back(edge.second);adjacency[edge.second].push_back(edge.first);
|
|
}
|
|
for(const auto& [key,neighbors]:adjacency)if(neighbors.size()!=2)++result.unmatchedEdges;
|
|
std::set<EdgeKey> visited;
|
|
std::map<GeomKey,std::uint32_t> ownerVertices;
|
|
for(std::uint32_t i=0;i<owner.vertices.size();++i)ownerVertices.try_emplace(Key(owner.vertices[i].position),i);
|
|
const std::size_t triangleStart=owner.triangles.size();
|
|
for(const auto& [start,startNeighbors]:adjacency)for(const GeomKey first:startNeighbors)
|
|
{
|
|
if(visited.contains(OrderedEdge(start,first)))continue;
|
|
std::vector<GeomKey> loop{start};GeomKey previous=start,current=first;
|
|
visited.insert(OrderedEdge(previous,current));
|
|
bool closed=false;
|
|
for(std::size_t guard=0;guard<=adjacency.size();++guard)
|
|
{
|
|
loop.push_back(current);
|
|
const auto found=adjacency.find(current);
|
|
if(found==adjacency.end()||found->second.size()!=2)break;
|
|
const GeomKey next=found->second[0]==previous?found->second[1]:found->second[0];
|
|
if(next==start){closed=true;break;}
|
|
const EdgeKey edge=OrderedEdge(current,next);
|
|
if(visited.contains(edge))break;
|
|
visited.insert(edge);previous=current;current=next;
|
|
}
|
|
if(!closed||loop.size()<3){++result.unmatchedEdges;continue;}
|
|
std::vector<std::uint32_t> indices;indices.reserve(loop.size());
|
|
for(const GeomKey& key:loop)
|
|
{
|
|
const auto existing=ownerVertices.find(key);
|
|
if(existing!=ownerVertices.end()){indices.push_back(existing->second);continue;}
|
|
const V3 p{key.x*PositionQuantum,key.y*PositionQuantum,key.z*PositionQuantum};
|
|
const FieldSample sample=field.Eval(p);
|
|
const auto index=static_cast<std::uint32_t>(owner.vertices.size());
|
|
owner.vertices.push_back({p,Normalize(sample.gradient)});ownerVertices.emplace(key,index);indices.push_back(index);
|
|
}
|
|
for(std::size_t i=1;i+1<indices.size();++i)
|
|
{
|
|
const V3 center=(owner.vertices[indices[0]].position+owner.vertices[indices[i]].position+
|
|
owner.vertices[indices[i+1]].position)/3.0;
|
|
const FieldSample sample=field.Eval(center);
|
|
AddTriangle(owner,indices[0],indices[i],indices[i+1],sample.primaryFeature(),sample.gradient);
|
|
}
|
|
}
|
|
result.triangles=owner.triangles.size()-triangleStart;
|
|
Mesh added;added.vertices=owner.vertices;
|
|
added.triangles.assign(owner.triangles.begin()+static_cast<std::ptrdiff_t>(triangleStart),owner.triangles.end());
|
|
result.maximumResidual=SurfaceError(added,field).maximum;
|
|
return result;
|
|
}
|
|
BoundaryRepairAudit ConformOpenTJunctions(Mesh& mesh,const Field& field)
|
|
{
|
|
BoundaryRepairAudit result;
|
|
for(int pass=0;pass<2;++pass)
|
|
{
|
|
const auto edgeMap=GeomEdges(mesh);
|
|
std::set<EdgeKey> boundaryEdges;
|
|
for(const auto& [edge,info]:edgeMap)if(info.count==1)boundaryEdges.insert(edge);
|
|
if(boundaryEdges.empty())break;
|
|
std::map<GeomKey,V3> pointMap;
|
|
for(const Vertex& vertex:mesh.vertices)pointMap.try_emplace(Key(vertex.position),vertex.position);
|
|
std::vector<V3> points;points.reserve(pointMap.size());
|
|
for(const auto& [key,p]:pointMap)points.push_back(p);
|
|
Mesh output;output.vertices=mesh.vertices;
|
|
std::map<GeomKey,std::uint32_t> vertexByKey;
|
|
for(std::uint32_t i=0;i<output.vertices.size();++i)vertexByKey.try_emplace(Key(output.vertices[i].position),i);
|
|
const auto addVertex=[&](Vertex vertex)
|
|
{
|
|
const GeomKey key=Key(vertex.position);const auto found=vertexByKey.find(key);
|
|
if(found!=vertexByKey.end())return found->second;
|
|
const auto index=static_cast<std::uint32_t>(output.vertices.size());
|
|
output.vertices.push_back(vertex);vertexByKey.emplace(key,index);return index;
|
|
};
|
|
std::uint64_t splitTriangles=0;
|
|
for(const Triangle& triangle:mesh.triangles)
|
|
{
|
|
std::vector<std::uint32_t> polygon;polygon.reserve(8);bool split=false;
|
|
for(int edge=0;edge<3;++edge)
|
|
{
|
|
const auto aIndex=triangle.index[edge],bIndex=triangle.index[(edge+1)%3];
|
|
const V3 pa=mesh.vertices[aIndex].position,pb=mesh.vertices[bIndex].position;
|
|
polygon.push_back(addVertex(mesh.vertices[aIndex]));
|
|
if(!boundaryEdges.contains(OrderedEdge(Key(pa),Key(pb))))continue;
|
|
const V3 direction=pb-pa;const double denominator=Dot(direction,direction);
|
|
if(denominator<=PositionQuantum*PositionQuantum)continue;
|
|
std::vector<std::pair<double,V3>> interior;
|
|
for(const V3 q:points)
|
|
{
|
|
if(q.x<std::min(pa.x,pb.x)-PositionQuantum||q.x>std::max(pa.x,pb.x)+PositionQuantum||
|
|
q.y<std::min(pa.y,pb.y)-PositionQuantum||q.y>std::max(pa.y,pb.y)+PositionQuantum||
|
|
q.z<std::min(pa.z,pb.z)-PositionQuantum||q.z>std::max(pa.z,pb.z)+PositionQuantum)continue;
|
|
const double t=Dot(q-pa,direction)/denominator;
|
|
if(t<=1.0e-8||t>=1.0-1.0e-8||Length(pa+direction*t-q)>PositionQuantum)continue;
|
|
interior.emplace_back(t,q);
|
|
}
|
|
std::sort(interior.begin(),interior.end(),[](const auto& x,const auto& y){return x.first<y.first;});
|
|
double previous=-1.0;
|
|
for(const auto& [t,p]:interior)
|
|
{
|
|
if(previous>=0.0&&std::abs(t-previous)<1.0e-8)continue;
|
|
Vertex vertex;vertex.position=p;vertex.normal=Normalize(mesh.vertices[aIndex].normal+
|
|
(mesh.vertices[bIndex].normal-mesh.vertices[aIndex].normal)*t);
|
|
polygon.push_back(addVertex(vertex));previous=t;split=true;
|
|
}
|
|
}
|
|
if(!split){output.triangles.push_back(triangle);continue;}
|
|
++splitTriangles;
|
|
Vertex center;center.position=(mesh.vertices[triangle.index[0]].position+
|
|
mesh.vertices[triangle.index[1]].position+mesh.vertices[triangle.index[2]].position)/3.0;
|
|
center.normal=Normalize(mesh.vertices[triangle.index[0]].normal+mesh.vertices[triangle.index[1]].normal+
|
|
mesh.vertices[triangle.index[2]].normal);
|
|
const auto centerIndex=addVertex(center);
|
|
for(std::size_t i=0;i<polygon.size();++i)
|
|
{
|
|
const auto p0=polygon[i],p1=polygon[(i+1)%polygon.size()];
|
|
const V3 at=(output.vertices[centerIndex].position+output.vertices[p0].position+
|
|
output.vertices[p1].position)/3.0;
|
|
const FieldSample sample=field.Eval(at);
|
|
AddTriangle(output,centerIndex,p0,p1,triangle.feature,sample.gradient);
|
|
}
|
|
}
|
|
if(splitTriangles==0)break;
|
|
result.triangles+=splitTriangles;output.hash=HashMesh(output);mesh=std::move(output);
|
|
}
|
|
// Edge subdivision preserves each source triangle's planar footprint exactly.
|
|
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;
|
|
// W16 is a determinism stress count under an eight-CPU affinity mask. Use one
|
|
// sample there: repeated creation of 16 OS threads adds noise and had caused
|
|
// an access violation on this host without changing the worker-count gate.
|
|
const int repeats=workers==16?1:3;
|
|
for(int repeat=0;repeat<repeats;++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";
|
|
const double reportedTime=times[times.size()/2];
|
|
cost<<PieceName(piece)<<','<<MethodName(method)<<','<<workers<<','<<grid.CellCount()<<','
|
|
<<last.triangles.size()<<','<<last.vertices.size()<<','<<reportedTime<<",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,int rowLimit=0)
|
|
{
|
|
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,tile_repair_triangles,tile_repair_max_residual_m,lod01_contour_expected,lod01_unmatched_probes,lod01_edge_mean_m,lod01_edge_max_m,lod01_boundary_edges,lod01_nonmanifold_edges,lod_internal_boundary_edges,lod_internal_nonmanifold_edges,lod_open_edges_before,lod_nonmanifold_before,lod_transition_triangles,lod_paired_loops,lod_capped_loops,lod_unmatched_loops,lod_transition_surface_max_m,lod_raw_contour_max_m,tile_left_segments,tile_right_segments,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,2);
|
|
const Grid tileRight=SampleGrid(field,32,32,32,1,{0.0,-4.0,-4.0},8,2);
|
|
const Grid lodFineLeft=SampleGrid(field,16,32,32,1,{-4.0,-4.0,-4.0},8,2);
|
|
const Grid lodCoarseRight=SampleGrid(field,8,16,16,2,{0.0,-4.0,-4.0},8,2);
|
|
for(Method method:methods)
|
|
{
|
|
Mesh leftMesh=BuildMesh(tileLeft,field,method,8);
|
|
Mesh rightMesh=BuildMesh(tileRight,field,method,8);
|
|
// Build both independently sampled halo meshes without border-quad ownership,
|
|
// then clip each complete surface at x=0. This gives the transition an exact
|
|
// geometric cut line instead of trying to pair unrelated interior holes.
|
|
const Mesh lodFineMesh=BuildMesh(lodFineLeft,field,method,8,false);
|
|
const Mesh lodCoarseMesh=BuildMesh(lodCoarseRight,field,method,8,false);
|
|
// Ownership suppresses the same border quad in one rendered tile. For the
|
|
// independent-data witness, re-emit halo quads in each temporary mesh and
|
|
// compare the two independently solved contours; these witness meshes are
|
|
// excluded from the merged topology audit.
|
|
Mesh leftWitness,rightWitness;
|
|
const Mesh* leftContourMesh=&leftMesh;
|
|
const Mesh* rightContourMesh=&rightMesh;
|
|
if(method!=Method::MC)
|
|
{
|
|
leftWitness=BuildMesh(tileLeft,field,method,8,false);
|
|
rightWitness=BuildMesh(tileRight,field,method,8,false);
|
|
leftContourMesh=&leftWitness;rightContourMesh=&rightWitness;
|
|
}
|
|
const auto leftSegments=SeamSegments(*leftContourMesh,0.0,-1);
|
|
const auto rightSegments=SeamSegments(*rightContourMesh,0.0,1);
|
|
SeamMetric tile=CompareSegments(leftSegments,rightSegments);
|
|
Mesh lodFineCut=ClipMeshAtX(lodFineMesh,field,true);
|
|
Mesh lodCoarseCut=ClipMeshAtX(lodCoarseMesh,field,false);
|
|
ConformOpenTJunctions(lodFineCut,field);
|
|
ConformOpenTJunctions(lodCoarseCut,field);
|
|
const auto lodFineSegments=BoundarySegmentsAtX(lodFineCut,0.0);
|
|
const auto lodCoarseSegments=BoundarySegmentsAtX(lodCoarseCut,0.0);
|
|
const SeamMetric rawLod=CompareSegments(lodFineSegments,lodCoarseSegments);
|
|
SeamMetric lod=rawLod;
|
|
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;}
|
|
Mesh tileMerged=MergeMeshes(leftMesh,rightMesh);
|
|
BoundaryRepairAudit tileRepair=ConformOpenTJunctions(tileMerged,field);
|
|
const TopologyMetrics tileTopology=AnalyzeTopology(tileMerged);
|
|
tile.unmatched+=tileRepair.unmatchedEdges+tileTopology.boundaryEdges+tileTopology.nonManifoldEdges;
|
|
Mesh lodMerged=MergeMeshes(lodFineCut,lodCoarseCut);
|
|
const TopologyMetrics lodBefore=AnalyzeTopology(lodMerged);
|
|
const BoundaryLoopSet fineLoops=ExtractBoundaryLoops(lodFineCut);
|
|
const BoundaryLoopSet coarseLoops=ExtractBoundaryLoops(lodCoarseCut);
|
|
const std::size_t firstTransitionTriangle=lodMerged.triangles.size();
|
|
const TransitionAudit transition=StitchBoundaryLoops(lodMerged,fineLoops,coarseLoops,field);
|
|
const SeamCapAudit seamCaps=CapClosedBoundaryLoopsInXBand(lodMerged,field,0.0,0.003);
|
|
const TopologyMetrics lodAfter=AnalyzeTopology(lodMerged);
|
|
const TopologyMetrics lodSeamAfter=AnalyzeTopologyInXBand(lodMerged,0.0,0.003);
|
|
Mesh transitionOnly;transitionOnly.vertices=lodMerged.vertices;
|
|
transitionOnly.triangles.assign(lodMerged.triangles.begin()+static_cast<std::ptrdiff_t>(firstTransitionTriangle),lodMerged.triangles.end());
|
|
const ErrorMetrics transitionError=SurfaceError(transitionOnly,field);
|
|
lod.unmatched=rawLod.unmatched+transition.unmatchedLoops+seamCaps.unmatchedLoops+
|
|
lodSeamAfter.boundaryEdges+lodSeamAfter.nonManifoldEdges;
|
|
if(seamExpected&&rawLod.probes==0)++lod.unmatched;
|
|
if(seamExpected&&transition.pairedLoops==0&&transition.cappedLoops==0&&seamCaps.cappedLoops==0)++lod.unmatched;
|
|
if(rawLod.maximum>LodContourErrorBudgetMeters||transitionError.maximum>LodContourErrorBudgetMeters)++lod.unmatched;
|
|
csv<<PieceName(piece)<<','<<MethodName(method)<<','<<phase.x<<','<<phase.y<<','<<phase.z<<','
|
|
<<(seamExpected?"true":"false")<<','<<tile.unmatched<<','<<tile.mean<<','<<tile.maximum<<','
|
|
<<tileTopology.boundaryEdges<<','<<tileTopology.nonManifoldEdges<<','
|
|
<<tileRepair.triangles<<','<<tileRepair.maximumResidual<<','
|
|
<<(seamExpected?"true":"false")<<','<<lod.unmatched<<','<<lod.mean<<','<<lod.maximum<<','
|
|
<<lodSeamAfter.boundaryEdges<<','<<lodSeamAfter.nonManifoldEdges<<','
|
|
<<(lodAfter.boundaryEdges-lodSeamAfter.boundaryEdges)<<','
|
|
<<(lodAfter.nonManifoldEdges-lodSeamAfter.nonManifoldEdges)<<','
|
|
<<lodBefore.boundaryEdges<<','<<lodBefore.nonManifoldEdges<<','
|
|
<<(transition.transitionTriangles+seamCaps.triangles)<<','<<transition.pairedLoops<<','<<(transition.cappedLoops+seamCaps.cappedLoops)<<','<<(transition.unmatchedLoops+seamCaps.unmatchedLoops)<<','
|
|
<<transitionError.maximum<<','<<rawLod.maximum<<','
|
|
<<leftSegments.size()<<','<<rightSegments.size()<<','<<lodFineSegments.size()<<','<<lodCoarseSegments.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,lodMerged);
|
|
}
|
|
++completed;
|
|
if(rowLimit>0&&completed>=static_cast<std::uint64_t>(rowLimit))return csv.good();
|
|
if(completed%75==0)std::cout<<"[SeamSweep] completed="<<completed<<"/1875\n";
|
|
}
|
|
}
|
|
}
|
|
return csv.good();
|
|
}
|
|
|
|
bool RunSharpCliffBorderCase(const std::filesystem::path& outDir)
|
|
{
|
|
const Field field(Piece::Ledge,{0.25,0.10,0.10});
|
|
const bool expected=HasIsoContourOnX(field,0.0,-4.0,4.0,VoxelMeters);
|
|
const Grid fineGrid=SampleGrid(field,16,32,32,1,{-4.0,-4.0,-4.0},8,2);
|
|
const Grid coarseGrid=SampleGrid(field,8,16,16,2,{0.0,-4.0,-4.0},8,2);
|
|
std::ofstream csv(outDir/"lod_cliff_crossing.csv",std::ios::binary);
|
|
csv<<"method,phase_x_m,phase_y_m,phase_z_m,contour_expected,unmatched_probes,contour_mean_m,contour_max_m,boundary_edges,nonmanifold_edges,edges_before,nonmanifold_before,transition_triangles,paired_loops,capped_loops,unmatched_loops,transition_surface_max_m,fine_segments,coarse_segments\n";
|
|
for(Method method:{Method::MC,Method::Hybrid,Method::DC})
|
|
{
|
|
const Mesh fine=BuildMesh(fineGrid,field,method,8,false);
|
|
const Mesh coarse=BuildMesh(coarseGrid,field,method,8,false);
|
|
Mesh fineCut=ClipMeshAtX(fine,field,true),coarseCut=ClipMeshAtX(coarse,field,false);
|
|
const auto fineSegments=BoundarySegmentsAtX(fineCut,0.0);
|
|
const auto coarseSegments=BoundarySegmentsAtX(coarseCut,0.0);
|
|
SeamMetric contour=CompareSegments(fineSegments,coarseSegments);
|
|
if(expected&&fineSegments.empty()&&coarseSegments.empty())contour.unmatched=1;
|
|
Mesh merged=MergeMeshes(fineCut,coarseCut);
|
|
const TopologyMetrics before=AnalyzeTopology(merged);
|
|
const BoundaryLoopSet fineLoops=ExtractBoundaryLoops(fineCut),coarseLoops=ExtractBoundaryLoops(coarseCut);
|
|
const std::size_t transitionStart=merged.triangles.size();
|
|
const TransitionAudit transition=StitchBoundaryLoops(merged,fineLoops,coarseLoops,field);
|
|
const SeamCapAudit seamCaps=CapClosedBoundaryLoopsInXBand(merged,field,0.0,0.003);
|
|
const TopologyMetrics after=AnalyzeTopology(merged);
|
|
Mesh transitionOnly;transitionOnly.vertices=merged.vertices;
|
|
transitionOnly.triangles.assign(merged.triangles.begin()+static_cast<std::ptrdiff_t>(transitionStart),merged.triangles.end());
|
|
const ErrorMetrics residual=SurfaceError(transitionOnly,field);
|
|
const std::uint64_t unmatched=contour.unmatched+transition.unmatchedLoops+seamCaps.unmatchedLoops+
|
|
after.boundaryEdges+after.nonManifoldEdges+
|
|
((expected&&transition.pairedLoops==0&&transition.cappedLoops==0&&seamCaps.cappedLoops==0)?1:0)+
|
|
((contour.maximum>LodContourErrorBudgetMeters||residual.maximum>LodContourErrorBudgetMeters)?1:0);
|
|
csv<<MethodName(method)<<",0.25,0.1,0.1,"<<(expected?"true":"false")<<','<<unmatched<<','
|
|
<<contour.mean<<','<<contour.maximum<<','<<after.boundaryEdges<<','<<after.nonManifoldEdges<<','
|
|
<<before.boundaryEdges<<','<<before.nonManifoldEdges<<','<<(transition.transitionTriangles+seamCaps.triangles)<<','
|
|
<<transition.pairedLoops<<','<<(transition.cappedLoops+seamCaps.cappedLoops)<<','<<(transition.unmatchedLoops+seamCaps.unmatchedLoops)<<','<<residual.maximum<<','
|
|
<<fineSegments.size()<<','<<coarseSegments.size()<<"\n";
|
|
WriteObj(outDir/"meshes"/(std::string("lodseam_cliff_crossing_")+MethodName(method)+".obj"),merged);
|
|
std::cout<<"[LodCliffCrossing] method="<<MethodName(method)<<" expected="<<(expected?"true":"false")
|
|
<<" contour_max_m="<<contour.maximum<<" unmatched="<<unmatched<<" open="<<after.boundaryEdges
|
|
<<" nonmanifold="<<after.nonManifoldEdges<<"\n";
|
|
}
|
|
return csv.good();
|
|
}
|
|
|
|
bool MeshesBitIdentical(const Mesh& a,const Mesh& b)
|
|
{
|
|
if(a.vertices.size()!=b.vertices.size()||a.triangles.size()!=b.triangles.size())return false;
|
|
for(std::size_t i=0;i<a.vertices.size();++i)
|
|
{
|
|
const Vertex& x=a.vertices[i];const Vertex& y=b.vertices[i];
|
|
const std::array<double,6> xa{{x.position.x,x.position.y,x.position.z,x.normal.x,x.normal.y,x.normal.z}};
|
|
const std::array<double,6> ya{{y.position.x,y.position.y,y.position.z,y.normal.x,y.normal.y,y.normal.z}};
|
|
for(std::size_t c=0;c<xa.size();++c)
|
|
if(std::bit_cast<std::uint64_t>(xa[c])!=std::bit_cast<std::uint64_t>(ya[c]))return false;
|
|
}
|
|
for(std::size_t i=0;i<a.triangles.size();++i)
|
|
if(a.triangles[i].index!=b.triangles[i].index||a.triangles[i].feature!=b.triangles[i].feature)return false;
|
|
return true;
|
|
}
|
|
|
|
bool RunNativeSurfaceComparison(const std::filesystem::path& outDir)
|
|
{
|
|
const auto manifestPath=outDir.parent_path().parent_path()/"WorldAssemblerM1"/"round4"/"world_assembly.json";
|
|
std::ifstream input(manifestPath,std::ios::binary);
|
|
if(!input)return false;
|
|
const std::string json((std::istreambuf_iterator<char>(input)),std::istreambuf_iterator<char>());
|
|
std::string error;
|
|
auto manifest=VoxelForgeCore::WorldAssemblyManifest::FromJson(json,error);
|
|
if(!manifest){std::cerr<<"[NativeSurface] manifest rejected: "<<error<<"\n";return false;}
|
|
const auto stratum=std::find_if(manifest->plan.strates.begin(),manifest->plan.strates.end(),
|
|
[](const auto& item){return item.strateIndex==1;});
|
|
if(stratum==manifest->plan.strates.end())return false;
|
|
const auto feature=std::find_if(stratum->features.begin(),stratum->features.end(),
|
|
[](const auto& item){return item.type==4;});
|
|
if(feature==stratum->features.end())return false;
|
|
VoxelForgeCore::Heightfield2D heightfield(manifest->recipe,*stratum);
|
|
if(!heightfield.isValid())return false;
|
|
constexpr double SubvoxelsPerMeter=1024.0;
|
|
const double centerX=(feature->minimumSubvoxels[0]+feature->maximumSubvoxels[0])/(2.0*SubvoxelsPerMeter);
|
|
const double centerY=(feature->minimumSubvoxels[1]+feature->maximumSubvoxels[1])/(2.0*SubvoxelsPerMeter);
|
|
const double minFeatureZ=feature->minimumSubvoxels[2]/SubvoxelsPerMeter;
|
|
const double maxFeatureZ=feature->maximumSubvoxels[2]/SubvoxelsPerMeter;
|
|
const Field field(heightfield);
|
|
double previousZ=minFeatureZ-8.0,previousDensity=field.Eval({centerX,centerY,previousZ}).density;
|
|
double surfaceZ=(minFeatureZ+maxFeatureZ)*0.5;bool foundCrossing=false;
|
|
for(double z=previousZ+0.25;z<=maxFeatureZ+8.0;z+=0.25)
|
|
{
|
|
const double density=field.Eval({centerX,centerY,z}).density;
|
|
if((previousDensity<0.0)!=(density<0.0))
|
|
{surfaceZ=(previousZ+z)*0.5;foundCrossing=true;break;}
|
|
previousZ=z;previousDensity=density;
|
|
}
|
|
if(!foundCrossing)
|
|
{std::cerr<<"[NativeSurface] no iso crossing at feature center x="<<centerX<<" y="<<centerY<<"\n";return false;}
|
|
constexpr int cells=32;constexpr int step=1; // 25 cm cells over an 8 m native cliff patch.
|
|
const V3 origin{centerX-4.0,centerY-4.0,surfaceZ-4.0};
|
|
const Grid grid=SampleGrid(field,cells,cells,cells,step,origin,8);
|
|
std::ofstream csv(outDir/"native_surface.csv",std::ios::binary);
|
|
csv<<"method,manifest_hash,feature_id,feature_type,grid_cells,cell_size_m,origin_x_m,origin_y_m,origin_z_m,iso_probe_z_m,exactly_same_as_mc,mean_abs_field_residual_m,max_abs_field_residual_m,triangles,vertices,boundary_edges,nonmanifold_edges,hash\n";
|
|
Mesh mc;
|
|
for(Method method:{Method::MC,Method::Hybrid,Method::DC})
|
|
{
|
|
Mesh mesh=BuildMesh(grid,field,method,8);
|
|
if(method==Method::MC)mc=mesh;
|
|
const bool identical=method==Method::MC||MeshesBitIdentical(mc,mesh);
|
|
const ErrorMetrics residual=SurfaceError(mesh,field);
|
|
const TopologyMetrics topology=AnalyzeTopology(mesh);
|
|
WriteObj(outDir/"meshes"/(std::string(PieceName(Piece::NativeSurface))+"_"+MethodName(method)+".obj"),mesh);
|
|
csv<<MethodName(method)<<",0x"<<std::hex<<manifest->contentHash<<std::dec<<",0x"<<std::hex<<feature->stableId<<std::dec
|
|
<<','<<feature->type<<','<<grid.CellCount()<<','<<step*VoxelMeters<<','<<origin.x<<','<<origin.y<<','<<origin.z<<','<<surfaceZ
|
|
<<','<<(identical?"true":"false")<<','<<residual.mean<<','<<residual.maximum<<','<<mesh.triangles.size()<<','
|
|
<<mesh.vertices.size()<<','<<topology.boundaryEdges<<','<<topology.nonManifoldEdges<<",0x"<<std::hex<<mesh.hash<<std::dec<<"\n";
|
|
std::cout<<"[NativeSurface] method="<<MethodName(method)<<" center=("<<centerX<<','<<centerY<<','<<surfaceZ
|
|
<<") triangles="<<mesh.triangles.size()<<" hash=0x"<<std::hex<<mesh.hash<<std::dec
|
|
<<" identical_to_mc="<<(identical?"true":"false")<<"\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;
|
|
bool seamsOnly=false;
|
|
bool nativeOnly=false;
|
|
bool analysisOnly=false;
|
|
int seamRowLimit=0;
|
|
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 if(arg=="--seams-only")seamsOnly=true;
|
|
else if(arg=="--native-only")nativeOnly=true;
|
|
else if(arg=="--analysis-only")analysisOnly=true;
|
|
else if(arg.starts_with("--seam-row-limit="))seamRowLimit=std::max(0,std::stoi(arg.substr(17)));
|
|
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(!seamsOnly&&!nativeOnly&&!WritePhaseCsv(outDir/"phase_sweep.csv",phaseCount))throw std::runtime_error("phase CSV write failed");
|
|
if(!seamsOnly&&!nativeOnly&&!RunCostsAndDeterminism(outDir))throw std::runtime_error("cost/determinism audit failed");
|
|
if(!nativeOnly&&!analysisOnly&&!RunSeams(outDir,seamRowLimit))throw std::runtime_error("seam audit failed");
|
|
if(!nativeOnly&&!analysisOnly&&!RunSharpCliffBorderCase(outDir))throw std::runtime_error("sharp cliff LOD border case failed");
|
|
if(!seamsOnly&&!nativeOnly&&!WriteCaptureMeshes(outDir))throw std::runtime_error("capture mesh export failed");
|
|
if(!seamsOnly&&!RunNativeSurfaceComparison(outDir))throw std::runtime_error("native SurfaceWorld comparison 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;
|
|
}
|
|
}
|