1 /**************************************************************************** 2 * VCGLib o o * 3 * Visual and Computer Graphics Library o o * 4 * _ O _ * 5 * Copyright(C) 2004-2016 \/)\/ * 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 __VCGLIB_POLY_MESH_ALGORITHM 24 #define __VCGLIB_POLY_MESH_ALGORITHM 25 26 #include <vcg/complex/algorithms/update/normal.h> 27 #include <vcg/complex/complex.h> 28 #include <vcg/space/polygon3.h> 29 #include <vcg/complex/algorithms/update/color.h> 30 #include <vcg/complex/algorithms/closest.h> 31 #include <vcg/complex/algorithms/point_sampling.h> 32 #include <vcg/complex/algorithms/update/quality.h> 33 #include <wrap/io_trimesh/export_obj.h> 34 35 //define a temporary triangle mesh type 36 class TempFace; 37 class TempVertex; 38 39 struct TempUsedTypes: public vcg::UsedTypes<vcg::Use<TempVertex>::AsVertexType, 40 vcg::Use<TempFace>::AsFaceType>{}; 41 42 class TempVertex:public vcg::Vertex<TempUsedTypes, 43 vcg::vertex::Coord3d, 44 vcg::vertex::Normal3d, 45 vcg::vertex::BitFlags> 46 {}; 47 48 class TempFace:public vcg::Face<TempUsedTypes, 49 vcg::face::VertexRef, 50 vcg::face::BitFlags, 51 vcg::face::FFAdj, 52 vcg::face::Mark, 53 vcg::face::Normal3d> 54 {}; 55 56 57 class TempMesh: public vcg::tri::TriMesh< std::vector<TempVertex>,std::vector<TempFace > > 58 {}; 59 60 namespace vcg{ 61 62 /*! 63 \ingroup PolyMeshType 64 65 \headerfile color.h vcg/complex/algorithms/polygonal_algorithms.h 66 67 \brief processing and optimization of generic polygonal meshes. 68 69 This class is used to performs varisous kind of geometric optimization on generic polygonal mesh such as flattengin or imptove the shape of polygons. 70 */ 71 72 template <class PolyMeshType> 73 class PolygonalAlgorithm 74 { 75 typedef typename PolyMeshType::FaceType FaceType; 76 typedef typename PolyMeshType::VertexType VertexType; 77 typedef typename PolyMeshType::VertexPointer VertexPointer; 78 typedef typename PolyMeshType::CoordType CoordType; 79 typedef typename PolyMeshType::ScalarType ScalarType; 80 typedef typename vcg::face::Pos<FaceType> PosType; 81 SetFacePos(PolyMeshType & poly_m,int IndexF,std::vector<CoordType> & Pos)82 static void SetFacePos(PolyMeshType &poly_m, 83 int IndexF,std::vector<CoordType> &Pos) 84 { 85 poly_m.face[IndexF].Dealloc(); 86 poly_m.face[IndexF].Alloc(Pos.size()); 87 //std::cout<<Pos.size()<<std::endl; 88 int sizeV=poly_m.vert.size(); 89 for (size_t i=0;i<Pos.size();i++) 90 vcg::tri::Allocator<PolyMeshType>::AddVertex(poly_m,Pos[i]); 91 92 for (size_t i=0;i<Pos.size();i++) 93 poly_m.face[IndexF].V(i)=&poly_m.vert[sizeV+i]; 94 } 95 96 public: 97 SubdivideStep(PolyMeshType & poly_m)98 static void SubdivideStep(PolyMeshType &poly_m) 99 { 100 //get the barycenters 101 std::vector<CoordType> Bary; 102 for (size_t i=0;i<poly_m.face.size();i++) 103 { 104 CoordType bary(0,0,0); 105 for (size_t j=0;j<poly_m.face[i].VN();j++) 106 bary+=poly_m.face[i].P(j); 107 108 bary/=poly_m.face[i].VN(); 109 Bary.push_back(bary); 110 } 111 112 //get center of edge 113 std::map<std::pair<CoordType,CoordType>, CoordType> EdgeVert; 114 for (size_t i=0;i<poly_m.face.size();i++) 115 for (size_t j=0;j<poly_m.face[i].VN();j++) 116 { 117 CoordType Pos0=poly_m.face[i].P0(j); 118 CoordType Pos1=poly_m.face[i].P1(j); 119 CoordType Avg=(Pos0+Pos1)/2; 120 std::pair<CoordType,CoordType> Key(std::min(Pos0,Pos1),std::max(Pos0,Pos1)); 121 EdgeVert[Key]=Avg; 122 } 123 124 int sizeF=poly_m.face.size(); 125 for (size_t i=0;i<sizeF;i++) 126 { 127 //retrieve the sequence of pos 128 std::vector<CoordType> Pos; 129 for (size_t j=0;j<poly_m.face[i].VN();j++) 130 { 131 CoordType Pos0=poly_m.face[i].P0(j); 132 CoordType Pos1=poly_m.face[i].P1(j); 133 std::pair<CoordType,CoordType> Key0(std::min(Pos0,Pos1),std::max(Pos0,Pos1)); 134 Pos0=EdgeVert[Key0]; 135 Pos.push_back(Pos0); 136 Pos.push_back(Pos1); 137 } 138 //get also the barycenter 139 CoordType BaryP=Bary[i]; 140 141 //then retrieve the face 142 std::vector<CoordType> PosQ; 143 PosQ.push_back(Pos[0]); 144 PosQ.push_back(Pos[1]); 145 PosQ.push_back(Pos[2]); 146 PosQ.push_back(BaryP); 147 SetFacePos(poly_m,i,PosQ); 148 149 int sizeV=Pos.size(); 150 //int start=0; 151 for (size_t j=2;j<sizeV;j+=2) 152 { 153 vcg::tri::Allocator<PolyMeshType>::AddFaces(poly_m,1); 154 std::vector<CoordType> PosQ; 155 PosQ.push_back(Pos[(j)%Pos.size()]); 156 PosQ.push_back(Pos[(j+1)%Pos.size()]); 157 PosQ.push_back(Pos[(j+2)%Pos.size()]); 158 PosQ.push_back(BaryP); 159 //start+=2; 160 SetFacePos(poly_m,poly_m.face.size()-1,PosQ); 161 //break; 162 } 163 } 164 vcg::tri::Clean<PolyMeshType>::RemoveDuplicateVertex(poly_m); 165 vcg::tri::Allocator<PolyMeshType>::CompactEveryVector(poly_m); 166 } 167 CollapseEdges(PolyMeshType & poly_m,const std::vector<PosType> & CollapsePos,const std::vector<CoordType> & InterpPos)168 static bool CollapseEdges(PolyMeshType &poly_m, 169 const std::vector<PosType> &CollapsePos, 170 const std::vector<CoordType> &InterpPos) 171 { 172 173 //this set how to remap the vertices after deletion 174 std::map<VertexType*,VertexType*> VertexRemap; 175 176 vcg::tri::UpdateFlags<PolyMeshType>::VertexClearS(poly_m); 177 178 bool collapsed=false; 179 //go over all faces and check the ones needed to be deleted 180 for (size_t i=0;i<CollapsePos.size();i++) 181 { 182 FaceType *currF=CollapsePos[i].F(); 183 int IndexE=CollapsePos[i].E(); 184 size_t NumV=currF->VN(); 185 VertexType *v0=currF->V(IndexE); 186 VertexType *v1=currF->V((IndexE+1)%NumV); 187 188 //safety check 189 assert(v0!=v1); 190 191 if (v0->IsS())continue; 192 if (v1->IsS())continue; 193 194 //put on the same position 195 v0->P()=InterpPos[i]; 196 v1->P()=InterpPos[i]; 197 198 //select the the two vertices 199 v0->SetS(); 200 v1->SetS(); 201 202 //set the remap 203 VertexRemap[v1]=v0; 204 205 collapsed=true; 206 } 207 208 //then remap vertices 209 for (size_t i=0;i<poly_m.face.size();i++) 210 { 211 int NumV=poly_m.face[i].VN(); 212 for (int j=0;j<NumV;j++) 213 { 214 //get the two vertices of the edge 215 VertexType *v0=poly_m.face[i].V(j); 216 //see if it must substituted or not 217 if (VertexRemap.count(v0)==0)continue; 218 //in that case remap to the new one 219 VertexType *newV=VertexRemap[v0]; 220 //assign new vertex 221 poly_m.face[i].V(j)=newV; 222 } 223 } 224 225 //then re-elaborate the face 226 for (size_t i=0;i<poly_m.face.size();i++) 227 { 228 //get vertices of the face 229 int NumV=poly_m.face[i].VN(); 230 std::vector<VertexType*> FaceV; 231 for (int j=0;j<NumV;j++) 232 { 233 VertexType *v0=poly_m.face[i].V(j); 234 VertexType *v1=poly_m.face[i].V((j+1)%NumV); 235 if(v0==v1)continue; 236 FaceV.push_back(v0); 237 } 238 239 //then deallocate face 240 if ((int)FaceV.size()==NumV)continue; 241 242 //otherwise deallocate and set new vertices 243 poly_m.face[i].Dealloc(); 244 poly_m.face[i].Alloc(FaceV.size()); 245 for (size_t j=0;j<FaceV.size();j++) 246 poly_m.face[i].V(j)=FaceV[j]; 247 } 248 249 //remove unreferenced vertices 250 vcg::tri::Clean<PolyMeshType>::RemoveUnreferencedVertex(poly_m); 251 252 //and compact them 253 vcg::tri::Allocator<PolyMeshType>::CompactEveryVector(poly_m); 254 255 return collapsed; 256 } 257 private: 258 259 CollapseBorderSmallEdgesStep(PolyMeshType & poly_m,const ScalarType edge_limit)260 static bool CollapseBorderSmallEdgesStep(PolyMeshType &poly_m, 261 const ScalarType edge_limit) 262 { 263 //update topology 264 vcg::tri::UpdateTopology<PolyMeshType>::FaceFace(poly_m); 265 266 //update border vertices 267 vcg::tri::UpdateFlags<PolyMeshType>::VertexBorderFromFaceAdj(poly_m); 268 269 270 vcg::tri::UpdateSelection<PolyMeshType>::VertexCornerBorder(poly_m,math::ToRad(150.0)); 271 272 std::vector<PosType> CollapsePos; 273 std::vector<CoordType> InterpPos; 274 275 //go over all faces and check the ones needed to be deleted 276 for (size_t i=0;i<poly_m.face.size();i++) 277 { 278 int NumV=poly_m.face[i].VN(); 279 for (int j=0;j<NumV;j++) 280 { 281 VertexType *v0=poly_m.face[i].V(j); 282 VertexType *v1=poly_m.face[i].V((j+1)%NumV); 283 assert(v0!=v1); 284 285 bool IsBV0=v0->IsB(); 286 bool IsBV1=v1->IsB(); 287 288 bool IsS0=v0->IsS(); 289 bool IsS1=v1->IsS(); 290 291 if ((IsS0)&&(IsS1))continue; 292 293 //in these cases is not possible to collapse 294 if ((!IsBV0)&&(!IsBV1))continue; 295 bool IsBorderE=(poly_m.face[i].FFp(j)==&poly_m.face[i]); 296 if ((!IsBorderE)&&(IsBV0)&&(IsBV1))continue; 297 298 assert((IsBV0)||(IsBV1)); 299 CoordType pos0=v0->P(); 300 CoordType pos1=v1->P(); 301 ScalarType currL=(pos0-pos1).Norm(); 302 if (currL>edge_limit)continue; 303 304 //then collapse the point 305 CoordType CurrInterpPos; 306 if ((IsBV0)&&(!IsBV1))CurrInterpPos=pos0; 307 if ((!IsBV0)&&(IsBV1))CurrInterpPos=pos1; 308 if ((IsBV0)&&(IsBV1)) 309 { 310 if ((!IsS0)&&(!IsS1)) 311 CurrInterpPos=(pos0+pos1)/2.0; 312 else 313 { 314 if ((!IsS0)&&(IsS1)) 315 CurrInterpPos=pos1; 316 else 317 { 318 assert((IsS0)&&(!IsS1)); 319 CurrInterpPos=pos0; 320 } 321 } 322 } 323 CollapsePos.push_back(PosType(&poly_m.face[i],j)); 324 InterpPos.push_back(CurrInterpPos); 325 } 326 } 327 328 return CollapseEdges(poly_m,CollapsePos,InterpPos); 329 } 330 LaplacianPos(PolyMeshType & poly_m,std::vector<CoordType> & AvVert)331 static void LaplacianPos(PolyMeshType &poly_m,std::vector<CoordType> &AvVert) 332 { 333 //cumulate step 334 AvVert.clear(); 335 AvVert.resize(poly_m.vert.size(),CoordType(0,0,0)); 336 std::vector<ScalarType> AvSum(poly_m.vert.size(),0); 337 for (size_t i=0;i<poly_m.face.size();i++) 338 for (size_t j=0;j<(size_t)poly_m.face[i].VN();j++) 339 { 340 //get current vertex 341 VertexType *currV=poly_m.face[i].V(j); 342 //and its position 343 CoordType currP=currV->P(); 344 //cumulate over other positions 345 ScalarType W=vcg::PolyArea(poly_m.face[i]); 346 //assert(W!=0); 347 for (size_t k=0;k<(size_t)poly_m.face[i].VN();k++) 348 { 349 if (k==j) continue; 350 int IndexV=vcg::tri::Index(poly_m,poly_m.face[i].V(k)); 351 AvVert[IndexV]+=currP*W; 352 AvSum[IndexV]+=W; 353 } 354 } 355 356 //average step 357 for (size_t i=0;i<poly_m.vert.size();i++) 358 { 359 if (AvSum[i]==0)continue; 360 AvVert[i]/=AvSum[i]; 361 } 362 } 363 364 365 UpdateNormal(FaceType & F)366 static void UpdateNormal(FaceType &F) 367 { 368 F.N()=vcg::PolygonNormal(F); 369 } 370 UpdateNormalByFitting(FaceType & F)371 static void UpdateNormalByFitting(FaceType &F) 372 { 373 UpdateNormal(F); 374 vcg::Plane3<ScalarType> PlF; 375 PlF=PolyFittingPlane(F); 376 if ((PlF.Direction()*F.N())<0) 377 F.N()=-PlF.Direction(); 378 else 379 F.N()=PlF.Direction(); 380 } 381 382 public: 383 SelectIrregularInternal(PolyMeshType & poly_m)384 static void SelectIrregularInternal(PolyMeshType &poly_m) 385 { 386 vcg::tri::UpdateQuality<PolyMeshType>::VertexValence(poly_m); 387 vcg::tri::UpdateSelection<PolyMeshType>::VertexClear(poly_m); 388 for (size_t i=0;i<poly_m.vert.size();i++) 389 { 390 if (poly_m.vert[i].IsB())continue; 391 if (poly_m.vert[i].Q()==4)continue; 392 poly_m.vert[i].SetS(); 393 } 394 } 395 SelectIrregularBorder(PolyMeshType & poly_m)396 static void SelectIrregularBorder(PolyMeshType &poly_m) 397 { 398 vcg::tri::UpdateQuality<PolyMeshType>::VertexValence(poly_m); 399 for (size_t i=0;i<poly_m.vert.size();i++) 400 { 401 if (!poly_m.vert[i].IsB())continue; 402 if (poly_m.vert[i].Q()==2)continue; 403 poly_m.vert[i].SetS(); 404 } 405 } 406 GetFaceGetBary(FaceType & F)407 static CoordType GetFaceGetBary(FaceType &F) 408 { 409 CoordType bary=PolyBarycenter(F); 410 return bary; 411 } 412 413 /*! \brief update the face normal by averaging among vertex's 414 * normals computed between adjacent edges 415 */ UpdateFaceNormals(PolyMeshType & poly_m)416 static void UpdateFaceNormals(PolyMeshType &poly_m) 417 { 418 for (size_t i=0;i<poly_m.face.size();i++) 419 UpdateNormal(poly_m.face[i]); 420 } 421 422 /*! \brief update the face normal by fitting a plane 423 */ UpdateFaceNormalByFitting(PolyMeshType & poly_m)424 static void UpdateFaceNormalByFitting(PolyMeshType &poly_m) 425 { 426 for (size_t i=0;i<poly_m.face.size();i++) 427 UpdateNormalByFitting(poly_m.face[i]); 428 } 429 430 431 enum PolyQualityType{QAngle,QPlanar,QTemplate}; 432 433 /*! \brief update the quality of the faces by considering different possibilities 434 * QAngle = consider the angle deviation from ideal one (ex 90° quad, 60° triangle...) 435 * QPlanar = consider the difference wrt interpolating plane 436 * QTemplate= consider the difference wrt template polygon as in "Statics Aware Grid Shells" 437 */ UpdateQuality(PolyMeshType & poly_m,const PolyQualityType & QType)438 static void UpdateQuality(PolyMeshType &poly_m, 439 const PolyQualityType &QType) 440 { 441 for (size_t i=0;i<poly_m.face.size();i++) 442 { 443 if (poly_m.face[i].IsD())continue; 444 switch (QType) 445 { 446 case QAngle: 447 ScalarType AvgDev,WorstDev; 448 vcg::PolyAngleDeviation(poly_m.face[i],AvgDev,WorstDev); 449 poly_m.face[i].Q()=WorstDev; 450 break; 451 case QPlanar: 452 poly_m.face[i].Q()=vcg::PolyFlatness(poly_m.face[i]); 453 break; 454 default: 455 poly_m.face[i].Q()=vcg::PolyAspectRatio(poly_m.face[i],true); 456 break; 457 } 458 } 459 } 460 461 /*! \brief given a face this function returns the template positions as in "Statics Aware Grid Shells" 462 */ GetRotatedTemplatePos(FaceType & f,std::vector<CoordType> & TemplatePos)463 static void GetRotatedTemplatePos(FaceType &f, 464 std::vector<CoordType> &TemplatePos) 465 { 466 vcg::GetPolyTemplatePos(f,TemplatePos,true); 467 468 CoordType NormT=Normal(TemplatePos); 469 470 //get the normal of vertices 471 //CoordType AVN(0,0,0); 472 //CoordType AVN0(0,0,0); 473 CoordType Origin(0,0,0); 474 // for (int j=0;j<f.VN();j++) 475 // AVN0=AVN0+f.V(j)->N(); 476 477 CoordType AVN=vcg::PolygonNormal(f); 478 //AVN0.Normalize(); 479 // std::cout<<"AVN "<<AVN.X()<<","<<AVN.Y()<<","<<AVN.Z()<<std::endl; 480 // std::cout<<"AVN0 "<<AVN0.X()<<","<<AVN0.Y()<<","<<AVN0.Z()<<std::endl; 481 // std::cout<<"NormT "<<NormT.X()<<","<<NormT.Y()<<","<<NormT.Z()<<std::endl; 482 483 for (size_t j=0;j<TemplatePos.size();j++) 484 Origin+=TemplatePos[j]; 485 486 Origin/=(ScalarType)TemplatePos.size(); 487 AVN.Normalize(); 488 489 //find rotation matrix 490 vcg::Matrix33<ScalarType> Rot=vcg::RotationMatrix(NormT,AVN); 491 492 //apply transformation 493 for (size_t j=0;j<TemplatePos.size();j++) 494 { 495 TemplatePos[j]=TemplatePos[j]-Origin; 496 TemplatePos[j]=Rot*TemplatePos[j]; 497 TemplatePos[j]=TemplatePos[j]+Origin; 498 } 499 } 500 501 /*! \brief This function performs the polygon regularization as in "Statics Aware Grid Shells" 502 */ 503 static void SmoothPCA(PolyMeshType &poly_m, 504 int relax_step=10, 505 ScalarType Damp=0.5, 506 bool FixS=false, 507 bool isotropic=true, 508 ScalarType smoothTerm=0.1, 509 bool fixB=true, 510 bool WeightByQuality=false, 511 const std::vector<bool> *IgnoreF=NULL) 512 { 513 (void)isotropic; 514 typedef typename PolyMeshType::FaceType PolygonType; 515 // // select irregular ones 516 // if (fixIrr) 517 // poly_m.NumIrregular(true); 518 // compute the average edge 519 ScalarType MeshArea=0; 520 for (size_t i=0;i<poly_m.face.size();i++) 521 MeshArea+=vcg::PolyArea(poly_m.face[i]); 522 523 ScalarType AvgArea=MeshArea/(ScalarType)poly_m.face.size(); 524 525 if (WeightByQuality) 526 UpdateQuality(poly_m,QTemplate); 527 528 if (IgnoreF!=NULL){assert((*IgnoreF).size()==poly_m.face.size());} 529 for (size_t s=0;s<(size_t)relax_step;s++) 530 { 531 //initialize the accumulation vector 532 std::vector<CoordType> avgPos(poly_m.vert.size(),CoordType(0,0,0)); 533 std::vector<ScalarType> weightSum(poly_m.vert.size(),0); 534 //then compute the templated positions 535 536 for (size_t i=0;i<poly_m.face.size();i++) 537 { 538 if ((IgnoreF!=NULL)&&((*IgnoreF)[i]))continue; 539 std::vector<typename PolygonType::CoordType> TemplatePos; 540 GetRotatedTemplatePos(poly_m.face[i],TemplatePos); 541 //then cumulate the position per vertex 542 ScalarType val=vcg::PolyArea(poly_m.face[i]); 543 if (val<(AvgArea*0.00001)) 544 val=(AvgArea*0.00001); 545 546 ScalarType W=1.0/val; 547 if (WeightByQuality) 548 W=poly_m.face[i].Q()+0.00001; 549 550 for (size_t j=0;j<TemplatePos.size();j++) 551 { 552 int IndexV=vcg::tri::Index(poly_m,poly_m.face[i].V(j)); 553 CoordType Pos=TemplatePos[j]; 554 //sum up contributes 555 avgPos[IndexV]+=Pos*W; 556 weightSum[IndexV]+=W; 557 558 } 559 560 } 561 562 //get the laplacian contribute 563 std::vector<CoordType> AvVert; 564 LaplacianPos(poly_m,AvVert); 565 //then update the position 566 for (size_t i=0;i<poly_m.vert.size();i++) 567 { 568 ScalarType alpha=smoothTerm;//PolyNormDeviation(poly_m.face[i]); 569 // if (alpha<0)alpha=0; 570 // if (alpha>1)alpha=1; 571 // if (isnan(alpha))alpha=1; 572 573 CoordType newP=poly_m.vert[i].P(); 574 //safety checks 575 if (weightSum[i]>0) 576 newP=avgPos[i]/weightSum[i]; 577 if (isnan(newP.X())||isnan(newP.Y())||isnan(newP.Z())) 578 newP=poly_m.vert[i].P(); 579 if ((newP-poly_m.vert[i].P()).Norm()>poly_m.bbox.Diag()) 580 newP=poly_m.vert[i].P(); 581 //std::cout<<"W "<<weightSum[i]<<std::endl; 582 newP=newP*(1-alpha)+AvVert[i]*alpha; 583 //newP=AvVert[i]; 584 if ((fixB)&&(poly_m.vert[i].IsB()))continue; 585 if ((FixS)&&(poly_m.vert[i].IsS()))continue; 586 poly_m.vert[i].P()=poly_m.vert[i].P()*Damp+ 587 newP*(1-Damp); 588 } 589 } 590 } 591 592 template <class TriMeshType> 593 static void ReprojectBorder(PolyMeshType &poly_m, 594 TriMeshType &tri_mesh, 595 bool FixS=true) 596 { 597 //then reproject on border 598 for (size_t i=0;i<poly_m.vert.size();i++) 599 { 600 if (!poly_m.vert[i].IsB())continue; 601 if (FixS && poly_m.vert[i].IsS())continue; 602 603 CoordType testPos=poly_m.vert[i].P(); 604 ScalarType minD=std::numeric_limits<ScalarType>::max(); 605 CoordType closPos; 606 for (size_t j=0;j<tri_mesh.face.size();j++) 607 for (size_t k=0;k<3;k++) 608 { 609 //check if border edge 610 if (tri_mesh.face[j].FFp(k)!=(&tri_mesh.face[j]))continue; 611 612 CoordType P0,P1; 613 P0.Import(tri_mesh.face[j].cP0(k)); 614 P1.Import(tri_mesh.face[j].cP1(k)); 615 vcg::Segment3<ScalarType> Seg(P0,P1); 616 ScalarType testD; 617 CoordType closTest; 618 vcg::SegmentPointDistance(Seg,testPos,closTest,testD); 619 if (testD>minD)continue; 620 minD=testD; 621 closPos=closTest; 622 } 623 poly_m.vert[i].P()=closPos; 624 } 625 } 626 627 /*! \brief This function smooth the borders of the polygonal mesh and reproject back to the triangolar one 628 * except the vertices that are considered as corner wrt the angleDeg threshold 629 */ 630 template <class TriMeshType> 631 static void LaplacianReprojectBorder(PolyMeshType &poly_m, 632 TriMeshType &tri_mesh, 633 int nstep=100, 634 ScalarType Damp=0.5, 635 ScalarType angleDeg=100) 636 { 637 //first select corners 638 vcg::tri::UpdateFlags<PolyMeshType>::VertexClearS(poly_m); 639 640 //update topology 641 vcg::tri::UpdateTopology<PolyMeshType>::FaceFace(poly_m); 642 643 //update border vertices 644 vcg::tri::UpdateFlags<PolyMeshType>::VertexBorderFromFaceAdj(poly_m); 645 646 //select corner vertices on the border 647 ScalarType angleRad=angleDeg * M_PI / 180; 648 vcg::tri::UpdateSelection<PolyMeshType>::VertexCornerBorder(poly_m,angleRad); 649 650 for (int s=0;s<nstep;s++) 651 { 652 std::vector<CoordType> AvVert; 653 LaplacianPos(poly_m,AvVert); 654 655 for (size_t i=0;i<poly_m.vert.size();i++) 656 { 657 if (!poly_m.vert[i].IsB())continue; 658 if (poly_m.vert[i].IsS())continue; 659 poly_m.vert[i].P()=poly_m.vert[i].P()*Damp+ 660 AvVert[i]*(1-Damp); 661 } 662 663 // //then reproject on border 664 // for (size_t i=0;i<poly_m.vert.size();i++) 665 // { 666 // if (!poly_m.vert[i].IsB())continue; 667 // if (poly_m.vert[i].IsS())continue; 668 669 // CoordType testPos=poly_m.vert[i].P(); 670 // ScalarType minD=std::numeric_limits<ScalarType>::max(); 671 // CoordType closPos; 672 // for (size_t j=0;j<tri_mesh.face.size();j++) 673 // for (size_t k=0;k<3;k++) 674 // { 675 // if (tri_mesh.face[j].FFp(k)!=(&tri_mesh.face[j]))continue; 676 677 // CoordType P0,P1; 678 // P0.Import(tri_mesh.face[j].cP0(k)); 679 // P1.Import(tri_mesh.face[j].cP1(k)); 680 // vcg::Segment3<ScalarType> Seg(P0,P1); 681 // ScalarType testD; 682 // CoordType closTest; 683 // vcg::SegmentPointDistance(Seg,testPos,closTest,testD); 684 // if (testD>minD)continue; 685 // minD=testD; 686 // closPos=closTest; 687 // } 688 // poly_m.vert[i].P()=closPos; 689 // } 690 ReprojectBorder(poly_m,tri_mesh); 691 } 692 } 693 694 /*! \brief This function smooth the borders of the polygonal mesh and reproject back to its border 695 */ 696 static void LaplacianReprojectBorder(PolyMeshType &poly_m, 697 int nstep=100, 698 ScalarType Damp=0.5, 699 ScalarType Angle=100) 700 { 701 //transform into triangular 702 TempMesh GuideSurf; 703 vcg::tri::PolygonSupport<TempMesh,PolyMeshType>::ImportFromPolyMesh(GuideSurf,poly_m); 704 vcg::tri::UpdateBounding<TempMesh>::Box(GuideSurf); 705 vcg::tri::UpdateNormal<TempMesh>::PerVertexNormalizedPerFace(GuideSurf); 706 vcg::tri::UpdateTopology<TempMesh>::FaceFace(GuideSurf); 707 vcg::tri::UpdateFlags<TempMesh>::FaceBorderFromFF(GuideSurf); 708 709 LaplacianReprojectBorder<TempMesh>(poly_m,GuideSurf,nstep,Damp,Angle); 710 } 711 712 /*! \brief This function performs the reprojection of the polygonal mesh onto a triangular one passed as input parameter 713 */ 714 template <class TriMeshType> 715 static void LaplacianReproject(PolyMeshType &poly_m, 716 TriMeshType &tri_mesh, 717 int nstep=100, 718 ScalarType DampS=0.5, 719 ScalarType DampR=0.5, 720 bool OnlyOnSelected=false) 721 { 722 typedef typename TriMeshType::FaceType TriFaceType; 723 typedef typename TriMeshType::ScalarType TriScalarType; 724 typedef typename TriMeshType::CoordType TriCoordType; 725 typedef vcg::GridStaticPtr<TriFaceType, TriScalarType> TriMeshGrid; 726 TriMeshGrid grid; 727 tri::MeshAssert<TriMeshType>::VertexNormalNormalized(tri_mesh); 728 //initialize the grid 729 grid.Set(tri_mesh.face.begin(),tri_mesh.face.end()); 730 731 TriScalarType MaxD=tri_mesh.bbox.Diag(); 732 733 for (int s=0;s<nstep;s++) 734 { 735 std::vector<CoordType> AvVert; 736 LaplacianPos(poly_m,AvVert); 737 738 for (size_t i=0;i<poly_m.vert.size();i++) 739 { 740 if (poly_m.vert[i].IsB()) continue; 741 if (OnlyOnSelected && !poly_m.vert[i].IsS()) continue; 742 743 poly_m.vert[i].P()=poly_m.vert[i].P()*DampS+ 744 AvVert[i]*(1-DampS); 745 } 746 747 748 for (size_t i=0;i<poly_m.vert.size();i++) 749 { 750 if(OnlyOnSelected && !poly_m.vert[i].IsS()) continue; 751 TriCoordType testPos; 752 testPos.Import(poly_m.vert[i].P()); 753 TriCoordType closestPt; 754 TriScalarType minDist; 755 TriFaceType *f=NULL; 756 TriCoordType norm,ip; 757 f=vcg::tri::GetClosestFaceBase(tri_mesh,grid,testPos,MaxD,minDist,closestPt,norm,ip); 758 CoordType closestImp; 759 closestImp.Import(closestPt); 760 poly_m.vert[i].P()=poly_m.vert[i].P()*DampR+ 761 closestImp*(1-DampR); 762 CoordType normalImp; 763 normalImp.Import(norm); 764 poly_m.vert[i].N()=normalImp; 765 } 766 } 767 768 } 769 770 static void LaplacianReproject(PolyMeshType &poly_m, 771 int nstep=100, 772 ScalarType Damp=0.5, 773 bool OnlyOnSelected=false) 774 { 775 //transform into triangular 776 TempMesh GuideSurf; 777 //vcg::tri::PolygonSupport<TempMesh,PolyMeshType>:(GuideSurf,poly_m); 778 TriangulateToTriMesh<TempMesh>(poly_m,GuideSurf); 779 vcg::tri::UpdateBounding<TempMesh>::Box(GuideSurf); 780 vcg::tri::UpdateNormal<TempMesh>::PerVertexNormalizedPerFace(GuideSurf); 781 vcg::tri::UpdateTopology<TempMesh>::FaceFace(GuideSurf); 782 vcg::tri::UpdateFlags<TempMesh>::FaceBorderFromFF(GuideSurf); 783 LaplacianReproject<TempMesh>(poly_m,GuideSurf,nstep,Damp,0.5,OnlyOnSelected); 784 } 785 786 static void Laplacian(PolyMeshType &poly_m, 787 bool FixS=false, 788 int nstep=10, 789 ScalarType Damp=0.5) 790 { 791 for (int s=0;s<nstep;s++) 792 { 793 std::vector<CoordType> AvVert; 794 LaplacianPos(poly_m,AvVert); 795 796 for (size_t i=0;i<poly_m.vert.size();i++) 797 { 798 if ((FixS) && (poly_m.vert[i].IsS()))continue; 799 poly_m.vert[i].P()=poly_m.vert[i].P()*Damp+ 800 AvVert[i]*(1-Damp); 801 } 802 } 803 804 } 805 806 /*! \brief This function performs the polygon regularization as in "Statics Aware Grid Shells" 807 * followed by a reprojection step on the triangle mesh passed as parameter 808 */ 809 template <class TriMeshType> 810 static void SmoothReprojectPCA(PolyMeshType &poly_m, 811 TriMeshType &tri_mesh, 812 int relaxStep=100, 813 bool fixS=false, 814 ScalarType Damp=0.5, 815 ScalarType SharpDeg=0, 816 bool WeightByQuality=false, 817 bool FixB=true) 818 { 819 //vcg::tri::UpdateFlags<PolyMeshType>::VertexClearS(poly_m); 820 821 vcg::tri::UpdateTopology<PolyMeshType>::FaceFace(poly_m); 822 823 //UpdateBorderVertexFromPFFAdj(poly_m); 824 vcg::tri::UpdateFlags<PolyMeshType>::VertexBorderFromFaceAdj(poly_m); 825 826 std::vector<std::vector<vcg::Line3<ScalarType> > > SharpEdge(poly_m.vert.size()); 827 //first select sharp features 828 if (SharpDeg>0) 829 { 830 for (int i=0;i<(int)poly_m.face.size();i++) 831 for (int j=0;j<(int)poly_m.face[i].VN();j++) 832 { 833 //check only one side 834 if ((&poly_m.face[i])>=poly_m.face[i].FFp(j))continue; 835 836 CoordType N0=poly_m.face[i].N(); 837 CoordType N1=poly_m.face[i].FFp(j)->N(); 838 839 ScalarType Angle=vcg::Angle(N0,N1); 840 if (fabs(Angle)>(SharpDeg* (M_PI / 180.0))) 841 { 842 CoordType Pos0=poly_m.face[i].V0(j)->P(); 843 CoordType Pos1=poly_m.face[i].V1(j)->P(); 844 CoordType Ori=Pos0; 845 CoordType Dir=Pos1-Pos0; 846 Dir.Normalize(); 847 vcg::Line3<ScalarType> L(Ori,Dir); 848 int Index0=vcg::tri::Index(poly_m,poly_m.face[i].V0(j)); 849 int Index1=vcg::tri::Index(poly_m,poly_m.face[i].V1(j)); 850 SharpEdge[Index0].push_back(L); 851 SharpEdge[Index1].push_back(L); 852 } 853 } 854 for (size_t i=0;i<poly_m.vert.size();i++) 855 { 856 if (SharpEdge[i].size()==0)continue; 857 if (SharpEdge[i].size()>2)poly_m.vert[i].SetS(); 858 } 859 } 860 // if (fixIrr) 861 // { 862 // vcg::tri::UpdateQuality<PolyMeshType>::VertexValence(poly_m); 863 // for (size_t i=0;i<poly_m.vert.size();i++) 864 // { 865 // if (poly_m.vert[i].IsB())continue; 866 // if (poly_m.vert[i].Q()==4)continue; 867 // poly_m.vert[i].SetS(); 868 // } 869 // } 870 871 872 typedef typename TriMeshType::FaceType FaceType; 873 typedef vcg::GridStaticPtr<FaceType, typename TriMeshType::ScalarType> TriMeshGrid; 874 TriMeshGrid grid; 875 876 //initialize the grid 877 grid.Set(tri_mesh.face.begin(),tri_mesh.face.end()); 878 879 ScalarType MaxD=tri_mesh.bbox.Diag(); 880 881 // //update quality as area 882 // for (size_t i=0;i<poly_m.face.size();i++) 883 // poly_m.face[i].Q()=vcg::PolyArea(poly_m.face[i]); 884 885 // for (size_t i=0;i<poly_m.vert.size();i++) 886 // { 887 // typename TriMeshType::CoordType testPos; 888 // testPos.Import(poly_m.vert[i].P()); 889 // typename TriMeshType::CoordType closestPt; 890 // typename TriMeshType::ScalarType minDist; 891 // typename TriMeshType::FaceType *f=NULL; 892 // typename TriMeshType::CoordType norm,ip; 893 // f=vcg::tri::GetClosestFaceBase(tri_mesh,grid,testPos,MaxD,minDist,closestPt,norm,ip); 894 // //poly_m.vert[i].N().Import(norm); 895 // } 896 897 for(int k=0;k<relaxStep;k++) 898 { 899 //smooth PCA step 900 SmoothPCA(poly_m,1,Damp,fixS,true,0.1,FixB,WeightByQuality); 901 //reprojection step 902 //laplacian smooth step 903 //Laplacian(poly_m,Damp,1); 904 905 for (size_t i=0;i<poly_m.vert.size();i++) 906 { 907 typename TriMeshType::CoordType testPos; 908 testPos.Import(poly_m.vert[i].P()); 909 typename TriMeshType::CoordType closestPt; 910 typename TriMeshType::ScalarType minDist; 911 if ((FixB)&&(poly_m.vert[i].IsB())) 912 {continue;} 913 else 914 if (SharpEdge[i].size()==0)//reproject onto original mesh 915 { 916 FaceType *f=NULL; 917 typename TriMeshType::CoordType norm,ip; 918 f=vcg::tri::GetClosestFaceBase(tri_mesh,grid,testPos,MaxD,minDist,closestPt,norm,ip); 919 poly_m.vert[i].P().Import(testPos*Damp+closestPt*(1-Damp)); 920 //poly_m.vert[i].N().Import(norm); 921 } 922 else //reproject onto segments 923 { 924 CoordType av_closest(0,0,0); 925 size_t sum=0; 926 for (size_t j=0;j<SharpEdge[i].size();j++) 927 { 928 CoordType currPos; 929 currPos.Import(testPos); 930 CoordType closest; 931 ScalarType dist; 932 vcg::LinePointDistance(SharpEdge[i][j],currPos,closest,dist); 933 av_closest+=closest; 934 sum++; 935 } 936 assert(sum>0); 937 poly_m.vert[i].P()=av_closest/sum; 938 } 939 } 940 if (!FixB) 941 ReprojectBorder(poly_m,tri_mesh,true); 942 UpdateFaceNormals(poly_m); 943 vcg::tri::UpdateNormal<PolyMeshType>::PerVertexFromCurrentFaceNormal(poly_m); 944 } 945 946 } 947 948 949 template <class TriMeshType> 950 static void TriangulateToTriMesh(PolyMeshType &poly_m,TriMeshType &triangle_mesh, bool alsoTriangles = true) 951 { 952 triangle_mesh.Clear(); 953 954 PolyMeshType PolySwap; 955 vcg::tri::Append<PolyMeshType,PolyMeshType>::Mesh(PolySwap,poly_m); 956 Triangulate(PolySwap, alsoTriangles); 957 958 //then copy onto the triangle mesh 959 vcg::tri::Append<TriMeshType,PolyMeshType>::Mesh(triangle_mesh,PolySwap); 960 } 961 962 /*! \brief This function performs the polygon regularization as in "Statics Aware Grid Shells" 963 * followed by a reprojection step on the original mesh 964 */ 965 static void SmoothReprojectPCA(PolyMeshType &poly_m, 966 int relaxStep=100, 967 bool fixS=false, 968 ScalarType Damp=0.5, 969 ScalarType SharpDeg=0, 970 bool WeightByQuality=false, 971 bool FixB=true) 972 { 973 //transform into triangular 974 TempMesh GuideSurf; 975 976 //vcg::tri::PolygonSupport<TempMesh,PolyMeshType>:(GuideSurf,poly_m); 977 TriangulateToTriMesh<TempMesh>(poly_m,GuideSurf); 978 vcg::tri::UpdateBounding<TempMesh>::Box(GuideSurf); 979 vcg::tri::UpdateNormal<TempMesh>::PerVertexNormalizedPerFace(GuideSurf); 980 vcg::tri::UpdateTopology<TempMesh>::FaceFace(GuideSurf); 981 vcg::tri::UpdateFlags<TempMesh>::FaceBorderFromFF(GuideSurf); 982 983 //optimize it 984 vcg::PolygonalAlgorithm<PolyMeshType>::SmoothReprojectPCA<TempMesh>(poly_m,GuideSurf,relaxStep,fixS,Damp,SharpDeg,WeightByQuality,FixB); 985 } 986 Reproject(PolyMeshType & poly_m,PolyMeshType & target)987 static void Reproject(PolyMeshType &poly_m, 988 PolyMeshType &target) 989 { 990 vcg::tri::UpdateTopology<PolyMeshType>::FaceFace(poly_m); 991 vcg::tri::UpdateFlags<PolyMeshType>::VertexBorderFromFaceAdj(poly_m); 992 993 //transform into triangular 994 TempMesh GuideSurf; 995 //vcg::tri::PolygonSupport<TempMesh,PolyMeshType>:(GuideSurf,poly_m); 996 TriangulateToTriMesh<TempMesh>(target,GuideSurf); 997 vcg::tri::UpdateBounding<TempMesh>::Box(GuideSurf); 998 vcg::tri::UpdateNormal<TempMesh>::PerVertexNormalizedPerFace(GuideSurf); 999 vcg::tri::UpdateTopology<TempMesh>::FaceFace(GuideSurf); 1000 vcg::tri::UpdateFlags<TempMesh>::FaceBorderFromFF(GuideSurf); 1001 1002 //initialize the grid 1003 typedef typename TempMesh::FaceType FaceType; 1004 typedef vcg::GridStaticPtr<FaceType, typename TempMesh::ScalarType> TriMeshGrid; 1005 TriMeshGrid grid; 1006 grid.Set(GuideSurf.face.begin(),GuideSurf.face.end()); 1007 1008 ScalarType MaxD=GuideSurf.bbox.Diag(); 1009 1010 for (size_t i=0;i<poly_m.vert.size();i++) 1011 { 1012 //reproject on border later 1013 if (poly_m.vert[i].IsB())continue; 1014 typename TempMesh::CoordType testPos; 1015 testPos.Import(poly_m.vert[i].P()); 1016 typename TempMesh::CoordType closestPt; 1017 typename TempMesh::ScalarType minDist; 1018 typename TempMesh::FaceType *f=NULL; 1019 typename TempMesh::CoordType norm,ip; 1020 f=vcg::tri::GetClosestFaceBase(GuideSurf,grid,testPos,MaxD,minDist,closestPt,norm,ip); 1021 poly_m.vert[i].P()=closestPt; 1022 } 1023 //then reprojec the border 1024 ReprojectBorder(poly_m,GuideSurf); 1025 } 1026 1027 template <class TriMesh> ReprojectonTriMesh(PolyMeshType & poly_m,TriMesh & target)1028 static void ReprojectonTriMesh(PolyMeshType &poly_m, 1029 TriMesh &target) 1030 { 1031 vcg::tri::UpdateTopology<PolyMeshType>::FaceFace(poly_m); 1032 vcg::tri::UpdateFlags<PolyMeshType>::VertexBorderFromFaceAdj(poly_m); 1033 1034 1035 //initialize the grid 1036 typedef typename TriMesh::FaceType FaceType; 1037 typedef vcg::GridStaticPtr<FaceType, typename TriMesh::ScalarType> TriMeshGrid; 1038 TriMeshGrid grid; 1039 grid.Set(target.face.begin(),target.face.end()); 1040 1041 ScalarType MaxD=target.bbox.Diag(); 1042 1043 for (size_t i=0;i<poly_m.vert.size();i++) 1044 { 1045 //reproject on border later 1046 if (poly_m.vert[i].IsB())continue; 1047 typename TriMesh::CoordType testPos; 1048 testPos.Import(poly_m.vert[i].P()); 1049 typename TriMesh::CoordType closestPt; 1050 typename TriMesh::ScalarType minDist; 1051 typename TriMesh::FaceType *f=NULL; 1052 typename TriMesh::CoordType norm,ip; 1053 f=vcg::tri::GetClosestFaceBase(target,grid,testPos,MaxD,minDist,closestPt,norm,ip); 1054 poly_m.vert[i].P()=closestPt; 1055 } 1056 //then reprojec the border 1057 ReprojectBorder(poly_m,target); 1058 } 1059 1060 /*! \brief This function return average edge size 1061 */ AverageEdge(const PolyMeshType & poly_m)1062 static ScalarType AverageEdge(const PolyMeshType &poly_m) 1063 { 1064 ScalarType AvL=0; 1065 size_t numE=0; 1066 for (size_t i=0;i<poly_m.face.size();i++) 1067 { 1068 int NumV=poly_m.face[i].VN(); 1069 for (int j=0;j<NumV;j++) 1070 { 1071 CoordType pos0=poly_m.face[i].cV(j)->P(); 1072 CoordType pos1=poly_m.face[i].cV((j+1)%NumV)->P(); 1073 AvL+=(pos0-pos1).Norm(); 1074 numE++; 1075 } 1076 } 1077 AvL/=numE; 1078 return AvL; 1079 } 1080 1081 1082 /*! \brief This function remove valence 2 faces from the mesh 1083 */ RemoveValence2Faces(PolyMeshType & poly_m)1084 static void RemoveValence2Faces(PolyMeshType &poly_m) 1085 { 1086 for (size_t i=0;i<poly_m.face.size();i++) 1087 { 1088 if (poly_m.face[i].VN()>=3)continue; 1089 vcg::tri::Allocator<PolyMeshType>::DeleteFace(poly_m,poly_m.face[i]); 1090 } 1091 1092 //then remove unreferenced vertices 1093 vcg::tri::Clean<PolyMeshType>::RemoveUnreferencedVertex(poly_m); 1094 vcg::tri::Allocator<PolyMeshType>::CompactEveryVector(poly_m); 1095 1096 } 1097 1098 /*! \brief This function remove valence 2 vertices on the border by considering the degree threshold 1099 * bacause there could be eventually some corner that should be preserved 1100 */ 1101 static void RemoveValence2Vertices(PolyMeshType &poly_m, 1102 ScalarType corner_degree=25) 1103 { 1104 //update topology 1105 vcg::tri::UpdateTopology<PolyMeshType>::FaceFace(poly_m); 1106 1107 //update border vertices 1108 //UpdateBorderVertexFromPFFAdj(poly_m); 1109 vcg::tri::UpdateFlags<PolyMeshType>::VertexBorderFromFaceAdj(poly_m); 1110 1111 vcg::tri::UpdateFlags<PolyMeshType>::VertexClearS(poly_m); 1112 1113 //select corners 1114 for (size_t i=0;i<poly_m.face.size();i++) 1115 { 1116 if (poly_m.face[i].IsD())continue; 1117 1118 //get vertices of the face 1119 int NumV=poly_m.face[i].VN(); 1120 1121 for (int j=0;j<NumV;j++) 1122 { 1123 VertexType *v0=poly_m.face[i].V((j+NumV-1)%NumV); 1124 VertexType *v1=poly_m.face[i].V(j); 1125 VertexType *v2=poly_m.face[i].V((j+1)%NumV); 1126 //must be 3 borders 1127 bool IsB=((v0->IsB())&&(v1->IsB())&&(v2->IsB())); 1128 CoordType dir0=(v0->P()-v1->P()); 1129 CoordType dir1=(v2->P()-v1->P()); 1130 dir0.Normalize(); 1131 dir1.Normalize(); 1132 ScalarType testDot=(dir0*dir1); 1133 if ((IsB)&&(testDot>(-cos(corner_degree* (M_PI / 180.0))))) 1134 v1->SetS(); 1135 } 1136 } 1137 1138 typename PolyMeshType::template PerVertexAttributeHandle<size_t> valenceVertH = 1139 vcg::tri::Allocator<PolyMeshType>:: template GetPerVertexAttribute<size_t> (poly_m); 1140 1141 //initialize to zero 1142 for (size_t i=0;i<poly_m.vert.size();i++) 1143 valenceVertH[i]=0; 1144 1145 //then sum up the valence 1146 for (size_t i=0;i<poly_m.face.size();i++) 1147 for (int j=0;j<poly_m.face[i].VN();j++) 1148 valenceVertH[poly_m.face[i].V(j)]++; 1149 1150 //then re-elaborate the faces 1151 for (size_t i=0;i<poly_m.face.size();i++) 1152 { 1153 if (poly_m.face[i].IsD())continue; 1154 1155 //get vertices of the face 1156 int NumV=poly_m.face[i].VN(); 1157 1158 std::vector<VertexType*> FaceV; 1159 for (int j=0;j<NumV;j++) 1160 { 1161 VertexType *v=poly_m.face[i].V(j); 1162 assert(!v->IsD()); 1163 //if ((!v->IsS()) && (v->IsB()) && (valenceVertH[v]==1)) continue; 1164 if ((!v->IsS()) && (v->IsB()) && (valenceVertH[v]==1)) continue; 1165 if ((!v->IsB()) && (valenceVertH[v]<3)) continue; 1166 //if (!v->IsS()) continue; 1167 FaceV.push_back(v); 1168 } 1169 1170 //then deallocate face 1171 if ((int)FaceV.size()==NumV)continue; 1172 1173 //otherwise deallocate and set new vertices 1174 poly_m.face[i].Dealloc(); 1175 poly_m.face[i].Alloc(FaceV.size()); 1176 for (size_t j=0;j<FaceV.size();j++) 1177 poly_m.face[i].V(j)=FaceV[j]; 1178 } 1179 1180 //then remove unreferenced vertices 1181 vcg::tri::Clean<PolyMeshType>::RemoveUnreferencedVertex(poly_m); 1182 vcg::tri::Allocator<PolyMeshType>::CompactEveryVector(poly_m); 1183 1184 vcg::tri::Allocator<PolyMeshType>::DeletePerVertexAttribute(poly_m,valenceVertH); 1185 } 1186 1187 /*! \brief This function collapse small edges which are on the boundary of the mesh 1188 * this is sometimes useful to remove small edges coming out from a quadrangulation which is not 1189 * aligned to boundaries 1190 */ 1191 static bool CollapseBorderSmallEdges(PolyMeshType &poly_m, 1192 const ScalarType perc_average=0.3) 1193 { 1194 //compute the average edge 1195 ScalarType AvEdge=AverageEdge(poly_m); 1196 ScalarType minLimit=AvEdge*perc_average; 1197 bool collapsed=false; 1198 while(CollapseBorderSmallEdgesStep(poly_m,minLimit)){collapsed=true;}; 1199 1200 RemoveValence2Faces(poly_m); 1201 1202 //RemoveValence2BorderVertices(poly_m); 1203 RemoveValence2Vertices(poly_m); 1204 return collapsed; 1205 } 1206 1207 /*! \brief This function use a local global approach to flatten polygonal faces 1208 * the approach is similar to "Shape-Up: Shaping Discrete Geometry with Projections" 1209 */ 1210 static ScalarType FlattenFaces(PolyMeshType &poly_m, size_t steps=100,bool OnlySFaces=false) 1211 { 1212 ScalarType MaxDispl=0; 1213 for (size_t s=0;s<steps;s++) 1214 { 1215 std::vector<std::vector<CoordType> > VertPos(poly_m.vert.size()); 1216 1217 for (size_t i=0;i<poly_m.face.size();i++) 1218 { 1219 if (poly_m.face[i].IsD())continue; 1220 1221 if (OnlySFaces && (!poly_m.face[i].IsS()))continue; 1222 //get vertices of the face 1223 int NumV=poly_m.face[i].VN(); 1224 if (NumV<=3)continue; 1225 1226 //save vertice's positions 1227 std::vector<CoordType> FacePos; 1228 for (int j=0;j<NumV;j++) 1229 { 1230 VertexType *v=poly_m.face[i].V(j); 1231 assert(!v->IsD()); 1232 FacePos.push_back(v->P()); 1233 } 1234 1235 //then fit the plane 1236 vcg::Plane3<ScalarType> FitPl; 1237 vcg::FitPlaneToPointSet(FacePos,FitPl); 1238 1239 //project each point onto fitting plane 1240 for (int j=0;j<NumV;j++) 1241 { 1242 VertexType *v=poly_m.face[i].V(j); 1243 int IndexV=vcg::tri::Index(poly_m,v); 1244 CoordType ProjP=FitPl.Projection(v->P()); 1245 VertPos[IndexV].push_back(ProjP); 1246 } 1247 } 1248 1249 for (size_t i=0;i<poly_m.vert.size();i++) 1250 { 1251 CoordType AvgPos(0,0,0); 1252 1253 for (size_t j=0;j<VertPos[i].size();j++) 1254 AvgPos+=VertPos[i][j]; 1255 1256 if (VertPos[i].size()==0)continue; 1257 1258 AvgPos/=(ScalarType)VertPos[i].size(); 1259 1260 MaxDispl=std::max(MaxDispl,(poly_m.vert[i].P()-AvgPos).Norm()); 1261 poly_m.vert[i].P()=AvgPos; 1262 } 1263 } 1264 return MaxDispl; 1265 } 1266 Area(PolyMeshType & poly_m)1267 static ScalarType Area(PolyMeshType &poly_m) 1268 { 1269 ScalarType MeshArea=0; 1270 for (size_t i=0;i<poly_m.face.size();i++) 1271 MeshArea+=vcg::PolyArea(poly_m.face[i]); 1272 return MeshArea; 1273 } 1274 InitQualityVertVoronoiArea(PolyMeshType & poly_m)1275 static void InitQualityVertVoronoiArea(PolyMeshType &poly_m) 1276 { 1277 for (size_t i=0;i<poly_m.vert.size();i++) 1278 poly_m.vert[i].Q()=0; 1279 1280 for (size_t i=0;i<poly_m.face.size();i++) 1281 { 1282 // ScalarType AreaF=vcg::PolyArea(poly_m.face[i]); 1283 size_t sizeV=poly_m.face[i].VN()-1; 1284 CoordType baryF=vcg::PolyBarycenter(poly_m.face[i]); 1285 for (int j=0;j<poly_m.face[i].VN();j++) 1286 { 1287 CoordType P0=poly_m.face[i].P((j+sizeV-1)%sizeV); 1288 CoordType P1=poly_m.face[i].P(j); 1289 CoordType P2=poly_m.face[i].P1(j); 1290 vcg::Triangle3<ScalarType> T0(P1,(P0+P1)/2,baryF); 1291 vcg::Triangle3<ScalarType> T1(P1,(P1+P2)/2,baryF); 1292 1293 poly_m.face[i].V(j)->Q()+=vcg::DoubleArea(T0)/2; 1294 poly_m.face[i].V(j)->Q()+=vcg::DoubleArea(T1)/2; 1295 } 1296 } 1297 } 1298 InitQualityFaceTorsion(PolyMeshType & poly_m)1299 static ScalarType InitQualityFaceTorsion(PolyMeshType &poly_m) 1300 { 1301 UpdateFaceNormalByFitting(poly_m); 1302 vcg::tri::UpdateNormal<PolyMeshType>::PerVertexFromCurrentFaceNormal(poly_m); 1303 ScalarType MaxA=0; 1304 for (size_t i=0;i<poly_m.face.size();i++) 1305 { 1306 poly_m.face[i].Q()=PolygonTorsion(poly_m.face[i]); 1307 MaxA=std::max(MaxA,poly_m.face[i].Q()); 1308 } 1309 return MaxA; 1310 } 1311 InitQualityFaceBending(PolyMeshType & poly_m)1312 static ScalarType InitQualityFaceBending(PolyMeshType &poly_m) 1313 { 1314 UpdateFaceNormalByFitting(poly_m); 1315 vcg::tri::UpdateNormal<PolyMeshType>::PerVertexFromCurrentFaceNormal(poly_m); 1316 ScalarType MaxA=0; 1317 for (size_t i=0;i<poly_m.face.size();i++) 1318 { 1319 poly_m.face[i].Q()=PolygonBending(poly_m.face[i]); 1320 MaxA=std::max(MaxA,poly_m.face[i].Q()); 1321 } 1322 return MaxA; 1323 } 1324 1325 InitQualityVertEdgeLenght(PolyMeshType & poly_m)1326 static void InitQualityVertEdgeLenght(PolyMeshType &poly_m) 1327 { 1328 for (size_t i=0;i<poly_m.vert.size();i++) 1329 poly_m.vert[i].Q()=0; 1330 1331 for (size_t i=0;i<poly_m.face.size();i++) 1332 { 1333 for (int j=0;j<poly_m.face[i].VN();j++) 1334 { 1335 FaceType *f=&poly_m.face[i]; 1336 FaceType *f1=f->FFp(j); 1337 if (f>f1)continue; 1338 ScalarType L=(poly_m.face[i].P0(j)-poly_m.face[i].P1(j)).Norm(); 1339 poly_m.face[i].V0(j)->Q()+=L; 1340 poly_m.face[i].V1(j)->Q()+=L; 1341 } 1342 } 1343 } 1344 InterpolateQualityVertFormFaces(PolyMeshType & poly_m)1345 static void InterpolateQualityVertFormFaces(PolyMeshType &poly_m) 1346 { 1347 std::vector<ScalarType> SumW(poly_m.vert.size(),0); 1348 1349 for (size_t i=0;i<poly_m.vert.size();i++) 1350 poly_m.vert[i].Q()=0; 1351 1352 for (size_t i=0;i<poly_m.face.size();i++) 1353 { 1354 ScalarType AreaF=vcg::PolyArea(poly_m.face[i]); 1355 for (size_t j=0;j<poly_m.face[i].VN();j++) 1356 { 1357 poly_m.face[i].V(j)->Q()+=AreaF*(ScalarType)poly_m.face[i].Q(); 1358 size_t IndexV=vcg::tri::Index(poly_m,poly_m.face[i].V(j)); 1359 SumW[IndexV]+=AreaF; 1360 } 1361 } 1362 for (size_t i=0;i<poly_m.vert.size();i++) 1363 { 1364 if (SumW[i]>0) 1365 poly_m.vert[i].Q()/=SumW[i]; 1366 else 1367 poly_m.vert[i].Q()=0; 1368 } 1369 } 1370 1371 ClosestPoint(const PolyMeshType & poly_m,const CoordType & pos,int & CloseF,CoordType & ClosePos)1372 static void ClosestPoint(const PolyMeshType &poly_m,const CoordType &pos, 1373 int &CloseF,CoordType &ClosePos) 1374 { 1375 ScalarType minD=std::numeric_limits<ScalarType>::max(); 1376 CloseF=-1; 1377 for (size_t i=0;i<poly_m.face.size();i++) 1378 { 1379 CoordType closeTest; 1380 ScalarType currD=vcg::PolygonPointDistance(poly_m.face[i],pos,closeTest); 1381 if (currD>minD)continue; 1382 minD=currD; 1383 CloseF=i; 1384 ClosePos=closeTest; 1385 } 1386 } 1387 1388 /*! \brief Triangulate a polygonal face with a triangle fan. 1389 * \returns pointer to the newly added vertex. 1390 */ Triangulate(PolyMeshType & poly_m,size_t IndexF)1391 static VertexPointer Triangulate(PolyMeshType & poly_m, size_t IndexF) 1392 { 1393 1394 const CoordType bary = vcg::PolyBarycenter(poly_m.face[IndexF]); 1395 size_t sizeV = poly_m.face[IndexF].VN(); 1396 1397 //add the new vertex 1398 VertexPointer newV = &(*vcg::tri::Allocator<PolyMeshType>::AddVertex(poly_m,bary)); 1399 1400 //then reupdate the faces 1401 for (size_t j=0;j<(sizeV-1);j++) 1402 { 1403 VertexType * v0=poly_m.face[IndexF].V0(j); 1404 VertexType * v1=poly_m.face[IndexF].V1(j); 1405 VertexType * v2=newV; 1406 1407 vcg::tri::Allocator<PolyMeshType>::AddFaces(poly_m,1); 1408 1409 poly_m.face.back().Alloc(3); 1410 poly_m.face.back().V(0)=v0; 1411 poly_m.face.back().V(1)=v1; 1412 poly_m.face.back().V(2)=v2; 1413 poly_m.face.back().Q()=poly_m.face[IndexF].Q(); 1414 } 1415 1416 VertexType * v0=poly_m.face[IndexF].V0((sizeV-1)); 1417 VertexType * v1=poly_m.face[IndexF].V1((sizeV-1)); 1418 poly_m.face[IndexF].Dealloc(); 1419 poly_m.face[IndexF].Alloc(3); 1420 poly_m.face[IndexF].V(0)=v0; 1421 poly_m.face[IndexF].V(1)=v1; 1422 poly_m.face[IndexF].V(2)=newV; 1423 return newV; 1424 } 1425 ReorderFaceVert(FaceType & f,const size_t & StartI)1426 static void ReorderFaceVert(FaceType &f,const size_t &StartI) 1427 { 1428 if (StartI==0)return; 1429 size_t sizeN=f.VN(); 1430 assert(StartI>=0); 1431 assert(StartI<sizeN); 1432 std::vector<VertexType*> NewV; 1433 for (size_t i=0;i<sizeN;i++) 1434 { 1435 int IndexV=(i+StartI)%sizeN; 1436 NewV.push_back(f.V(IndexV)); 1437 } 1438 //then reset all vertices 1439 for (size_t i=0;i<sizeN;i++) 1440 f.V(i)=NewV[i]; 1441 } 1442 MergeAlongEdge(PolyMeshType & poly_m,FaceType & f,const size_t & EdgeI)1443 static void MergeAlongEdge(PolyMeshType &poly_m, 1444 FaceType &f, 1445 const size_t &EdgeI) 1446 { 1447 //cannot be a border 1448 assert(f.FFp(EdgeI)!=&f); 1449 FaceType *f1=f.FFp(EdgeI); 1450 int EdgeI1=f.FFi(EdgeI); 1451 1452 //sort first face 1453 int FirstV0=(EdgeI+1) % f.VN(); 1454 ReorderFaceVert(f,FirstV0); 1455 int FirstV1=(EdgeI1+1)%f1->VN(); 1456 ReorderFaceVert(*f1,FirstV1); 1457 1458 std::vector<VertexType*> NewV; 1459 for (size_t i=0;i<(f.VN()-1);i++) 1460 NewV.push_back(f.V(i)); 1461 1462 for (size_t i=0;i<(f1->VN()-1);i++) 1463 NewV.push_back(f1->V(i)); 1464 1465 f.Dealloc(); 1466 f.Alloc(NewV.size()); 1467 for (size_t i=0;i<NewV.size();i++) 1468 f.V(i)=NewV[i]; 1469 1470 vcg::tri::Allocator<PolyMeshType>::DeleteFace(poly_m,*f1); 1471 } 1472 MergeAlongEdges(PolyMeshType & poly_m,const std::vector<FaceType * > & PolyF,const std::vector<size_t> & EdgeI)1473 static void MergeAlongEdges(PolyMeshType &poly_m, 1474 const std::vector<FaceType*> &PolyF, 1475 const std::vector<size_t> &EdgeI) 1476 { 1477 //create a table with all edges that have to be merged 1478 std::set<std::pair<CoordType,CoordType> > NeedMerge; 1479 for (size_t i=0;i<PolyF.size();i++) 1480 { 1481 CoordType P0=PolyF[i]->P0(EdgeI[i]); 1482 CoordType P1=PolyF[i]->P1(EdgeI[i]); 1483 std::pair<CoordType,CoordType> key(std::min(P0,P1),std::max(P0,P1)); 1484 NeedMerge.insert(key); 1485 } 1486 1487 //then cycle and collapse 1488 do{ 1489 for (size_t i=0;i<poly_m.face.size();i++) 1490 { 1491 if (poly_m.face[i].IsD())continue; 1492 for (size_t j=0;j<poly_m.face[i].VN();j++) 1493 { 1494 CoordType P0=poly_m.face[i].P0(j); 1495 CoordType P1=poly_m.face[i].P1(j); 1496 std::pair<CoordType,CoordType> key(std::min(P0,P1),std::max(P0,P1)); 1497 if (NeedMerge.count(key)==0)continue; 1498 1499 //do the merge 1500 MergeAlongEdge(poly_m,poly_m.face[i],j); 1501 //remove it 1502 NeedMerge.erase(key); 1503 break; 1504 } 1505 } 1506 vcg::tri::Allocator<PolyMeshType>::CompactEveryVector(poly_m); 1507 }while (!NeedMerge.empty()); 1508 } 1509 1510 static void Triangulate(PolyMeshType &poly_m, 1511 bool alsoTriangles = true, 1512 bool OnlyS=false) 1513 { 1514 size_t size0 = poly_m.face.size(); 1515 if (alsoTriangles) 1516 { 1517 for (size_t i=0; i<size0; i++) 1518 { 1519 if ((OnlyS)&&(!poly_m.face[i].IsS()))continue; 1520 Triangulate(poly_m, i); 1521 } 1522 } 1523 else 1524 { 1525 for (size_t i=0; i<size0; i++) 1526 { 1527 if ((OnlyS)&&(!poly_m.face[i].IsS()))continue; 1528 if (poly_m.face[i].VN() > 3) 1529 { 1530 Triangulate(poly_m, i); 1531 } 1532 } 1533 } 1534 } 1535 1536 }; 1537 1538 }//end namespace vcg 1539 1540 #endif 1541