cleaning of isotropic remeshing and small fix to adaptivity

This commit is contained in:
korialis 2021-11-26 16:25:17 +01:00
parent 0c4e210bba
commit f2ba3e973e
1 changed files with 185 additions and 251 deletions

View File

@ -155,7 +155,7 @@ private:
int iter = 0; int iter = 0;
do do
{ {
vcg::tri::UpdateTopology<MeshType>::FaceFace(m); // vcg::tri::UpdateTopology<MeshType>::FaceFace(m);
vcg::tri::UnMarkAll(m); vcg::tri::UnMarkAll(m);
count = 0; count = 0;
@ -163,44 +163,49 @@ private:
{ {
FaceType & f = m.face[i]; FaceType & f = m.face[i];
ScalarType quality = vcg::QualityRadii(f.cP(0), f.cP(1), f.cP(2)); const ScalarType quality = vcg::QualityRadii(f.cP(0), f.cP(1), f.cP(2));
if (quality <= 0.001) if (quality <= 0.001)
{ {
//find longest edge //find longest edge
double edges[3]; const double edges[3] = {
edges[0] = vcg::Distance(f.cP(0), f.cP(1)); vcg::Distance(f.cP(0), f.cP(1)),
edges[1] = vcg::Distance(f.cP(1), f.cP(2)); vcg::Distance(f.cP(1), f.cP(2)),
edges[2] = vcg::Distance(f.cP(2), f.cP(0)); vcg::Distance(f.cP(2), f.cP(0))
};
ScalarType smallestEdge = std::min(edges[0], std::min(edges[1], edges[2])); const double smallestEdge = std::min(edges[0], std::min(edges[1], edges[2]));
int longestIdx = std::find(edges, edges+3, std::max(std::max(edges[0], edges[1]), edges[2])) - (edges); const int longestIdx = int(std::find(edges, edges+3, std::max(std::max(edges[0], edges[1]), edges[2])) - (edges));
if (vcg::tri::IsMarked(m, f.V2(longestIdx))) if (vcg::tri::IsMarked(m, f.V2(longestIdx)))
continue; continue;
auto f1 = f.cFFp(longestIdx); auto f1 = f.cFFp(longestIdx);
vcg::tri::Mark(m,f.V2(longestIdx)); vcg::tri::Mark(m, f.V2(longestIdx));
if (!vcg::face::IsBorder(f, longestIdx) && vcg::face::IsManifold(f, longestIdx) && vcg::face::checkFlipEdgeNotManifold<FaceType>(f, longestIdx)) { if (!vcg::face::IsBorder(f, longestIdx) && vcg::face::IsManifold(f, longestIdx) && vcg::face::checkFlipEdgeNotManifold<FaceType>(f, longestIdx)) {
// Check if EdgeFlipping improves quality // Check if EdgeFlipping improves quality
FacePointer g = f.FFp(longestIdx); int k = f.FFi(longestIdx); const FacePointer g = f.FFp(longestIdx);
vcg::Triangle3<ScalarType> t1(f.P(longestIdx), f.P1(longestIdx), f.P2(longestIdx)), t2(g->P(k), g->P1(k), g->P2(k)), const int k = f.FFi(longestIdx);
t3(f.P(longestIdx), g->P2(k), f.P2(longestIdx)), t4(g->P(k), f.P2(longestIdx), g->P2(k));
auto n1 = vcg::TriangleNormal(t1); const vcg::Triangle3<ScalarType> t1(f.P(longestIdx), f.P1(longestIdx), f.P2(longestIdx));
auto n2 = vcg::TriangleNormal(t2); const vcg::Triangle3<ScalarType> t2(g->P(k), g->P1(k), g->P2(k));
auto n3 = vcg::TriangleNormal(t3); const vcg::Triangle3<ScalarType> t3(f.P(longestIdx), g->P2(k), f.P2(longestIdx));
auto n4 = vcg::TriangleNormal(t4); const vcg::Triangle3<ScalarType> t4(g->P(k), f.P2(longestIdx), g->P2(k));
auto biggestSmallest = vcg::DoubleArea(t1) > vcg::DoubleArea(t2) ? std::make_pair(t1, t2) : std::make_pair(t2, t1); const auto n1 = vcg::TriangleNormal(t1);
auto areaRatio = vcg::DoubleArea(biggestSmallest.first) / vcg::DoubleArea(biggestSmallest.second); const auto n2 = vcg::TriangleNormal(t2);
const auto n3 = vcg::TriangleNormal(t3);
const auto n4 = vcg::TriangleNormal(t4);
const auto biggestSmallest = vcg::DoubleArea(t1) > vcg::DoubleArea(t2) ? std::make_pair(t1, t2) : std::make_pair(t2, t1);
const auto areaRatio = vcg::DoubleArea(biggestSmallest.first) / vcg::DoubleArea(biggestSmallest.second);
bool normalCheck = true; bool normalCheck = true;
// if (n1.Norm() > 0.001 && n2.Norm() > 0.001) // if (n1.Norm() > 0.001 && n2.Norm() > 0.001)
{ {
auto referenceNormal = vcg::NormalizedTriangleNormal(biggestSmallest.first); const auto referenceNormal = vcg::NormalizedTriangleNormal(biggestSmallest.first);
normalCheck &= vcg::NormalizedTriangleNormal(t3) * referenceNormal >= 0.95; normalCheck &= vcg::NormalizedTriangleNormal(t3) * referenceNormal >= 0.95;
normalCheck &= vcg::NormalizedTriangleNormal(t4) * referenceNormal >= 0.95; normalCheck &= vcg::NormalizedTriangleNormal(t4) * referenceNormal >= 0.95;
@ -240,9 +245,9 @@ private:
vcg::tri::Clean<MeshType>::RemoveUnreferencedVertex(m); vcg::tri::Clean<MeshType>::RemoveUnreferencedVertex(m);
vcg::tri::Allocator<MeshType>::CompactEveryVector(m); vcg::tri::Allocator<MeshType>::CompactEveryVector(m);
vcg::tri::UpdateTopology<MeshType>::FaceFace(m); // vcg::tri::UpdateTopology<MeshType>::FaceFace(m);
removeColinearFaces(m, params); removeColinearFaces(m, params);
vcg::tri::UpdateTopology<MeshType>::FaceFace(m); // vcg::tri::UpdateTopology<MeshType>::FaceFace(m);
} }
public: public:
@ -262,7 +267,6 @@ public:
assert(&toRemesh != &toProject); assert(&toRemesh != &toProject);
params.stat.Reset(); params.stat.Reset();
tri::UpdateBounding<MeshType>::Box(toRemesh); tri::UpdateBounding<MeshType>::Box(toRemesh);
{ {
@ -287,15 +291,15 @@ public:
{ {
if(cb) cb(100*i/params.iter, "Remeshing"); if(cb) cb(100*i/params.iter, "Remeshing");
if (params.adapt) if (params.adapt)
{ {
computeQualityDistFromRadii(toRemesh); computeQualityDistFromRadii(toRemesh);
tri::Smooth<MeshType>::VertexQualityLaplacian(toRemesh, 2); vcg::tri::Smooth<MeshType>::VertexQualityLaplacian(toRemesh, 2);
} }
if(params.splitFlag) if(params.splitFlag)
SplitLongEdges(toRemesh, params); SplitLongEdges(toRemesh, params);
#ifdef DEBUG_CREASE #ifdef DEBUG_CREASE
debug_crease(toRemesh, std::string("after_ref"), i); debug_crease(toRemesh, std::string("after_ref"), i);
#endif #endif
@ -334,25 +338,24 @@ public:
if (p.IsBorder()) if (p.IsBorder())
p.F()->SetFaceEdgeS(p.E()); p.F()->SetFaceEdgeS(p.E());
// if((p.F1Flip() > p.F())) const FaceType *ff = p.F();
{ const FaceType *ffAdj = p.FFlip();
FaceType *ff = p.F();
FaceType *ffAdj = p.FFlip();
double quality = vcg::QualityRadii(ff->cP(0), ff->cP(1), ff->cP(2)); const double quality = vcg::QualityRadii(ff->cP(0), ff->cP(1), ff->cP(2));
double qualityAdj = vcg::QualityRadii(ffAdj->cP(0), ffAdj->cP(1), ffAdj->cP(2)); const double qualityAdj = vcg::QualityRadii(ffAdj->cP(0), ffAdj->cP(1), ffAdj->cP(2));
bool qualityCheck = quality > 0.00000001 && qualityAdj > 0.00000001; const bool qualityCheck = quality > 0.00000001 && qualityAdj > 0.00000001;
// bool areaCheck = vcg::DoubleArea(*ff) > 0.000001 && vcg::DoubleArea(*ffAdj) > 0.000001; // bool areaCheck = vcg::DoubleArea(*ff) > 0.000001 && vcg::DoubleArea(*ffAdj) > 0.000001;
if ((forceTag || !params.userSelectedCreases) && (testCreaseEdge(p, params.creaseAngleCosThr) /*&& areaCheck*//* && qualityCheck*/) || p.IsBorder()) if ((forceTag || !params.userSelectedCreases) && (testCreaseEdge(p, params.creaseAngleCosThr) /*&& areaCheck*/ /* && qualityCheck*/) || p.IsBorder())
{ {
PosType pp = p; PosType pp = p;
std::vector<FacePointer> faces; std::vector<FacePointer> faces;
std::vector<int> edges; std::vector<int> edges;
bool allOk = true; bool allOk = true;
do { do
{
faces.push_back(pp.F()); faces.push_back(pp.F());
edges.push_back(pp.E()); edges.push_back(pp.E());
// pp.F()->SetFaceEdgeS(pp.E()); // pp.F()->SetFaceEdgeS(pp.E());
@ -374,7 +377,7 @@ public:
creaseQueue.push(p); creaseQueue.push(p);
} }
}
}); });
return count; return count;
} }
@ -392,14 +395,14 @@ private:
*/ */
IsotropicRemeshing() {} IsotropicRemeshing() {}
// this returns the value of cos(a) where a is the angle between n0 and n1. (scalar prod is cos(a)) // this returns the value of cos(a) where a is the angle between n0 and n1. (scalar prod is cos(a))
static inline ScalarType fastAngle(Point3<ScalarType> n0, Point3<ScalarType> n1) static inline ScalarType fastAngle(const Point3<ScalarType> & n0, const Point3<ScalarType> & n1)
{ {
return math::Clamp(n0*n1,(ScalarType)-1.0,(ScalarType)1.0); return math::Clamp(n0*n1,(ScalarType)-1.0,(ScalarType)1.0);
} }
// compare the value of the scalar prod with the cos of the crease threshold // compare the value of the scalar prod with the cos of the crease threshold
static inline bool testCreaseEdge(PosType &p, ScalarType creaseCosineThr) static inline bool testCreaseEdge(PosType &p, const ScalarType creaseCosineThr)
{ {
ScalarType angle = fastAngle(NormalizedTriangleNormal(*(p.F())), NormalizedTriangleNormal(*(p.FFlip()))); const ScalarType angle = fastAngle(NormalizedTriangleNormal(*(p.F())), NormalizedTriangleNormal(*(p.FFlip())));
return angle <= creaseCosineThr && angle >= -0.98; return angle <= creaseCosineThr && angle >= -0.98;
// return (angle <= creaseCosineThr && angle >= -creaseCosineThr); // return (angle <= creaseCosineThr && angle >= -creaseCosineThr);
} }
@ -432,16 +435,15 @@ private:
std::vector<FaceType*> ff; std::vector<FaceType*> ff;
face::VFExtendedStarVF(&*vi, 2, ff); face::VFExtendedStarVF(&*vi, 2, ff);
ScalarType tot = 0.f; assert(ff.size() > 0);
auto it = ff.begin();
Point3<ScalarType> fNormal = NormalizedTriangleNormal(**it); const Point3<ScalarType> & fNormal = NormalizedTriangleNormal(**it);
++it;
while(it != ff.end()) const auto tot = std::accumulate(++ff.begin(), ff.end(), 0.d, [&](const Scalartype acc, const FaceType * f) {
{ return acc + (1 - math::Abs(fastAngle(n, NormalizedTriangleNormal(*f))));
tot+= 1-math::Abs(fastAngle(fNormal, NormalizedTriangleNormal(**it))); });
++it;
} vi->Q() = tot / (std::max(1, ((int)ff.size()-1)));
vi->Q() = tot / (ScalarType)(std::max(1, ((int)ff.size()-1)));
vi->SetV(); vi->SetV();
} }
tri::Smooth<MeshType>::VertexQualityLaplacian(m, 3); tri::Smooth<MeshType>::VertexQualityLaplacian(m, 3);
@ -452,8 +454,10 @@ private:
tri::RequirePerVertexQuality(m); tri::RequirePerVertexQuality(m);
tri::UpdateTopology<MeshType>::FaceFace(m); tri::UpdateTopology<MeshType>::FaceFace(m);
// tri::UpdateFlags<MeshType>::VertexClearV(m); // tri::UpdateFlags<MeshType>::VertexClearV(m);
for (size_t i=0;i<m.vert.size();i++)
m.vert[i].IMark()=0; ForEachVertex(m, [&] (VertexType & v) {
v.IMark() = 0;
});
std::vector<VertexPointer> seeds; std::vector<VertexPointer> seeds;
ForEachFace(m, [&] (FaceType & f) { ForEachFace(m, [&] (FaceType & f) {
@ -485,8 +489,8 @@ private:
tri::RequirePerVertexQuality(m); tri::RequirePerVertexQuality(m);
tri::RequirePerFaceQuality(m); tri::RequirePerFaceQuality(m);
ScalarType maxV = 0; ScalarType maxV = std::numeric_limits<ScalarType>::lowest();
ScalarType minV = 10; ScalarType minV = std::numeric_limits<ScalarType>::max();
ForEachFace(m, [&] (FaceType & f) { ForEachFace(m, [&] (FaceType & f) {
f.Q() = 1. - vcg::QualityRadii(f.cP(0), f.cP(1), f.cP(2)); f.Q() = 1. - vcg::QualityRadii(f.cP(0), f.cP(1), f.cP(2));
@ -494,29 +498,22 @@ private:
minV = std::min(minV, f.Q()); minV = std::min(minV, f.Q());
}); });
//normalize vcg::tri::UpdateQuality<MeshType>::VertexFromFace(m);
ForEachFace(m, [&] (FaceType & f) {
f.Q() = std::pow((f.Q() - minV) / (maxV - minV), 2.); maxV = std::numeric_limits<ScalarType>::lowest();
minV = std::numeric_limits<ScalarType>::max();
//normalize in [0,1] with square reshape
ForEachVertex(m, [&] (const VertexType & v) {
maxV = std::max(maxV, v.Q());
minV = std::min(minV, v.Q());
}); });
std::vector<ScalarType> vertMax(m.VN(), 0); const ScalarType vRange = maxV - minV + 0.000001;
std::vector<ScalarType> vertMin(m.VN(), 10);
ForEachFace(m, [&] (FaceType & f) { ForEachVertex(m, [&] (VertexType & v) {
for (int i = 0; i < 3; ++i) v.Q() = std::pow((v.Q() - minV) / vRange, 2.);
{
auto vidx = vcg::tri::Index(m, f.V(i));
vertMax[vidx] = std::max(vertMax[vidx], f.Q());
vertMin[vidx] = std::min(vertMin[vidx], f.Q());
}
}); });
for (size_t v = 0; v < m.VN(); ++v)
{
m.vert[v].Q() = vertMax[v] - vertMin[v];
}
// tri::UpdateQuality<MeshType>::VertexFromFace(m);
} }
static void computeQualityDistFromHeight(MeshType & m, const ScalarType lowerBound, const ScalarType higherBound) static void computeQualityDistFromHeight(MeshType & m, const ScalarType lowerBound, const ScalarType higherBound)
@ -524,10 +521,6 @@ private:
tri::RequirePerVertexQuality(m); tri::RequirePerVertexQuality(m);
tri::RequirePerFaceQuality(m); tri::RequirePerFaceQuality(m);
ScalarType maxV = 0;
ScalarType minV = 10;
ForEachFace(m, [&] (FaceType & f) { ForEachFace(m, [&] (FaceType & f) {
ScalarType minH = std::numeric_limits<ScalarType>::max(); ScalarType minH = std::numeric_limits<ScalarType>::max();
for (int i = 0; i < 3; ++i) for (int i = 0; i < 3; ++i)
@ -549,36 +542,35 @@ private:
4 for border vertices 4 for border vertices
6 for internal vertices 6 for internal vertices
*/ */
static inline int idealValence(PosType &p) static inline int idealValence(const PosType &p)
{ {
if(p.IsBorder()) return 4; if(p.IsBorder()) return 4;
return 6; return 6;
} }
static inline int idealValence(VertexType &v) static inline int idealValence(const VertexType &v)
{ {
if(v.IsB()) return 4; if(v.IsB()) return 4;
return 6; return 6;
} }
static inline int idealValenceSlow(PosType &p) static inline int idealValenceSlow(const PosType &p)
{ {
std::vector<PosType> posVec; std::vector<PosType> posVec;
VFOrderedStarFF(p,posVec); VFOrderedStarFF(p,posVec);
float angleSumRad =0;
for(PosType &ip : posVec) const auto angleSumRad = std::accumulate(posVec.begin, posVec.end(), 0, [](const ScalarType acc, const PosType & p) {
{ return acc + p.AngleRad();
angleSumRad += ip.AngleRad(); });
}
return (int)(std::ceil(angleSumRad / (M_PI/3.0f))); return (int)(std::ceil(angleSumRad / (M_PI/3.0f)));
} }
static bool testHausdorff (MeshType & m, StaticGrid & grid, const std::vector<CoordType> & verts, const ScalarType maxD, const CoordType checkOrientation = CoordType(0,0,0)) static bool testHausdorff (MeshType & m, StaticGrid & grid, const std::vector<CoordType> & verts, const ScalarType maxD, const CoordType & checkOrientation = CoordType(0,0,0))
{ {
for (CoordType v : verts) for (CoordType v : verts)
{ {
CoordType closest, normal, ip; CoordType closest, normal, ip;
ScalarType dist = 0; ScalarType dist = 0;
FaceType* fp = GetClosestFaceBase(m, grid, v, maxD, dist, closest); const FaceType* fp = GetClosestFaceBase(m, grid, v, maxD, dist, closest);
//you can't use this kind of orientation check, since when you stand on edges it fails //you can't use this kind of orientation check, since when you stand on edges it fails
if (fp == NULL || (checkOrientation != CoordType(0,0,0) && checkOrientation * fp->N() < 0.7)) if (fp == NULL || (checkOrientation != CoordType(0,0,0) && checkOrientation * fp->N() < 0.7))
@ -612,52 +604,54 @@ private:
v3 v3 v3 v3
Before Swap After Swap Before Swap After Swap
*/ */
static bool testSwap(PosType p, ScalarType creaseAngleCosThr)
//TODO: check if you can optimize using posType valence counting functions
static bool testSwap(const PosType & p, const ScalarType creaseAngleCosThr)
{ {
//if border or feature, do not swap //if border or feature, do not swap
if (/*p.IsBorder() || */p.IsEdgeS()) return false; if (/*p.IsBorder() || */p.IsEdgeS()) return false;
int oldDist = 0, newDist = 0, idealV, actualV; int oldDist = 0, newDist = 0, idealV = 0, actualV = 0;
PosType tp=p; PosType tp=p;
VertexType *v0=tp.V(); const VertexType *v0=tp.V();
std::vector<VertexType*> incident; std::vector<VertexType*> incident;
vcg::face::VVStarVF<FaceType>(tp.V(), incident); // vcg::face::VVStarVF<FaceType>(tp.V(), incident);
idealV = idealValence(tp); actualV = incident.size(); idealV = idealValence(tp); actualV = tp.NumberOfIncidentVertices();//int(incident.size());
oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV - 1)); oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV - 1));
tp.NextF();tp.FlipE();tp.FlipV(); tp.NextF();tp.FlipE();tp.FlipV();
VertexType *v1=tp.V(); const VertexType *v1=tp.V();
vcg::face::VVStarVF<FaceType>(tp.V(), incident); // vcg::face::VVStarVF<FaceType>(tp.V(), incident);
idealV = idealValence(tp); actualV = incident.size(); idealV = idealValence(tp); actualV = tp.NumberOfIncidentVertices();//int(incident.size());
oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV + 1)); oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV + 1));
tp.FlipE();tp.FlipV();tp.FlipE(); tp.FlipE();tp.FlipV();tp.FlipE();
VertexType *v2=tp.V(); const VertexType *v2=tp.V();
vcg::face::VVStarVF<FaceType>(tp.V(), incident); // vcg::face::VVStarVF<FaceType>(tp.V(), incident);
idealV = idealValence(tp); actualV = incident.size(); idealV = idealValence(tp); actualV = tp.NumberOfIncidentVertices();//int(incident.size());
oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV - 1)); oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV - 1));
tp.NextF();tp.FlipE();tp.FlipV(); tp.NextF();tp.FlipE();tp.FlipV();
VertexType *v3=tp.V(); const VertexType *v3=tp.V();
vcg::face::VVStarVF<FaceType>(tp.V(), incident); // vcg::face::VVStarVF<FaceType>(tp.V(), incident);
idealV = idealValence(tp); actualV = incident.size(); idealV = idealValence(tp); actualV = tp.NumberOfIncidentVertices();//int(incident.size());
oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV + 1)); oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV + 1));
ScalarType qOld = std::min(Quality(v0->P(),v2->P(),v3->P()),Quality(v0->P(),v1->P(),v2->P())); const ScalarType qOld = std::min(Quality(v0->P(),v2->P(),v3->P()),Quality(v0->P(),v1->P(),v2->P()));
ScalarType qNew = std::min(Quality(v0->P(),v1->P(),v3->P()),Quality(v2->P(),v3->P(),v1->P())); const ScalarType qNew = std::min(Quality(v0->P(),v1->P(),v3->P()),Quality(v2->P(),v3->P(),v1->P()));
return (newDist < oldDist && qNew >= qOld * 0.50f) || return (newDist < oldDist && qNew >= qOld * 0.50f) ||
(newDist == oldDist && qNew > qOld * 1.f) || qNew > 1.5f * qOld; (newDist == oldDist && qNew > qOld * 1.f) || qNew > 1.5f * qOld;
} }
static bool checkManifoldness(FaceType & f, int z) static bool checkManifoldness(const FaceType & f, const int z)
{ {
PosType pos(&f, (z+2)%3, f.V2(z)); PosType pos(&f, (z+2)%3, f.V2(z));
PosType start = pos; const PosType start = pos;
do { do {
pos.FlipE(); pos.FlipE();
@ -672,7 +666,7 @@ private:
// Edge swap step: edges are flipped in order to optimize valence and triangle quality across the mesh // Edge swap step: edges are flipped in order to optimize valence and triangle quality across the mesh
static void ImproveValence(MeshType &m, Params &params) static void ImproveValence(MeshType &m, Params &params)
{ {
static ScalarType foldCheckRad = math::ToRad(5.); const static ScalarType foldCheckRad = math::ToRad(5.);
tri::UpdateTopology<MeshType>::FaceFace(m); tri::UpdateTopology<MeshType>::FaceFace(m);
tri::UpdateTopology<MeshType>::VertexFace(m); tri::UpdateTopology<MeshType>::VertexFace(m);
ForEachFace(m, [&] (FaceType & f) { ForEachFace(m, [&] (FaceType & f) {
@ -681,10 +675,8 @@ private:
{ {
if (&f > f.cFFp(i)) if (&f > f.cFFp(i))
{ {
PosType pi(&f, i); const PosType pi(&f, i);
CoordType swapEdgeMidPoint = (f.cP2(i) + f.cFFp(i)->cP2(f.cFFi(i))) / 2.; const CoordType swapEdgeMidPoint = (f.cP2(i) + f.cFFp(i)->cP2(f.cFFi(i))) / 2.;
std::vector<CoordType> toCheck(1, swapEdgeMidPoint);
if(((!params.selectedOnly) || (f.IsS() && f.cFFp(i)->IsS())) && if(((!params.selectedOnly) || (f.IsS() && f.cFFp(i)->IsS())) &&
!face::IsBorder(f, i) && !face::IsBorder(f, i) &&
@ -692,15 +684,15 @@ private:
face::checkFlipEdgeNotManifold(f, i) && face::checkFlipEdgeNotManifold(f, i) &&
testSwap(pi, params.creaseAngleCosThr) && testSwap(pi, params.creaseAngleCosThr) &&
// face::CheckFlipEdge(f, i) && // face::CheckFlipEdge(f, i) &&
face::CheckFlipEdgeNormal(f, i, params.creaseAngleRadThr) && //vcg::math::ToRad(5.)) && face::CheckFlipEdgeNormal(f, i, float(vcg::math::ToRad(5.))) &&
(!params.surfDistCheck || testHausdorff(*params.mProject, params.grid, toCheck, params.maxSurfDist))) (!params.surfDistCheck || testHausdorff(*params.mProject, params.grid, { swapEdgeMidPoint }, params.maxSurfDist)))
{ {
//When doing the swap we need to preserve and update the crease info accordingly //When doing the swap we need to preserve and update the crease info accordingly
FaceType* g = f.cFFp(i); FaceType* g = f.cFFp(i);
int w = f.FFi(i); const int w = f.FFi(i);
bool creaseF = g->IsFaceEdgeS((w + 1) % 3); const bool creaseF = g->IsFaceEdgeS((w + 1) % 3);
bool creaseG = f.IsFaceEdgeS((i + 1) % 3); const bool creaseG = f.IsFaceEdgeS((i + 1) % 3);
face::FlipEdgeNotManifold(f, i); face::FlipEdgeNotManifold(f, i);
@ -732,9 +724,9 @@ private:
bool operator()(PosType &ep) bool operator()(PosType &ep)
{ {
ScalarType quality = (((math::Abs(ep.V()->Q())+math::Abs(ep.VFlip()->Q()))/(ScalarType)2.0)-minQ)/(maxQ-minQ); const ScalarType quality = ((ep.V()->Q()+ ep.VFlip()->Q())/(ScalarType)2.0);
ScalarType mult = computeLengthThrMult(params, quality); const ScalarType mult = computeLengthThrMult(params, quality);
ScalarType dist = Distance(ep.V()->P(), ep.VFlip()->P()); const ScalarType dist = Distance(ep.V()->P(), ep.VFlip()->P());
if(dist > mult * length) if(dist > mult * length)
{ {
++count; ++count;
@ -769,7 +761,7 @@ private:
tri::UpdateTopology<MeshType>::FaceFace(m); tri::UpdateTopology<MeshType>::FaceFace(m);
tri::MidPoint<MeshType> midFunctor(&m); tri::MidPoint<MeshType> midFunctor(&m);
ScalarType minQ,maxQ; ScalarType minQ = 0, maxQ = 0;
if(params.adapt){ if(params.adapt){
computeVQualityDistrMinMax(m, minQ, maxQ); computeVQualityDistrMinMax(m, minQ, maxQ);
EdgeSplitAdaptPred ep(params); EdgeSplitAdaptPred ep(params);
@ -790,20 +782,20 @@ private:
static int VtoE(const int v0, const int v1) static int VtoE(const int v0, const int v1)
{ {
static /*constexpr*/ int Vmat[3][3] = { -1, 0, 2, static constexpr int Vmat[3][3] = { -1, 0, 2,
0, -1, 1, 0, -1, 1,
2, 1, -1}; 2, 1, -1};
return Vmat[v0][v1]; return Vmat[v0][v1];
} }
static bool checkCanMoveOnCollapse(PosType p, std::vector<FaceType*> & faces, std::vector<int> & vIdxes, Params &params) static bool checkCanMoveOnCollapse(const PosType & p, const std::vector<FaceType*> & faces, const std::vector<int> & vIdxes, const Params &params)
{ {
bool allIncidentFaceSelected = true; bool allIncidentFaceSelected = true;
PosType pi = p; PosType pi = p;
CoordType dEdgeVector = (p.V()->cP() - p.VFlip()->cP()).Normalize(); const CoordType dEdgeVector = (p.V()->cP() - p.VFlip()->cP()).Normalize();
int incidentFeatures = 0; int incidentFeatures = 0;
@ -815,7 +807,7 @@ private:
{ {
vcg::tri::Mark(*params.m,faces[i]->V1(vIdxes[i])); vcg::tri::Mark(*params.m,faces[i]->V1(vIdxes[i]));
incidentFeatures++; incidentFeatures++;
CoordType movingEdgeVector0 = (faces[i]->cP1(vIdxes[i]) - faces[i]->cP(vIdxes[i])).Normalize(); const CoordType movingEdgeVector0 = (faces[i]->cP1(vIdxes[i]) - faces[i]->cP(vIdxes[i])).Normalize();
if (std::fabs(movingEdgeVector0 * dEdgeVector) < .9f || !p.IsEdgeS()) if (std::fabs(movingEdgeVector0 * dEdgeVector) < .9f || !p.IsEdgeS())
return false; return false;
} }
@ -823,7 +815,7 @@ private:
{ {
vcg::tri::Mark(*params.m,faces[i]->V2(vIdxes[i])); vcg::tri::Mark(*params.m,faces[i]->V2(vIdxes[i]));
incidentFeatures++; incidentFeatures++;
CoordType movingEdgeVector1 = (faces[i]->cP2(vIdxes[i]) - faces[i]->cP(vIdxes[i])).Normalize(); const CoordType movingEdgeVector1 = (faces[i]->cP2(vIdxes[i]) - faces[i]->cP(vIdxes[i])).Normalize();
if (std::fabs(movingEdgeVector1 * dEdgeVector) < .9f || !p.IsEdgeS()) if (std::fabs(movingEdgeVector1 * dEdgeVector) < .9f || !p.IsEdgeS())
return false; return false;
} }
@ -836,61 +828,51 @@ private:
return params.selectedOnly ? allIncidentFaceSelected : true; return params.selectedOnly ? allIncidentFaceSelected : true;
} }
static bool checkFacesAfterCollapse (std::vector<FaceType*> & faces, PosType p, const Point3<ScalarType> &mp, Params &params, bool relaxed) static bool checkFacesAfterCollapse (const std::vector<FaceType*> & faces, const PosType & p, const Point3<ScalarType> &mp, Params &params, bool relaxed)
{ {
for (FaceType* f : faces) for (FaceType* f : faces)
{ {
if(!(*f).IsD() && f != p.F()) //i'm not a deleted face if(!(*f).IsD() && f != p.F()) //i'm not a deleted face
{ {
PosType pi(f, p.V()); //same vertex const PosType pi(f, p.V()); //same vertex
VertexType *v0 = pi.V(); const auto v0 = pi.V();
VertexType *v1 = pi.F()->V1(pi.VInd()); const auto v1 = pi.F()->V1(pi.VInd());
VertexType *v2 = pi.F()->V2(pi.VInd()); const auto v2 = pi.F()->V2(pi.VInd());
if( v1 == p.VFlip() || v2 == p.VFlip()) //i'm the other deleted face if( v1 == p.VFlip() || v2 == p.VFlip()) //i'm the other deleted face
continue; continue;
//check on new face quality //check on new face quality
{ {
ScalarType newQ = Quality(mp, v1->P(), v2->P()); const auto newQ = Quality(mp, v1->P(), v2->P());
ScalarType oldQ = Quality(v0->P(), v1->P(), v2->P()); const auto oldQ = Quality(v0->P(), v1->P(), v2->P());
if(newQ <= 0.5*oldQ) if(newQ <= 0.5*oldQ)
return false; return false;
} }
// we prevent collapse that makes edges too long (except for cross) // we prevent collapse that makes edges too long (except for cross)
if(!relaxed) if(!relaxed && (Distance(mp, v1->P()) > params.maxLength || Distance(mp, v2->P()) > params.maxLength))
if((Distance(mp, v1->P()) > params.maxLength || Distance(mp, v2->P()) > params.maxLength))
return false; return false;
Point3<ScalarType> oldN = NormalizedTriangleNormal(*(pi.F())); const auto oldN = NormalizedTriangleNormal(*(pi.F()));
Point3<ScalarType> newN = Normal(mp, v1->P(), v2->P()).Normalize(); const auto newN = Normal(mp, v1->P(), v2->P()).Normalize();
// if (oldN * newN < 0.5f) if (oldN * newN < 0.7f)
// return false;
std::vector<CoordType> baryP(1);
baryP[0] = (v1->cP() + v2->cP() + mp) / 3.;
if (!testHausdorff(*(params.mProject), params.grid, baryP, params.maxSurfDist, newN))
return false; return false;
//check on new face distance from original mesh //check on new face distance from original mesh
if (params.surfDistCheck) if (params.surfDistCheck)
{ {
std::vector<CoordType> points(3); const auto points = {
std::vector<CoordType> baryP(1); (v1->cP() + mp) / 2.,
(v2->cP() + mp) / 2.,
mp,
};
baryP[0] = (v1->cP() + v2->cP() + mp) / 3.; if (!testHausdorff(*(params.mProject), params.grid, points, params.maxSurfDist) ||
!testHausdorff(*(params.mProject), params.grid, { (v1->cP() + v2->cP() + mp) / 3. }, params.maxSurfDist, newN))
points[0] = (v1->cP() + mp) / 2.;
points[1] = (v2->cP() + mp) / 2.;
points[2] = mp;
if (!testHausdorff(*(params.mProject), params.grid, points, params.maxSurfDist))// ||
// !testHausdorff(*(params.mProject), params.grid, baryP, params.maxSurfDist, newN))
return false; return false;
} }
} }
@ -900,10 +882,9 @@ private:
//TODO: Refactor code and implement the correct set up of crease info when collapsing towards a crease edge //TODO: Refactor code and implement the correct set up of crease info when collapsing towards a crease edge
static bool checkCollapseFacesAroundVert1(PosType &p, VertexPair & pair, Point3<ScalarType> &mp, Params &params, bool relaxed) static bool checkCollapseFacesAroundVert1(const PosType &p, VertexPair & pair, Point3<ScalarType> &mp, Params &params, bool relaxed)
{ {
PosType p0 = p, p1 = p; PosType p0 = p, p1 = p;
p1.FlipV(); p1.FlipV();
std::vector<int> vi0, vi1; std::vector<int> vi0, vi1;
@ -913,8 +894,8 @@ private:
face::VFStarVF<FaceType>(p1.V(), ff1, vi1); face::VFStarVF<FaceType>(p1.V(), ff1, vi1);
//check crease-moveability //check crease-moveability
bool moveable0 = checkCanMoveOnCollapse(p0, ff0, vi0, params) && !p0.V()->IsS(); const bool moveable0 = checkCanMoveOnCollapse(p0, ff0, vi0, params) && !p0.V()->IsS();
bool moveable1 = checkCanMoveOnCollapse(p1, ff1, vi1, params) && !p1.V()->IsS(); const bool moveable1 = checkCanMoveOnCollapse(p1, ff1, vi1, params) && !p1.V()->IsS();
//if both moveable => go to midpoint //if both moveable => go to midpoint
// else collapse on movable one // else collapse on movable one
@ -923,9 +904,6 @@ private:
pair = moveable0 ? VertexPair(p0.V(), p1.V()) : VertexPair(p1.V(), p0.V()); pair = moveable0 ? VertexPair(p0.V(), p1.V()) : VertexPair(p1.V(), p0.V());
//casting int(true) is always 1 and int(false) = =0
assert(int(true) == 1);
assert(int(false) == 0);
mp = (p0.V()->cP() * int(moveable1) + p1.V()->cP() * int(moveable0)) / (int(moveable0) + int(moveable1)); mp = (p0.V()->cP() * int(moveable1) + p1.V()->cP() * int(moveable0)) / (int(moveable0) + int(moveable1));
if (checkFacesAfterCollapse(ff0, p0, mp, params, relaxed)) if (checkFacesAfterCollapse(ff0, p0, mp, params, relaxed))
@ -934,14 +912,15 @@ private:
return false; return false;
} }
static bool testCollapse1(PosType &p, VertexPair & pair, Point3<ScalarType> &mp, ScalarType minQ, ScalarType maxQ, Params &params, bool relaxed = false) static bool testCollapse1(const PosType &p, VertexPair & pair, Point3<ScalarType> &mp, ScalarType minQ, ScalarType maxQ, Params &params, bool relaxed = false)
{ {
ScalarType quality = (((math::Abs(p.V()->Q())+math::Abs(p.VFlip()->Q()))/(ScalarType)2.0)-minQ)/(maxQ-minQ); const ScalarType quality = params.adapt ? ((p.V()->Q()+ p.VFlip()->Q())/(ScalarType)2.0) : 0;
ScalarType mult = computeLengthThrMult(params, quality);
ScalarType thr = mult*params.minLength;
ScalarType dist = Distance(p.V()->P(), p.VFlip()->P()); const ScalarType mult = computeLengthThrMult(params, quality);
ScalarType area = DoubleArea(*(p.F()))/2.f; const ScalarType thr = mult*params.minLength;
const ScalarType dist = Distance(p.V()->P(), p.VFlip()->P());
const ScalarType area = DoubleArea(*(p.F()))/2.f;
if(relaxed || (dist < thr || area < params.minLength*params.minLength/100.f))//if to collapse if(relaxed || (dist < thr || area < params.minLength*params.minLength/100.f))//if to collapse
{ {
return checkCollapseFacesAroundVert1(p, pair, mp, params, relaxed); return checkCollapseFacesAroundVert1(p, pair, mp, params, relaxed);
@ -969,9 +948,9 @@ private:
if(!(*v).IsD() && (*v).IsB() && v != p.V()) //ignore non border if(!(*v).IsD() && (*v).IsB() && v != p.V()) //ignore non border
collapsedNV1 = ((*v).P() - p.V()->P()).normalized(); //edge vector after collapse collapsedNV1 = ((*v).P() - p.V()->P()).normalized(); //edge vector after collapse
float cosine = cos(math::ToRad(1.5f)); const float cosine = cos(math::ToRad(1.5f));
float angle0 = fabs(fastAngle(collapseNV, collapsedNV0)); const float angle0 = fabs(fastAngle(collapseNV, collapsedNV0));
float angle1 = fabs(fastAngle(collapseNV, collapsedNV1)); const float angle1 = fabs(fastAngle(collapseNV, collapsedNV1));
//if on both sides we deviate too much after collapse => don't collapse //if on both sides we deviate too much after collapse => don't collapse
if(angle0 <= cosine && angle1 <= cosine) if(angle0 <= cosine && angle1 <= cosine)
return false; return false;
@ -984,7 +963,7 @@ private:
// the linkConditions are preserved // the linkConditions are preserved
static void CollapseShortEdges(MeshType &m, Params &params) static void CollapseShortEdges(MeshType &m, Params &params)
{ {
ScalarType minQ, maxQ; ScalarType minQ = 0, maxQ = 0;
int candidates = 0; int candidates = 0;
if(params.adapt) if(params.adapt)
@ -998,17 +977,16 @@ private:
ss.push(); ss.push();
{ {
tri::UpdateTopology<MeshType>::FaceFace(m); // tri::UpdateTopology<MeshType>::FaceFace(m);
Clean<MeshType>::CountNonManifoldVertexFF(m,true); Clean<MeshType>::CountNonManifoldVertexFF(m,true);
//FROM NOW ON VSelection is NotManifold //FROM NOW ON VSelection is NotManifold
ForEachFace(m, [&](FaceType &f) {
for(auto fi=m.face.begin(); fi!=m.face.end(); ++fi) if(!f.IsD() && (params.selectedOnly == false || f.IsS()))
if(!(*fi).IsD() && (params.selectedOnly == false || fi->IsS()))
{ {
for(auto i=0; i<3; ++i) for(auto i=0; i<3; ++i)
{ {
PosType pi(&*fi, i); PosType pi(&f, i);
++candidates; ++candidates;
VertexPair bp = VertexPair(pi.V(), pi.VFlip()); VertexPair bp = VertexPair(pi.V(), pi.VFlip());
Point3<ScalarType> mp = (pi.V()->P()+pi.VFlip()->P())/2.f; Point3<ScalarType> mp = (pi.V()->P()+pi.VFlip()->P())/2.f;
@ -1022,6 +1000,7 @@ private:
} }
} }
});
} }
ss.pop(); ss.pop();
} }
@ -1103,52 +1082,8 @@ private:
// v2 = (fv1 == v1) ? fv2 : fv1; // v2 = (fv1 == v1) ? fv2 : fv1;
} }
} }
face::VVStarVF<FaceType>(v0, vv0);
face::VVStarVF<FaceType>(v1, vv1);
face::VVStarVF<FaceType>(v2, vv2);
face::VVStarVF<FaceType>(v3, vv3);
int nv0 = vv0.size(), nv1 = vv1.size();
int nv2 = vv2.size(), nv3 = vv3.size();
int delta1 = (idealValence(*v0) - nv0) + (idealValence(*v2) - nv2);
int delta2 = (idealValence(*v1) - nv1) + (idealValence(*v3) - nv3);
ScalarType Q1 = std::min(Quality(v0->P(), v1->P(), v3->P()), Quality(v1->P(), v2->P(), v3->P()));
ScalarType Q2 = std::min(Quality(v0->P(), v1->P(), v2->P()), Quality(v2->P(), v3->P(), v0->P()));
if (crease[0] || crease[1] || crease[2] || crease[3])
return false;
// if (crease[0] && crease[1] && crease[2] && crease[3])
// {
// return false;
// }
// if (crease[0] || crease[2])
// {
// bp = VertexPair(p.V(), v0);
// return true;
// }
// if (crease[1] || crease[3])
// {
// bp = VertexPair(p.V(), v1);
// return true;
// }
//no crease
if(delta1 < delta2 && Q1 >= 0.6f*Q2)
{
bp = VertexPair(p.V(), v1);
return true;
}
else
{
bp = VertexPair(p.V(), v0);
return true;
}
} }
//Cross Collapse pass: This pass cleans the mesh from cross vertices, keeping in mind the link conditions //Cross Collapse pass: This pass cleans the mesh from cross vertices, keeping in mind the link conditions
//and feature preservations tests. //and feature preservations tests.
static void CollapseCrosses(MeshType &m , Params &params) static void CollapseCrosses(MeshType &m , Params &params)
@ -1160,19 +1095,17 @@ private:
SelectionStack<MeshType> ss(m); SelectionStack<MeshType> ss(m);
ss.push(); ss.push();
{ {
tri::UpdateTopology<MeshType>::FaceFace(m); tri::UpdateTopology<MeshType>::FaceFace(m);
Clean<MeshType>::CountNonManifoldVertexFF(m,true); Clean<MeshType>::CountNonManifoldVertexFF(m,true);
//From now on Selection on vertices is not manifoldness //From now on Selection on vertices is not manifoldness
ForEachFace(m, [&](FaceType &f) {
for(auto fi=m.face.begin(); fi!=m.face.end(); ++fi) if(!f.IsD() && (params.selectedOnly == false || f.IsS()))
if(!(*fi).IsD() && (!params.selectedOnly || fi->IsS()))
{ {
for(auto i=0; i<3; ++i) for(auto i=0; i<3; ++i)
{ {
PosType pi(&*fi, i); PosType pi(&f, i);
if(!pi.V()->IsB()) if(!pi.V()->IsB())
{ {
std::vector<FaceType*> ff; std::vector<FaceType*> ff;
@ -1198,6 +1131,7 @@ private:
} }
} }
} }
});
} }
ss.pop(); ss.pop();