1 /**************************************************************************** 2 * VCGLib o o * 3 * Visual and Computer Graphics Library o o * 4 * _ O _ * 5 * Copyright(C) 2004-2017 \/)\/ * 6 * Visual Computing Lab /\/| * 7 * ISTI - Italian National Research Council | * 8 * \ * 9 * All rights reserved. * 10 * * 11 * This program is free software; you can redistribute it and/or modify * 12 * it under the terms of the GNU General Public License as published by * 13 * the Free Software Foundation; either version 2 of the License, or * 14 * (at your option) any later version. * 15 * * 16 * This program is distributed in the hope that it will be useful, * 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * 20 * for more details. * 21 * * 22 ****************************************************************************/ 23 #ifndef _VCG_ISOTROPICREMESHING_H 24 #define _VCG_ISOTROPICREMESHING_H 25 26 #include <vcg/complex/algorithms/update/quality.h> 27 #include <vcg/complex/algorithms/update/curvature.h> 28 #include <vcg/complex/algorithms/update/normal.h> 29 #include <vcg/complex/algorithms/refine.h> 30 #include <vcg/complex/algorithms/stat.h> 31 #include <vcg/complex/algorithms/smooth.h> 32 #include <vcg/complex/algorithms/local_optimization/tri_edge_collapse.h> 33 #include <vcg/space/index/spatial_hashing.h> 34 #include <vcg/complex/append.h> 35 #include <vcg/complex/allocate.h> 36 #include <wrap/io_trimesh/export.h> 37 38 namespace vcg { 39 namespace tri { 40 template<class TRI_MESH_TYPE> 41 class IsotropicRemeshing 42 { 43 public: 44 typedef TRI_MESH_TYPE MeshType; 45 typedef typename MeshType::FaceType FaceType; 46 typedef typename MeshType::FacePointer FacePointer; 47 typedef typename FaceType::VertexType VertexType; 48 typedef typename FaceType::VertexPointer VertexPointer; 49 typedef typename VertexType::ScalarType ScalarType; 50 typedef typename VertexType::CoordType CoordType; 51 typedef typename face::Pos<FaceType> PosType; 52 typedef BasicVertexPair<VertexType> VertexPair; 53 typedef EdgeCollapser<MeshType, VertexPair> Collapser; 54 typedef GridStaticPtr<FaceType, ScalarType> StaticGrid; 55 56 57 typedef struct Params { 58 typedef struct Stat { 59 int splitNum; 60 int collapseNum; 61 int flipNum; 62 ResetParams::Stat63 void Reset() { 64 splitNum=0; 65 collapseNum=0; 66 flipNum=0; 67 } 68 } Stat; 69 70 71 ScalarType minLength; // minimal admitted length: no edge should be shorter than this value (used when collapsing) 72 ScalarType maxLength; // maximal admitted length: no edge should be longer than this value (used when refining) 73 ScalarType lengthThr; 74 75 ScalarType minimalAdmittedArea; 76 ScalarType maxSurfDist; 77 78 ScalarType aspectRatioThr = 0.05; //min aspect ratio: during relax bad triangles will be relaxed 79 ScalarType foldAngleCosThr = cos(math::ToRad(140.)); //min angle to be considered folding: during relax folded triangles will be relaxed 80 81 ScalarType creaseAngleRadThr = math::ToRad(10.0); 82 ScalarType creaseAngleCosThr = cos(math::ToRad(10.0)); //min angle to be considered crease: two faces with normals diverging more than thr share a crease edge 83 84 bool splitFlag = true; 85 bool swapFlag = true; 86 bool collapseFlag = true; 87 bool smoothFlag=true; 88 bool projectFlag=true; 89 bool selectedOnly = false; 90 bool cleanFlag = true; 91 92 bool userSelectedCreases = false; 93 bool surfDistCheck = true; 94 95 bool adapt=false; 96 int iter=1; 97 Stat stat; 98 SetTargetLenParams99 void SetTargetLen(const ScalarType len) 100 { 101 minLength=len*4./5.; 102 maxLength=len*4./3.; 103 lengthThr=len*4./3.; 104 minimalAdmittedArea = (minLength * minLength)/1000.0; 105 } 106 SetFeatureAngleDegParams107 void SetFeatureAngleDeg(const ScalarType angle) 108 { 109 creaseAngleRadThr = math::ToRad(angle); 110 creaseAngleCosThr = cos(creaseAngleRadThr); 111 } 112 113 StaticGrid grid; 114 MeshType* m; 115 MeshType* mProject; 116 117 } Params; 118 119 private: debug_crease(MeshType & toRemesh,std::string prepend,int i)120 static void debug_crease (MeshType & toRemesh, std::string prepend, int i) 121 { 122 ForEachVertex(toRemesh, [] (VertexType & v) { 123 v.C() = Color4b::Gray; 124 v.Q() = 0; 125 }); 126 127 ForEachFacePos(toRemesh, [&](PosType &p){ 128 if (p.F()->IsFaceEdgeS(p.E())) 129 { 130 p.V()->Q() += 1; 131 p.VFlip()->Q() += 1; 132 } 133 }); 134 135 ForEachVertex(toRemesh, [] (VertexType & v) { 136 if (v.Q() >= 4) 137 v.C() = Color4b::Green; 138 else if (v.Q() >= 2) 139 v.C() = Color4b::Red; 140 }); 141 prepend += "_creases" + std::to_string(i) + ".ply"; 142 vcg::tri::io::Exporter<MeshType>::Save(toRemesh, prepend.c_str(), vcg::tri::io::Mask::IOM_ALL); 143 } 144 145 removeColinearFaces(MeshType & m,Params & params)146 static void removeColinearFaces(MeshType & m, Params & params) 147 { 148 vcg::tri::UpdateTopology<MeshType>::FaceFace(m); 149 150 int count = 0; 151 int iter = 0; 152 do 153 { 154 vcg::tri::UpdateTopology<MeshType>::FaceFace(m); 155 vcg::tri::UnMarkAll(m); 156 157 count = 0; 158 for (size_t i = 0; i < size_t(m.FN()); ++i) 159 { 160 FaceType & f = m.face[i]; 161 162 ScalarType quality = vcg::QualityRadii(f.cP(0), f.cP(1), f.cP(2)); 163 164 if (quality <= 0.001) 165 { 166 //find longest edge 167 double edges[3]; 168 edges[0] = vcg::Distance(f.cP(0), f.cP(1)); 169 edges[1] = vcg::Distance(f.cP(1), f.cP(2)); 170 edges[2] = vcg::Distance(f.cP(2), f.cP(0)); 171 172 ScalarType smallestEdge = std::min(edges[0], std::min(edges[1], edges[2])); 173 int longestIdx = std::find(edges, edges+3, std::max(std::max(edges[0], edges[1]), edges[2])) - (edges); 174 175 if (vcg::tri::IsMarked(m, f.V2(longestIdx))) 176 continue; 177 178 179 auto f1 = f.cFFp(longestIdx); 180 vcg::tri::Mark(m,f.V2(longestIdx)); 181 if (!vcg::face::IsBorder(f, longestIdx) && vcg::face::IsManifold(f, longestIdx) && vcg::face::checkFlipEdgeNotManifold<FaceType>(f, longestIdx)) { 182 183 // Check if EdgeFlipping improves quality 184 FacePointer g = f.FFp(longestIdx); int k = f.FFi(longestIdx); 185 vcg::Triangle3<ScalarType> t1(f.P(longestIdx), f.P1(longestIdx), f.P2(longestIdx)), t2(g->P(k), g->P1(k), g->P2(k)), 186 t3(f.P(longestIdx), g->P2(k), f.P2(longestIdx)), t4(g->P(k), f.P2(longestIdx), g->P2(k)); 187 188 if ( std::min( QualityFace(t1), QualityFace(t2) ) <= std::min( QualityFace(t3), QualityFace(t4) )) 189 { 190 ScalarType dist; 191 CoordType closest; 192 auto fp0 = vcg::tri::GetClosestFaceBase(*params.mProject, params.grid, vcg::Barycenter(t3), smallestEdge/4., dist, closest); 193 if (fp0 == NULL) 194 continue; 195 196 auto fp1 = vcg::tri::GetClosestFaceBase(*params.mProject, params.grid, vcg::Barycenter(t4), smallestEdge/4., dist, closest); 197 if (fp1 == NULL) 198 continue; 199 200 vcg::face::FlipEdgeNotManifold<FaceType>(f, longestIdx); 201 ++count; 202 } 203 } 204 } 205 } 206 } while (count && ++iter < 50); 207 } 208 cleanMesh(MeshType & m,Params & params)209 static void cleanMesh(MeshType & m, Params & params) 210 { 211 vcg::tri::Clean<MeshType>::RemoveDuplicateFace(m); 212 vcg::tri::Clean<MeshType>::RemoveUnreferencedVertex(m); 213 vcg::tri::Allocator<MeshType>::CompactEveryVector(m); 214 215 vcg::tri::UpdateTopology<MeshType>::FaceFace(m); 216 removeColinearFaces(m, params); 217 vcg::tri::UpdateTopology<MeshType>::FaceFace(m); 218 } 219 220 public: 221 222 static void Do(MeshType &toRemesh, Params & params, vcg::CallBackPos * cb=0) 223 { 224 MeshType toProjectCopy; 225 tri::UpdateBounding<MeshType>::Box(toRemesh); 226 tri::UpdateNormal<MeshType>::PerVertexNormalizedPerFaceNormalized(toRemesh); 227 228 tri::Append<MeshType,MeshType>::MeshCopy(toProjectCopy, toRemesh); 229 230 Do(toRemesh,toProjectCopy,params,cb); 231 } 232 static void Do(MeshType &toRemesh, MeshType &toProject, Params & params, vcg::CallBackPos * cb=0) 233 { 234 assert(&toRemesh != &toProject); 235 params.stat.Reset(); 236 237 238 tri::UpdateBounding<MeshType>::Box(toRemesh); 239 240 { 241 tri::UpdateBounding<MeshType>::Box(toProject); 242 tri::UpdateNormal<MeshType>::PerFaceNormalized(toProject); 243 params.m = &toRemesh; 244 params.mProject = &toProject; 245 params.grid.Set(toProject.face.begin(), toProject.face.end()); 246 } 247 248 if (params.cleanFlag) 249 cleanMesh(toRemesh, params); 250 251 tri::UpdateTopology<MeshType>::FaceFace(toRemesh); 252 tri::UpdateFlags<MeshType>::VertexBorderFromFaceAdj(toRemesh); 253 tri::UpdateTopology<MeshType>::VertexFace(toRemesh); 254 255 // computeQuality(toRemesh); 256 // tri::UpdateQuality<MeshType>::VertexSaturate(toRemesh); 257 258 tagCreaseEdges(toRemesh, params); 259 260 for(int i=0; i < params.iter; ++i) 261 { 262 // params.stat.Reset(); 263 if(cb) cb(100*i/params.iter, "Remeshing"); 264 265 if(params.splitFlag) 266 SplitLongEdges(toRemesh, params); 267 #ifdef DEBUG_CREASE 268 debug_crease(toRemesh, std::string("after_ref"), i); 269 #endif 270 271 if(params.collapseFlag) 272 { 273 CollapseShortEdges(toRemesh, params); 274 CollapseCrosses(toRemesh, params); 275 } 276 277 if(params.swapFlag) 278 ImproveValence(toRemesh, params); 279 280 if(params.smoothFlag) 281 ImproveByLaplacian(toRemesh, params); 282 283 if(params.projectFlag) 284 ProjectToSurface(toRemesh, params); 285 } 286 } 287 288 private: 289 /* 290 TODO: Add better crease support: detect all creases at starting time, saving it on facedgesel flags 291 All operations must then preserve the faceedgesel flag accordingly: 292 Refinement -> Check that refiner propagates faceedgesel [should be doing it] 293 Collapse -> Implement 1D edge collapse and better check on corners and creases 294 Swap -> Totally avoid swapping crease edges [ok] 295 Smooth -> Apply 1D smoothing to crease vertices + check on 296 (http://www.cs.ubc.ca/labs/imager/tr/2009/eltopo/sisc2009.pdf) 297 */ IsotropicRemeshing()298 IsotropicRemeshing() {} 299 // this returns the value of cos(a) where a is the angle between n0 and n1. (scalar prod is cos(a)) fastAngle(Point3<ScalarType> n0,Point3<ScalarType> n1)300 static inline ScalarType fastAngle(Point3<ScalarType> n0, Point3<ScalarType> n1) 301 { 302 return math::Clamp(n0*n1,(ScalarType)-1.0,(ScalarType)1.0); 303 } 304 // compare the value of the scalar prod with the cos of the crease threshold testCreaseEdge(PosType & p,ScalarType creaseCosineThr)305 static inline bool testCreaseEdge(PosType &p, ScalarType creaseCosineThr) 306 { 307 ScalarType angle = fastAngle(NormalizedTriangleNormal(*(p.F())), NormalizedTriangleNormal(*(p.FFlip()))); 308 return angle <= creaseCosineThr && angle >= -0.98; 309 // return (angle <= creaseCosineThr && angle >= -creaseCosineThr); 310 } 311 // this stores in minQ the value of the 10th percentile of the VertQuality distribution and in 312 // maxQ the value of the 90th percentile. computeVQualityDistrMinMax(MeshType & m,ScalarType & minQ,ScalarType & maxQ)313 static inline void computeVQualityDistrMinMax(MeshType &m, ScalarType &minQ, ScalarType &maxQ) 314 { 315 Distribution<ScalarType> distr; 316 tri::Stat<MeshType>::ComputePerVertexQualityDistribution(m,distr); 317 318 maxQ = distr.Percentile(0.9f); 319 minQ = distr.Percentile(0.1f); 320 } 321 322 //Computes PerVertexQuality as a function of the 'deviation' of the normals taken from 323 //the faces incident to each vertex computeQuality(MeshType & m)324 static void computeQuality(MeshType &m) 325 { 326 tri::RequirePerVertexQuality(m); 327 tri::UpdateFlags<MeshType>::VertexClearV(m); 328 329 for(auto vi=m.vert.begin(); vi!=m.vert.end(); ++vi) 330 if(!(*vi).IsD()) 331 { 332 vector<FaceType*> ff; 333 face::VFExtendedStarVF(&*vi, 0, ff); 334 335 ScalarType tot = 0.f; 336 auto it = ff.begin(); 337 Point3<ScalarType> fNormal = NormalizedTriangleNormal(**it); 338 ++it; 339 while(it != ff.end()) 340 { 341 tot+= 1-math::Abs(fastAngle(fNormal, NormalizedTriangleNormal(**it))); 342 ++it; 343 } 344 vi->Q() = tot / (ScalarType)(std::max(1, ((int)ff.size()-1))); 345 vi->SetV(); 346 } 347 } 348 349 /* 350 Computes the ideal valence for the vertex in pos p: 351 4 for border vertices 352 6 for internal vertices 353 */ idealValence(PosType & p)354 static inline int idealValence(PosType &p) 355 { 356 if(p.IsBorder()) return 4; 357 return 6; 358 } idealValence(VertexType & v)359 static inline int idealValence(VertexType &v) 360 { 361 if(v.IsB()) return 4; 362 return 6; 363 } idealValenceSlow(PosType & p)364 static inline int idealValenceSlow(PosType &p) 365 { 366 std::vector<PosType> posVec; 367 VFOrderedStarFF(p,posVec); 368 float angleSumRad =0; 369 for(PosType &ip : posVec) 370 { 371 angleSumRad += ip.AngleRad(); 372 } 373 374 return (int)(std::ceil(angleSumRad / (M_PI/3.0f))); 375 } 376 testHausdorff(MeshType & m,StaticGrid & grid,const std::vector<CoordType> & verts,const ScalarType maxD)377 static bool testHausdorff (MeshType & m, StaticGrid & grid, const std::vector<CoordType> & verts, const ScalarType maxD) 378 { 379 for (CoordType v : verts) 380 { 381 CoordType closest; 382 ScalarType dist = 0; 383 FaceType* fp = GetClosestFaceBase(m, grid, v, maxD, dist, closest); 384 385 if (fp == NULL) 386 { 387 return false; 388 } 389 } 390 return true; 391 } 392 tagCreaseEdges(MeshType & m,Params & params)393 static int tagCreaseEdges(MeshType &m, Params & params) 394 { 395 int count = 0; 396 std::vector<char> creaseVerts(m.VN(), 0); 397 398 vcg::tri::UpdateFlags<MeshType>::VertexClearV(m); 399 std::queue<PosType> creaseQueue; 400 ForEachFacePos(m, [&](PosType &p){ 401 402 if (p.IsBorder()) 403 p.F()->SetFaceEdgeS(p.E()); 404 405 // if((p.FFlip() > p.F())) 406 { 407 FaceType *ff = p.F(); 408 FaceType *ffAdj = p.FFlip(); 409 410 double quality = vcg::QualityRadii(ff->cP(0), ff->cP(1), ff->cP(2)); 411 double qualityAdj = vcg::QualityRadii(ffAdj->cP(0), ffAdj->cP(1), ffAdj->cP(2)); 412 413 bool qualityCheck = quality > 0.00000001 && qualityAdj > 0.00000001; 414 // bool areaCheck = vcg::DoubleArea(*ff) > 0.000001 && vcg::DoubleArea(*ffAdj) > 0.000001; 415 416 if (!params.userSelectedCreases && (testCreaseEdge(p, params.creaseAngleCosThr) /*&& areaCheck*//* && qualityCheck*/) || p.IsBorder()) 417 { 418 PosType pp = p; 419 do { 420 pp.F()->SetFaceEdgeS(pp.E()); 421 pp.NextF(); 422 } while (pp != p); 423 424 creaseQueue.push(p); 425 } 426 } 427 }); 428 429 // //now all creases are checked... 430 // //prune false positive (too small) (count + scale?) 431 432 // while (!creaseQueue.empty()) 433 // { 434 // PosType & p = creaseQueue.front(); 435 // creaseQueue.pop(); 436 437 // std::stack<PosType> chainQueue; 438 // std::vector<size_t> chainVerts; 439 440 // if (!p.V()->IsV()) 441 // { 442 // chainQueue.push(p); 443 // } 444 445 // p.FlipV(); 446 // p.NextEdgeS(); 447 448 // if (!p.V()->IsV()) 449 // { 450 // chainQueue.push(p); 451 // } 452 453 // while (!chainQueue.empty()) 454 // { 455 // PosType p = chainQueue.top(); 456 // chainQueue.pop(); 457 458 // p.V()->SetV(); 459 // chainVerts.push_back(vcg::tri::Index(m, p.V())); 460 461 // PosType pp = p; 462 463 // //circle around vert in search for new crease edges 464 // do { 465 // pp.NextF(); //jump adj face 466 // pp.FlipE(); // this edge is already ok => jump to next 467 // if (pp.IsEdgeS()) 468 // { 469 // PosType nextPos = pp; 470 // nextPos.FlipV(); // go to next vert in the chain 471 // if (!nextPos.V()->IsV()) // if already visited...ignore 472 // { 473 // chainQueue.push(nextPos); 474 // } 475 // } 476 // } 477 // while (pp != p); 478 479 // } 480 481 // if (chainVerts.size() > 5) 482 // { 483 // for (auto vp : chainVerts) 484 // { 485 // creaseVerts[vp] = 1; 486 // } 487 // } 488 // } 489 // //store crease on V() 490 491 // //this aspect ratio check doesn't work on cadish meshes (long thin triangles spanning whole mesh) 492 // ForEachFace(m, [&] (FaceType & f) { 493 // if (vcg::QualityRadii(f.cP(0), f.cP(1), f.cP(2)) < params.aspectRatioThr) 494 // { 495 // if (creaseVerts[vcg::tri::Index(m, f.V(0))] == 0) 496 // f.V(0)->SetS(); 497 // if (creaseVerts[vcg::tri::Index(m, f.V(1))] == 0) 498 // f.V(1)->SetS(); 499 // if (creaseVerts[vcg::tri::Index(m, f.V(2))] == 0) 500 // f.V(2)->SetS(); 501 // } 502 // }); 503 504 // ForEachFace(m, [&] (FaceType & f) { 505 // for (int i = 0; i < 3; ++i) 506 // { 507 // if (f.FFp(i) > &f) 508 // { 509 // ScalarType angle = fastAngle(NormalizedTriangleNormal(f), NormalizedTriangleNormal(*(f.FFp(i)))); 510 // if (angle <= params.foldAngleCosThr) 511 // { 512 // // if (creaseVerts[vcg::tri::Index(m, f.V0(i))] == 0) 513 // f.V0(i)->SetS(); 514 // // if (creaseVerts[vcg::tri::Index(m, f.V1(i))] == 0) 515 // f.V1(i)->SetS(); 516 // // if (creaseVerts[vcg::tri::Index(m, f.V2(i))] == 0) 517 // f.V2(i)->SetS(); 518 // // if (creaseVerts[vcg::tri::Index(m, f.FFp(i)->V2(f.FFi(i)))] == 0) 519 // f.FFp(i)->V2(f.FFi(i))->SetS(); 520 // } 521 // } 522 // } 523 // }); 524 525 return count; 526 } 527 528 529 /* 530 Edge Swap Step: 531 This method optimizes the valence of each vertex. 532 oldDist is the sum of the absolute distance of each vertex from its ideal valence 533 newDist is the sum of the absolute distance of each vertex from its ideal valence after 534 the edge swap. 535 If the swap decreases the total absolute distance, then it's applied, preserving the triangle 536 quality. +1 537 v1 v1 538 / \ /|\ 539 / \ / | \ 540 / \ / | \ 541 / _*p\ -1/ | \ -1 542 v2--------v0 ========> v2 | v0 543 \ / \ | / 544 \ / \ | / 545 \ / \ | / 546 \ / \|/ +1 547 v3 v3 548 Before Swap After Swap 549 */ testSwap(PosType p,ScalarType creaseAngleCosThr)550 static bool testSwap(PosType p, ScalarType creaseAngleCosThr) 551 { 552 //if border or feature, do not swap 553 if (/*p.IsBorder() || */p.IsEdgeS()) return false; 554 555 int oldDist = 0, newDist = 0, idealV, actualV; 556 557 PosType tp=p; 558 559 VertexType *v0=tp.V(); 560 561 std::vector<VertexType*> incident; 562 563 vcg::face::VVStarVF<FaceType>(tp.V(), incident); 564 idealV = idealValence(tp); actualV = incident.size(); 565 oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV - 1)); 566 567 tp.NextF();tp.FlipE();tp.FlipV(); 568 VertexType *v1=tp.V(); 569 vcg::face::VVStarVF<FaceType>(tp.V(), incident); 570 idealV = idealValence(tp); actualV = incident.size(); 571 oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV + 1)); 572 573 tp.FlipE();tp.FlipV();tp.FlipE(); 574 VertexType *v2=tp.V(); 575 vcg::face::VVStarVF<FaceType>(tp.V(), incident); 576 idealV = idealValence(tp); actualV = incident.size(); 577 oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV - 1)); 578 579 tp.NextF();tp.FlipE();tp.FlipV(); 580 VertexType *v3=tp.V(); 581 vcg::face::VVStarVF<FaceType>(tp.V(), incident); 582 idealV = idealValence(tp); actualV = incident.size(); 583 oldDist += abs(idealV - actualV); newDist += abs(idealV - (actualV + 1)); 584 585 ScalarType qOld = std::min(Quality(v0->P(),v2->P(),v3->P()),Quality(v0->P(),v1->P(),v2->P())); 586 ScalarType qNew = std::min(Quality(v0->P(),v1->P(),v3->P()),Quality(v2->P(),v3->P(),v1->P())); 587 588 return (newDist < oldDist && qNew >= qOld * 0.50f) || 589 (newDist == oldDist && qNew > qOld * 1.f) || qNew > 1.5f * qOld; 590 } 591 checkManifoldness(FaceType & f,int z)592 static bool checkManifoldness(FaceType & f, int z) 593 { 594 PosType pos(&f, (z+2)%3, f.V2(z)); 595 PosType start = pos; 596 597 do { 598 pos.FlipE(); 599 if (!face::IsManifold(*pos.F(), pos.E())) 600 break; 601 pos.FlipF(); 602 } while (pos!=start); 603 604 return pos == start; 605 } 606 607 // Edge swap step: edges are flipped in order to optimize valence and triangle quality across the mesh ImproveValence(MeshType & m,Params & params)608 static void ImproveValence(MeshType &m, Params ¶ms) 609 { 610 static ScalarType foldCheckRad = math::ToRad(5.); 611 tri::UpdateTopology<MeshType>::FaceFace(m); 612 tri::UpdateTopology<MeshType>::VertexFace(m); 613 ForEachFace(m, [&] (FaceType & f) { 614 // if (face::IsManifold(f, 0) && face::IsManifold(f, 1) && face::IsManifold(f, 2)) 615 for (int i = 0; i < 3; ++i) 616 { 617 if (&f > f.cFFp(i)) 618 { 619 PosType pi(&f, i); 620 CoordType swapEdgeMidPoint = (f.cP2(i) + f.cFFp(i)->cP2(f.cFFi(i))) / 2.; 621 std::vector<CoordType> toCheck(1, swapEdgeMidPoint); 622 623 624 if(((!params.selectedOnly) || (f.IsS() && f.cFFp(i)->IsS())) && 625 !face::IsBorder(f, i) && 626 face::IsManifold(f, i) && /*checkManifoldness(f, i) &&*/ 627 face::checkFlipEdgeNotManifold(f, i) && 628 testSwap(pi, params.creaseAngleCosThr) && 629 (!params.surfDistCheck || testHausdorff(*params.mProject, params.grid, toCheck, params.maxSurfDist)) && 630 face::CheckFlipEdgeNormal(f, i, vcg::math::ToRad(5.))) 631 { 632 //When doing the swap we need to preserve and update the crease info accordingly 633 FaceType* g = f.cFFp(i); 634 int w = f.FFi(i); 635 636 bool creaseF = g->IsFaceEdgeS((w + 1) % 3); 637 bool creaseG = f.IsFaceEdgeS((i + 1) % 3); 638 639 face::FlipEdgeNotManifold(f, i); 640 641 f.ClearFaceEdgeS((i + 1) % 3); 642 g->ClearFaceEdgeS((w + 1) % 3); 643 644 if (creaseF) 645 f.SetFaceEdgeS(i); 646 if (creaseG) 647 g->SetFaceEdgeS(w); 648 649 ++params.stat.flipNum; 650 break; 651 } 652 } 653 } 654 }); 655 } 656 657 // The predicate that defines which edges should be split 658 class EdgeSplitAdaptPred 659 { 660 public: 661 int count = 0; 662 ScalarType length, lengthThr, minQ, maxQ; operator()663 bool operator()(PosType &ep) 664 { 665 ScalarType mult = math::ClampedLerp((ScalarType)0.5,(ScalarType)1.5, (((math::Abs(ep.V()->Q())+math::Abs(ep.VFlip()->Q()))/(ScalarType)2.0)/(maxQ-minQ))); 666 ScalarType dist = Distance(ep.V()->P(), ep.VFlip()->P()); 667 if(dist > std::max(mult*length,lengthThr*2)) 668 { 669 ++count; 670 return true; 671 } 672 else 673 return false; 674 } 675 }; 676 677 class EdgeSplitLenPred 678 { 679 public: 680 int count = 0; 681 ScalarType squaredlengthThr; operator()682 bool operator()(PosType &ep) 683 { 684 if(SquaredDistance(ep.V()->P(), ep.VFlip()->P()) > squaredlengthThr) 685 { 686 ++count; 687 return true; 688 } 689 else 690 return false; 691 } 692 }; 693 694 //Split pass: This pass uses the tri::RefineE from the vcglib to implement 695 //the refinement step, using EdgeSplitPred as a predicate to decide whether to split or not SplitLongEdges(MeshType & m,Params & params)696 static void SplitLongEdges(MeshType &m, Params ¶ms) 697 { 698 tri::UpdateTopology<MeshType>::FaceFace(m); 699 tri::MidPoint<MeshType> midFunctor(&m); 700 701 ScalarType minQ,maxQ; 702 if(params.adapt){ 703 computeVQualityDistrMinMax(m, minQ, maxQ); 704 EdgeSplitAdaptPred ep; 705 ep.minQ = minQ; 706 ep.maxQ = maxQ; 707 ep.length = params.maxLength; 708 ep.lengthThr = params.lengthThr; 709 tri::RefineE(m,midFunctor,ep); 710 params.stat.splitNum+=ep.count; 711 } 712 else { 713 EdgeSplitLenPred ep; 714 ep.squaredlengthThr = params.maxLength*params.maxLength; 715 tri::RefineMidpoint(m, ep, params.selectedOnly); 716 params.stat.splitNum+=ep.count; 717 } 718 } 719 VtoE(const int v0,const int v1)720 static int VtoE(const int v0, const int v1) 721 { 722 static /*constexpr*/ int Vmat[3][3] = { -1, 0, 2, 723 0, -1, 1, 724 2, 1, -1}; 725 return Vmat[v0][v1]; 726 } 727 728 checkCanMoveOnCollapse(PosType p,std::vector<FaceType * > & faces,std::vector<int> & vIdxes,Params & params)729 static bool checkCanMoveOnCollapse(PosType p, std::vector<FaceType*> & faces, std::vector<int> & vIdxes, Params ¶ms) 730 { 731 bool allIncidentFaceSelected = true; 732 733 PosType pi = p; 734 735 CoordType dEdgeVector = (p.V()->cP() - p.VFlip()->cP()).Normalize(); 736 737 int incidentFeatures = 0; 738 739 for (size_t i = 0; i < faces.size(); ++i) 740 // if (faces[i] != p.F() && faces[i] != p.FFlip()) 741 { 742 if (faces[i]->IsFaceEdgeS(VtoE(vIdxes[i], (vIdxes[i]+1)%3))) 743 { 744 incidentFeatures++; 745 CoordType movingEdgeVector0 = (faces[i]->cP1(vIdxes[i]) - faces[i]->cP(vIdxes[i])).Normalize(); 746 if (std::fabs(movingEdgeVector0 * dEdgeVector) < .9f || !p.IsEdgeS()) 747 return false; 748 } 749 if (faces[i]->IsFaceEdgeS(VtoE(vIdxes[i], (vIdxes[i]+2)%3))) 750 { 751 incidentFeatures++; 752 CoordType movingEdgeVector1 = (faces[i]->cP2(vIdxes[i]) - faces[i]->cP(vIdxes[i])).Normalize(); 753 if (std::fabs(movingEdgeVector1 * dEdgeVector) < .9f || !p.IsEdgeS()) 754 return false; 755 } 756 allIncidentFaceSelected &= faces[i]->IsS(); 757 } 758 759 if (incidentFeatures > 4) 760 return false; 761 762 return params.selectedOnly ? allIncidentFaceSelected : true; 763 } 764 checkFacesAfterCollapse(std::vector<FaceType * > & faces,PosType p,const Point3<ScalarType> & mp,Params & params,bool relaxed)765 static bool checkFacesAfterCollapse (std::vector<FaceType*> & faces, PosType p, const Point3<ScalarType> &mp, Params ¶ms, bool relaxed) 766 { 767 for (FaceType* f : faces) 768 { 769 if(!(*f).IsD() && f != p.F()) //i'm not a deleted face 770 { 771 PosType pi(f, p.V()); //same vertex 772 773 VertexType *v0 = pi.V(); 774 VertexType *v1 = pi.F()->V1(pi.VInd()); 775 VertexType *v2 = pi.F()->V2(pi.VInd()); 776 777 if( v1 == p.VFlip() || v2 == p.VFlip()) //i'm the other deleted face 778 continue; 779 780 //check on new face quality 781 { 782 ScalarType newQ = Quality(mp, v1->P(), v2->P()); 783 ScalarType oldQ = Quality(v0->P(), v1->P(), v2->P()); 784 785 if( newQ <= 0.5*oldQ ) 786 return false; 787 } 788 789 // we prevent collapse that makes edges too long (except for cross) 790 if(!relaxed) 791 if((Distance(mp, v1->P()) > params.maxLength || Distance(mp, v2->P()) > params.maxLength)) 792 return false; 793 794 Point3<ScalarType> oldN = NormalizedTriangleNormal(*(pi.F())); 795 Point3<ScalarType> newN = Normal(mp, v1->P(), v2->P()).Normalize(); 796 797 float div = fastAngle(oldN, newN); 798 if(div < .0f ) return false; 799 800 // // check on new face distance from original mesh 801 if (params.surfDistCheck) 802 { 803 std::vector<CoordType> points(4); 804 points[0] = (v1->cP() + v2->cP() + mp) / 3.; 805 points[1] = (v1->cP() + mp) / 2.; 806 points[2] = (v2->cP() + mp) / 2.; 807 points[3] = mp; 808 if (!testHausdorff(*(params.mProject), params.grid, points, params.maxSurfDist)) 809 return false; 810 } 811 } 812 } 813 return true; 814 } 815 816 817 //TODO: Refactor code and implement the correct set up of crease info when collapsing towards a crease edge checkCollapseFacesAroundVert1(PosType & p,Point3<ScalarType> & mp,Params & params,bool relaxed)818 static bool checkCollapseFacesAroundVert1(PosType &p, Point3<ScalarType> &mp, Params ¶ms, bool relaxed) 819 { 820 PosType p0 = p, p1 = p; 821 822 p1.FlipV(); 823 824 vector<int> vi0, vi1; 825 vector<FaceType*> ff0, ff1; 826 827 face::VFStarVF<FaceType>(p0.V(), ff0, vi0); 828 face::VFStarVF<FaceType>(p1.V(), ff1, vi1); 829 830 //check crease-moveability 831 bool moveable0 = checkCanMoveOnCollapse(p0, ff0, vi0, params) && !p0.V()->IsS(); 832 bool moveable1 = checkCanMoveOnCollapse(p1, ff1, vi1, params) && !p1.V()->IsS(); 833 834 //if both moveable => go to midpoint 835 // else collapse on movable one 836 if (!moveable0 && !moveable1) 837 return false; 838 839 //casting int(true) is always 1 and int(false) = =0 840 assert(int(true) == 1); 841 assert(int(false) == 0); 842 mp = (p0.V()->cP() * int(moveable1) + p1.V()->cP() * int(moveable0)) / (int(moveable0) + int(moveable1)); 843 844 if (!moveable0) 845 p = p0; 846 else 847 p = p1; 848 849 if (checkFacesAfterCollapse(ff0, p0, mp, params, relaxed)) 850 return checkFacesAfterCollapse(ff1, p1, mp, params, relaxed); 851 852 return false; 853 } 854 855 static bool testCollapse1(PosType &p, Point3<ScalarType> &mp, ScalarType minQ, ScalarType maxQ, Params ¶ms, bool relaxed = false) 856 { 857 ScalarType mult = (params.adapt) ? math::ClampedLerp((ScalarType)0.5,(ScalarType)1.5, (((math::Abs(p.V()->Q())+math::Abs(p.VFlip()->Q()))/(ScalarType)2.0)/(maxQ-minQ))) : (ScalarType)1; 858 ScalarType dist = Distance(p.V()->P(), p.VFlip()->P()); 859 ScalarType thr = mult*params.minLength; 860 ScalarType area = DoubleArea(*(p.F()))/2.f; 861 if(relaxed || (dist < thr || area < params.minLength*params.minLength/100.f))//if to collapse 862 { 863 return checkCollapseFacesAroundVert1(p, mp, params, relaxed); 864 } 865 return false; 866 } 867 868 //This function is especially useful to enforce feature preservation during collapses 869 //of boundary edges in planar or near planar section of the mesh chooseBoundaryCollapse(PosType & p,VertexPair & pair)870 static bool chooseBoundaryCollapse(PosType &p, VertexPair &pair) 871 { 872 Point3<ScalarType> collapseNV, collapsedNV0, collapsedNV1; 873 collapseNV = (p.V()->P() - p.VFlip()->P()).normalized(); 874 875 vector<VertexType*> vv; 876 face::VVStarVF<FaceType>(p.V(), vv); 877 878 for(VertexType *v: vv) 879 if(!(*v).IsD() && (*v).IsB() && v != p.VFlip()) //ignore non border 880 collapsedNV0 = ((*v).P() - p.VFlip()->P()).normalized(); //edge vector after collapse 881 882 face::VVStarVF<FaceType>(p.VFlip(), vv); 883 884 for(VertexType *v: vv) 885 if(!(*v).IsD() && (*v).IsB() && v != p.V()) //ignore non border 886 collapsedNV1 = ((*v).P() - p.V()->P()).normalized(); //edge vector after collapse 887 888 float cosine = cos(math::ToRad(1.5f)); 889 float angle0 = fabs(fastAngle(collapseNV, collapsedNV0)); 890 float angle1 = fabs(fastAngle(collapseNV, collapsedNV1)); 891 //if on both sides we deviate too much after collapse => don't collapse 892 if(angle0 <= cosine && angle1 <= cosine) 893 return false; 894 //choose the best collapse (the more parallel one to the previous edge..) 895 pair = (angle0 >= angle1) ? VertexPair(p.V(), p.VFlip()) : VertexPair(p.VFlip(), p.V()); 896 return true; 897 } 898 899 //The actual collapse step: foreach edge it is collapse iff TestCollapse returns true AND 900 // the linkConditions are preserved CollapseShortEdges(MeshType & m,Params & params)901 static void CollapseShortEdges(MeshType &m, Params ¶ms) 902 { 903 ScalarType minQ, maxQ; 904 int candidates = 0; 905 906 if(params.adapt) 907 computeVQualityDistrMinMax(m, minQ, maxQ); 908 909 tri::UpdateTopology<MeshType>::VertexFace(m); 910 tri::UpdateFlags<MeshType>::FaceBorderFromVF(m); 911 tri::UpdateFlags<MeshType>::VertexBorderFromFaceBorder(m); 912 913 SelectionStack<MeshType> ss(m); 914 ss.push(); 915 916 { 917 tri::UpdateTopology<MeshType>::FaceFace(m); 918 Clean<MeshType>::CountNonManifoldVertexFF(m,true); 919 920 //FROM NOW ON VSelection is NotManifold 921 922 for(auto fi=m.face.begin(); fi!=m.face.end(); ++fi) 923 if(!(*fi).IsD() && (params.selectedOnly == false || fi->IsS())) 924 { 925 for(auto i=0; i<3; ++i) 926 { 927 PosType pi(&*fi, i); 928 ++candidates; 929 VertexPair bp = VertexPair(pi.V(), pi.VFlip()); 930 Point3<ScalarType> mp = (pi.V()->P()+pi.VFlip()->P())/2.f; 931 932 if(testCollapse1(pi, mp, minQ, maxQ, params) && Collapser::LinkConditions(bp)) 933 { 934 //collapsing on pi.V() 935 bp = VertexPair(pi.VFlip(), pi.V()); 936 937 Collapser::Do(m, bp, mp, true); 938 ++params.stat.collapseNum; 939 break; 940 } 941 942 } 943 } 944 } 945 946 ss.pop(); 947 } 948 949 950 //Here I just need to check the faces of the cross, since the other faces are not 951 //affected by the collapse of the internal faces of the cross. testCrossCollapse(PosType & p,std::vector<FaceType * > ff,std::vector<int> vi,Point3<ScalarType> & mp,Params & params)952 static bool testCrossCollapse(PosType &p, std::vector<FaceType*> ff, std::vector<int> vi, Point3<ScalarType> &mp, Params ¶ms) 953 { 954 if(!checkFacesAfterCollapse(ff, p, mp, params, true)) 955 return false; 956 return true; 957 } 958 959 /* 960 *Choose the best way to collapse a cross based on the (external) cross vertices valence 961 *and resulting face quality 962 * +0 -1 963 * v1 v1 v1 964 * /| \ /|\ / \ 965 * / | \ / | \ / \ 966 * / | \ / | \ / \ 967 * / *p| \ -1/ | \ -1 +0/ \+0 968 * v0-------- v2 ========> v0 | v2 OR v0-------v2 969 * \ | / \ | / \ / 970 * \ | / \ | / \ / 971 * \ | / \ | / \ / 972 * \ | / \|/ +0 \ / -1 973 * v3 v3 v3 974 */ chooseBestCrossCollapse(PosType & p,VertexPair & bp,vector<FaceType * > & ff)975 static bool chooseBestCrossCollapse(PosType &p, VertexPair& bp, vector<FaceType*> &ff) 976 { 977 vector<VertexType*> vv0, vv1, vv2, vv3; 978 VertexType *v0, *v1, *v2, *v3; 979 980 v0 = p.F()->V1(p.VInd()); 981 v1 = p.F()->V2(p.VInd()); 982 983 984 bool crease[4] = {false, false, false, false}; 985 986 crease[0] = p.F()->IsFaceEdgeS(VtoE(p.VInd(), (p.VInd()+1)%3)); 987 crease[1] = p.F()->IsFaceEdgeS(VtoE(p.VInd(), (p.VInd()+2)%3)); 988 989 for(FaceType *f: ff) 990 if(!(*f).IsD() && f != p.F()) 991 { 992 PosType pi(f, p.V()); 993 VertexType *fv1 = pi.F()->V1(pi.VInd()); 994 VertexType *fv2 = pi.F()->V2(pi.VInd()); 995 996 if(fv1 == v0 || fv2 == v0) 997 { 998 if (fv1 == 0) 999 { 1000 v3 = fv2; 1001 crease[3] = f->IsFaceEdgeS(VtoE(pi.VInd(), (pi.VInd()+2)%3)); 1002 } 1003 else 1004 { 1005 v3 = fv1; 1006 crease[3] = f->IsFaceEdgeS(VtoE(pi.VInd(), (pi.VInd()+1)%3)); 1007 } 1008 // v3 = (fv1 == v0) ? fv2 : fv1; 1009 } 1010 1011 if(fv1 == v1 || fv2 == v1) 1012 { 1013 if (fv1 == v1) 1014 { 1015 v2 = fv2; 1016 crease[2] = f->IsFaceEdgeS(VtoE(pi.VInd(), (pi.VInd()+2)%3)); 1017 } 1018 else 1019 { 1020 v2 = fv1; 1021 crease[2] = f->IsFaceEdgeS(VtoE(pi.VInd(), (pi.VInd()+1)%3)); 1022 } 1023 // v2 = (fv1 == v1) ? fv2 : fv1; 1024 } 1025 } 1026 1027 face::VVStarVF<FaceType>(v0, vv0); 1028 face::VVStarVF<FaceType>(v1, vv1); 1029 face::VVStarVF<FaceType>(v2, vv2); 1030 face::VVStarVF<FaceType>(v3, vv3); 1031 1032 int nv0 = vv0.size(), nv1 = vv1.size(); 1033 int nv2 = vv2.size(), nv3 = vv3.size(); 1034 1035 int delta1 = (idealValence(*v0) - nv0) + (idealValence(*v2) - nv2); 1036 int delta2 = (idealValence(*v1) - nv1) + (idealValence(*v3) - nv3); 1037 1038 ScalarType Q1 = std::min(Quality(v0->P(), v1->P(), v3->P()), Quality(v1->P(), v2->P(), v3->P())); 1039 ScalarType Q2 = std::min(Quality(v0->P(), v1->P(), v2->P()), Quality(v2->P(), v3->P(), v0->P())); 1040 1041 if (crease[0] || crease[1] || crease[2] || crease[3]) 1042 return false; 1043 // if (crease[0] && crease[1] && crease[2] && crease[3]) 1044 // { 1045 // return false; 1046 // } 1047 1048 // if (crease[0] || crease[2]) 1049 // { 1050 // bp = VertexPair(p.V(), v0); 1051 // return true; 1052 // } 1053 1054 // if (crease[1] || crease[3]) 1055 // { 1056 // bp = VertexPair(p.V(), v1); 1057 // return true; 1058 // } 1059 1060 //no crease 1061 if(delta1 < delta2 && Q1 >= 0.6f*Q2) 1062 { 1063 bp = VertexPair(p.V(), v1); 1064 return true; 1065 } 1066 else 1067 { 1068 bp = VertexPair(p.V(), v0); 1069 return true; 1070 } 1071 } 1072 //Cross Collapse pass: This pass cleans the mesh from cross vertices, keeping in mind the link conditions 1073 //and feature preservations tests. CollapseCrosses(MeshType & m,Params & params)1074 static void CollapseCrosses(MeshType &m , Params ¶ms) 1075 { 1076 tri::UpdateTopology<MeshType>::VertexFace(m); 1077 tri::UpdateFlags<MeshType>::VertexBorderFromNone(m); 1078 int count = 0; 1079 1080 SelectionStack<MeshType> ss(m); 1081 ss.push(); 1082 1083 1084 { 1085 tri::UpdateTopology<MeshType>::FaceFace(m); 1086 Clean<MeshType>::CountNonManifoldVertexFF(m,true); 1087 1088 //From now on Selection on vertices is not manifoldness 1089 1090 for(auto fi=m.face.begin(); fi!=m.face.end(); ++fi) 1091 if(!(*fi).IsD() && (!params.selectedOnly || fi->IsS())) 1092 { 1093 for(auto i=0; i<3; ++i) 1094 { 1095 PosType pi(&*fi, i); 1096 if(!pi.V()->IsB()) 1097 { 1098 vector<FaceType*> ff; 1099 vector<int> vi; 1100 face::VFStarVF<FaceType>(pi.V(), ff, vi); 1101 1102 //if cross need to check what creases you have and decide where to collapse accordingly 1103 //if tricuspidis need whenever you have at least one crease => can't collapse anywhere 1104 if(ff.size() == 4 || ff.size() == 3) 1105 { 1106 // VertexPair bp; 1107 VertexPair bp = VertexPair(pi.V(), pi.VFlip()); 1108 Point3<ScalarType> mp = (pi.V()->P()+pi.VFlip()->P())/2.f; 1109 1110 if(testCollapse1(pi, mp, 0, 0, params, true) && Collapser::LinkConditions(bp)) 1111 { 1112 bp = VertexPair(pi.VFlip(), pi.V()); 1113 Collapser::Do(m, bp, mp, true); 1114 ++params.stat.collapseNum; 1115 ++count; 1116 break; 1117 } 1118 } 1119 } 1120 } 1121 } 1122 } 1123 1124 ss.pop(); 1125 Allocator<MeshType>::CompactEveryVector(m); 1126 } 1127 1128 // This function sets the selection bit on vertices that lie on creases selectVertexFromCrease(MeshType & m,ScalarType creaseThr)1129 static int selectVertexFromCrease(MeshType &m, ScalarType creaseThr) 1130 { 1131 int count = 0; 1132 Clean<MeshType>::CountNonManifoldVertexFF(m, true, false); 1133 1134 ForEachFacePos(m, [&](PosType &p){ 1135 if(p.IsBorder() || p.IsEdgeS()/*testCreaseEdge(p, creaseThr)*/) 1136 { 1137 p.V()->SetS(); 1138 p.VFlip()->SetS(); 1139 ++count; 1140 } 1141 }); 1142 return count; 1143 } 1144 selectVertexFromFold(MeshType & m,Params & params)1145 static int selectVertexFromFold(MeshType &m, Params & params) 1146 { 1147 std::vector<char> creaseVerts(m.VN(), 0); 1148 ForEachFacePos(m, [&] (PosType & p) { 1149 if (p.IsEdgeS()) 1150 { 1151 creaseVerts[vcg::tri::Index(m, p.V())] = 1; 1152 creaseVerts[vcg::tri::Index(m, p.VFlip())] = 1; 1153 } 1154 }); 1155 1156 1157 //this aspect ratio check doesn't work on cadish meshes (long thin triangles spanning whole mesh) 1158 ForEachFace(m, [&] (FaceType & f) { 1159 if (vcg::Quality(f.cP(0), f.cP(1), f.cP(2)) < params.aspectRatioThr || vcg::DoubleArea(f) < 0.00001) 1160 { 1161 if (creaseVerts[vcg::tri::Index(m, f.V(0))] == 0) 1162 f.V(0)->SetS(); 1163 if (creaseVerts[vcg::tri::Index(m, f.V(1))] == 0) 1164 f.V(1)->SetS(); 1165 if (creaseVerts[vcg::tri::Index(m, f.V(2))] == 0) 1166 f.V(2)->SetS(); 1167 } 1168 }); 1169 1170 1171 ForEachFace(m, [&] (FaceType & f) { 1172 for (int i = 0; i < 3; ++i) 1173 { 1174 if (f.FFp(i) > &f) 1175 { 1176 ScalarType angle = fastAngle(NormalizedTriangleNormal(f), NormalizedTriangleNormal(*(f.FFp(i)))); 1177 if (angle <= params.foldAngleCosThr) 1178 { 1179 if (creaseVerts[vcg::tri::Index(m, f.V0(i))] == 0) 1180 f.V0(i)->SetS(); 1181 if (creaseVerts[vcg::tri::Index(m, f.V1(i))] == 0) 1182 f.V1(i)->SetS(); 1183 if (creaseVerts[vcg::tri::Index(m, f.V2(i))] == 0) 1184 f.V2(i)->SetS(); 1185 if (creaseVerts[vcg::tri::Index(m, f.FFp(i)->V2(f.FFi(i)))] == 0) 1186 f.FFp(i)->V2(f.FFi(i))->SetS(); 1187 } 1188 } 1189 } 1190 }); 1191 1192 return 0; 1193 } 1194 1195 1196 1197 1198 static void FoldRelax(MeshType &m, Params params, const int step, const bool strict = true) 1199 { 1200 typename vcg::tri::Smooth<MeshType>::LaplacianInfo lpz(CoordType(0, 0, 0), 0); 1201 SimpleTempData<typename MeshType::VertContainer, typename vcg::tri::Smooth<MeshType>::LaplacianInfo> TD(m.vert, lpz); 1202 const ScalarType maxDist = (strict) ? params.maxSurfDist / 1000. : params.maxSurfDist; 1203 for (int i = 0; i < step; ++i) 1204 { 1205 TD.Init(lpz); 1206 vcg::tri::Smooth<MeshType>::AccumulateLaplacianInfo(m, TD, false); 1207 1208 for (auto fi = m.face.begin(); fi != m.face.end(); ++fi) 1209 { 1210 std::vector<CoordType> newPos(4); 1211 bool moving = false; 1212 1213 for (int j = 0; j < 3; ++j) 1214 { 1215 newPos[j] = fi->cP(j); 1216 if (!fi->V(j)->IsD() && TD[fi->V(j)].cnt > 0) 1217 { 1218 if (fi->V(j)->IsS()) 1219 { 1220 newPos[j] = (fi->V(j)->P() + TD[fi->V(j)].sum) / (TD[fi->V(j)].cnt + 1); 1221 moving = true; 1222 } 1223 } 1224 } 1225 1226 if (moving) 1227 { 1228 // const CoordType oldN = vcg::NormalizedTriangleNormal(*fi); 1229 // const CoordType newN = vcg::Normal(newPos[0], newPos[1], newPos[2]).Normalize(); 1230 1231 newPos[3] = (newPos[0] + newPos[1] + newPos[2]) / 3.; 1232 if (/*(strict || oldN * newN > 0.99) &&*/ (!params.surfDistCheck || testHausdorff(*params.mProject, params.grid, newPos, maxDist))) 1233 { 1234 for (int j = 0; j < 3; ++j) 1235 fi->V(j)->P() = newPos[j]; 1236 } 1237 } 1238 } 1239 } 1240 } 1241 1242 static void VertexCoordPlanarLaplacian(MeshType &m, Params & params, int step, ScalarType delta = 0.2) 1243 { 1244 typename vcg::tri::Smooth<MeshType>::LaplacianInfo lpz(CoordType(0, 0, 0), 0); 1245 SimpleTempData<typename MeshType::VertContainer, typename vcg::tri::Smooth<MeshType>::LaplacianInfo> TD(m.vert, lpz); 1246 for (int i = 0; i < step; ++i) 1247 { 1248 TD.Init(lpz); 1249 vcg::tri::Smooth<MeshType>::AccumulateLaplacianInfo(m, TD, false); 1250 // First normalize the AccumulateLaplacianInfo 1251 for (auto vi = m.vert.begin(); vi != m.vert.end(); ++vi) 1252 if (!(*vi).IsD() && TD[*vi].cnt > 0) 1253 { 1254 if ((*vi).IsS()) 1255 TD[*vi].sum = ((*vi).P() + TD[*vi].sum) / (TD[*vi].cnt + 1); 1256 } 1257 1258 for (auto fi = m.face.begin(); fi != m.face.end(); ++fi) 1259 { 1260 if (!(*fi).IsD()) 1261 { 1262 for (int j = 0; j < 3; ++j) 1263 { 1264 if (Angle(Normal(TD[(*fi).V0(j)].sum, (*fi).P1(j), (*fi).P2(j)), 1265 Normal((*fi).P0(j), (*fi).P1(j), (*fi).P2(j))) > M_PI/2.) 1266 TD[(*fi).V0(j)].sum = (*fi).P0(j); 1267 } 1268 } 1269 } 1270 for (auto fi = m.face.begin(); fi != m.face.end(); ++fi) 1271 { 1272 if (!(*fi).IsD()) 1273 { 1274 for (int j = 0; j < 3; ++j) 1275 { 1276 if (Angle(Normal(TD[(*fi).V0(j)].sum, TD[(*fi).V1(j)].sum, (*fi).P2(j)), 1277 Normal((*fi).P0(j), (*fi).P1(j), (*fi).P2(j))) > M_PI/2.) 1278 { 1279 TD[(*fi).V0(j)].sum = (*fi).P0(j); 1280 TD[(*fi).V1(j)].sum = (*fi).P1(j); 1281 } 1282 } 1283 } 1284 } 1285 1286 for (auto vi = m.vert.begin(); vi != m.vert.end(); ++vi) 1287 if (!(*vi).IsD() && TD[*vi].cnt > 0) 1288 { 1289 std::vector<CoordType> newPos(1, TD[*vi].sum); 1290 if ((*vi).IsS() && testHausdorff(*params.mProject, params.grid, newPos, params.maxSurfDist)) 1291 (*vi).P() = (*vi).P() * (1-delta) + TD[*vi].sum * (delta); 1292 } 1293 } // end step 1294 } 1295 1296 // static int 1297 /** 1298 * Simple Laplacian Smoothing step 1299 * Border and crease vertices are kept fixed. 1300 * If there are selected faces and the param.onlySelected is true we compute 1301 * the set of internal vertices to the selection and we combine it in and with 1302 * the vertexes not on border or creases 1303 */ ImproveByLaplacian(MeshType & m,Params params)1304 static void ImproveByLaplacian(MeshType &m, Params params) 1305 { 1306 SelectionStack<MeshType> ss(m); 1307 1308 if(params.selectedOnly) { 1309 ss.push(); 1310 tri::UpdateSelection<MeshType>::VertexFromFaceStrict(m); 1311 ss.push(); 1312 } 1313 tri::UpdateTopology<MeshType>::FaceFace(m); 1314 tri::UpdateFlags<MeshType>::VertexBorderFromFaceAdj(m); 1315 tri::UpdateSelection<MeshType>::VertexFromBorderFlag(m); 1316 selectVertexFromCrease(m, params.creaseAngleCosThr); 1317 tri::UpdateSelection<MeshType>::VertexInvert(m); 1318 if(params.selectedOnly) { 1319 ss.popAnd(); 1320 } 1321 1322 VertexCoordPlanarLaplacian(m, params, 1); 1323 1324 tri::UpdateSelection<MeshType>::VertexClear(m); 1325 1326 selectVertexFromFold(m, params); 1327 FoldRelax(m, params, 2); 1328 1329 tri::UpdateSelection<MeshType>::VertexClear(m); 1330 1331 if(params.selectedOnly) { 1332 ss.pop(); 1333 } 1334 } 1335 /* 1336 Reprojection step, this method reprojects each vertex on the original surface 1337 sampling the nearest Point3 onto it using a uniform grid StaticGrid t 1338 */ 1339 //TODO: improve crease reprojection: 1340 // crease verts should reproject only on creases. ProjectToSurface(MeshType & m,Params & params)1341 static void ProjectToSurface(MeshType &m, Params & params) 1342 { 1343 for(auto vi=m.vert.begin();vi!=m.vert.end();++vi) 1344 if(!(*vi).IsD()) 1345 { 1346 Point3<ScalarType> newP, normP, barP; 1347 ScalarType maxDist = params.maxSurfDist * 1.5f, minDist = 0.f; 1348 FaceType* fp = GetClosestFaceBase(*params.mProject, params.grid, vi->cP(), maxDist, minDist, newP, normP, barP); 1349 1350 if (fp != NULL) 1351 { 1352 vi->P() = newP; 1353 } 1354 } 1355 } 1356 }; 1357 } // end namespace tri 1358 } // end namespace vcg 1359 #endif 1360