1 /**************************************************************************** 2 * VCGLib o o * 3 * Visual and Computer Graphics Library o o * 4 * _ O _ * 5 * Copyright(C) 2014 \/)\/ * 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 24 #ifndef SMOOTHER_FIELD_H 25 #define SMOOTHER_FIELD_H 26 27 //eigen stuff 28 #include <eigenlib/Eigen/Sparse> 29 30 //vcg stuff 31 #include <vcg/complex/algorithms/update/color.h> 32 #include <vcg/complex/algorithms/update/quality.h> 33 #include <vcg/complex/algorithms/parametrization/tangent_field_operators.h> 34 #include <vcg/complex/algorithms/mesh_to_matrix.h> 35 36 //igl related stuff 37 #include <igl/n_polyvector.h> 38 #include <igl/principal_curvature.h> 39 #include <igl/igl_inline.h> 40 #include <igl/comiso/nrosy.h> 41 42 namespace vcg { 43 namespace tri { 44 45 enum SmoothMethod{SMMiq,SMNPoly}; 46 47 template <class MeshType> 48 class FieldSmoother 49 { 50 typedef typename MeshType::FaceType FaceType; 51 typedef typename MeshType::VertexType VertexType; 52 typedef typename MeshType::ScalarType ScalarType; 53 typedef typename MeshType::CoordType CoordType; 54 55 InitQualityByAnisotropyDir(MeshType & mesh)56 static void InitQualityByAnisotropyDir(MeshType &mesh) 57 { 58 std::vector<ScalarType> QVal; 59 for (size_t i=0;i<mesh.vert.size();i++) 60 { 61 ScalarType N1=fabs(mesh.vert[i].K1()); 62 ScalarType N2=fabs(mesh.vert[i].K2()); 63 QVal.push_back(N1); 64 QVal.push_back(N2); 65 } 66 67 std::sort(QVal.begin(),QVal.end()); 68 int percUp=int(floor(QVal.size()*0.95+0.5)); 69 ScalarType trimUp=QVal[percUp]; 70 71 for (size_t i=0;i<mesh.vert.size();i++) 72 { 73 ScalarType N1=(mesh.vert[i].K1()); 74 ScalarType N2=(mesh.vert[i].K2()); 75 76 ScalarType NMax=std::max(N1,N2)/trimUp; 77 ScalarType NMin=std::min(N1,N2)/trimUp; 78 79 if (NMax>1)NMax=1; 80 if (NMax<-1)NMax=-1; 81 82 if (NMin>1)NMin=1; 83 if (NMin<-1)NMin=-1; 84 85 ScalarType CurvAni=(NMax-NMin)/2; 86 mesh.vert[i].Q()=CurvAni; 87 } 88 vcg::tri::UpdateQuality<MeshType>::FaceFromVertex(mesh); 89 } 90 SetEdgeDirection(FaceType * f,int edge)91 static void SetEdgeDirection(FaceType *f,int edge) 92 { 93 CoordType dir=f->P0(edge)-f->P1(edge); 94 dir.Normalize(); 95 ScalarType prod1=fabs(dir*f->PD1()); 96 ScalarType prod2=fabs(dir*f->PD2()); 97 if (prod1>prod2) 98 { 99 f->PD1()=dir; 100 f->PD2()=f->N()^dir; 101 }else 102 { 103 f->PD2()=dir; 104 f->PD1()=f->N()^dir; 105 } 106 } 107 108 static void AddSharpEdgesConstraints(MeshType & mesh, 109 const ScalarType &thr=0.2) 110 { 111 for (size_t i=0;i<mesh.face.size();i++) 112 for (int j=0;j<mesh.face[i].VN();j++) 113 { 114 FaceType *f0=&mesh.face[i]; 115 FaceType *f1=f0->FFp(j); 116 if (f0==f1)continue; 117 CoordType N0=f0->N(); 118 CoordType N1=f1->N(); 119 if ((N0*N1)>thr)continue; 120 SetEdgeDirection(f0,j); 121 f0->SetS(); 122 } 123 } 124 AddBorderConstraints(MeshType & mesh)125 static void AddBorderConstraints(MeshType & mesh) 126 { 127 for (size_t i=0;i<mesh.face.size();i++) 128 for (int j=0;j<mesh.face[i].VN();j++) 129 { 130 FaceType *f0=&mesh.face[i]; 131 FaceType *f1=f0->FFp(j); 132 assert(f1!=NULL); 133 if (f0!=f1)continue; 134 SetEdgeDirection(f0,j); 135 f0->SetS(); 136 } 137 } 138 139 static void AddCurvatureConstraints(MeshType & mesh,const ScalarType &thr=0.5) 140 { 141 for (size_t i=0;i<mesh.face.size();i++) 142 if (mesh.face[i].Q()>thr)mesh.face[i].SetS(); 143 } 144 145 //hard constraints have selected face CollectHardConstraints(MeshType & mesh,Eigen::VectorXi & HardI,Eigen::MatrixXd & HardD,SmoothMethod SMethod,int Ndir)146 static void CollectHardConstraints( MeshType & mesh, 147 Eigen::VectorXi &HardI, 148 Eigen::MatrixXd &HardD, 149 SmoothMethod SMethod, 150 int Ndir) 151 { 152 //count number of hard constraints 153 int numS=vcg::tri::UpdateSelection<MeshType>::FaceCount(mesh); 154 HardI=Eigen::MatrixXi(numS,1); 155 if ((Ndir==2)||(SMethod==SMMiq)) 156 HardD=Eigen::MatrixXd(numS,3); 157 else 158 HardD=Eigen::MatrixXd(numS,6); 159 //then update them 160 int curr_index=0; 161 for (size_t i=0;i<mesh.face.size();i++) 162 { 163 if (!mesh.face[i].IsS())continue; 164 165 CoordType dir=mesh.face[i].PD1(); 166 dir.Normalize(); 167 168 HardI(curr_index,0)=i; 169 170 HardD(curr_index,0)=dir.X(); 171 HardD(curr_index,1)=dir.Y(); 172 HardD(curr_index,2)=dir.Z(); 173 if ((Ndir==4)&&(SMethod==SMNPoly)) 174 { 175 dir=mesh.face[i].PD2(); 176 HardD(curr_index,3)=dir.X(); 177 HardD(curr_index,4)=dir.Y(); 178 HardD(curr_index,5)=dir.Z(); 179 } 180 curr_index++; 181 182 } 183 184 } 185 186 //hard constraints have selected face CollectSoftConstraints(MeshType & mesh,Eigen::VectorXi & SoftI,Eigen::MatrixXd & SoftD,Eigen::VectorXd & SoftW)187 static void CollectSoftConstraints( MeshType & mesh, 188 Eigen::VectorXi &SoftI, 189 Eigen::MatrixXd &SoftD, 190 Eigen::VectorXd &SoftW) 191 { 192 //count number of soft constraints 193 int numS=vcg::tri::UpdateSelection<MeshType>::FaceCount(mesh); 194 numS=mesh.fn-numS; 195 //allocate eigen matrix 196 SoftI=Eigen::MatrixXi(numS,1); 197 SoftD=Eigen::MatrixXd(numS,3); 198 SoftW=Eigen::MatrixXd(numS,1); 199 200 //then update them 201 int curr_index=0; 202 for (size_t i=0;i<mesh.face.size();i++) 203 { 204 if (mesh.face[i].IsS())continue; 205 206 CoordType dir=mesh.face[i].PD1(); 207 dir.Normalize(); 208 209 SoftI(curr_index,0)=i; 210 211 SoftD(curr_index,0)=dir.X(); 212 SoftD(curr_index,1)=dir.Y(); 213 SoftD(curr_index,2)=dir.Z(); 214 215 SoftW(curr_index,0)=mesh.face[i].Q(); 216 217 curr_index++; 218 219 } 220 } 221 SmoothMIQ(MeshType & mesh,Eigen::VectorXi & HardI,Eigen::MatrixXd & HardD,Eigen::VectorXi & SoftI,Eigen::MatrixXd & SoftD,Eigen::VectorXd & SoftW,ScalarType alpha_soft,int Ndir)222 static void SmoothMIQ(MeshType &mesh, 223 Eigen::VectorXi &HardI, //hard constraints index 224 Eigen::MatrixXd &HardD, //hard directions 225 Eigen::VectorXi &SoftI, //soft constraints 226 Eigen::MatrixXd &SoftD, //weight of soft constraints 227 Eigen::VectorXd &SoftW, //soft directions 228 ScalarType alpha_soft, 229 int Ndir) 230 { 231 232 assert((Ndir==2)||(Ndir==4)); 233 Eigen::MatrixXi F; 234 Eigen::MatrixXd V; 235 236 MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V); 237 238 Eigen::MatrixXd output_field; 239 Eigen::VectorXd output_sing; 240 241 igl::nrosy(V,F,HardI,HardD,SoftI,SoftW,SoftD,Ndir,alpha_soft,output_field,output_sing); 242 243 //finally update the principal directions 244 for (size_t i=0;i<mesh.face.size();i++) 245 { 246 CoordType dir1; 247 dir1[0]=output_field(i,0); 248 dir1[1]=output_field(i,1); 249 dir1[2]=output_field(i,2); 250 251 dir1.Normalize(); 252 CoordType dir2=mesh.face[i].N()^dir1; 253 dir2.Normalize(); 254 255 ScalarType Norm1=mesh.face[i].PD1().Norm(); 256 ScalarType Norm2=mesh.face[i].PD2().Norm(); 257 258 mesh.face[i].PD1()=dir1*Norm1; 259 mesh.face[i].PD2()=dir2*Norm2; 260 } 261 } 262 SmoothNPoly(MeshType & mesh,Eigen::VectorXi & HardI,Eigen::MatrixXd & HardD,int Ndir)263 static void SmoothNPoly(MeshType &mesh, 264 Eigen::VectorXi &HardI, //hard constraints index 265 Eigen::MatrixXd &HardD, //hard directions 266 int Ndir) 267 { 268 assert((Ndir==2)||(Ndir==4)); 269 270 Eigen::MatrixXi F; 271 Eigen::MatrixXd V; 272 273 MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V); 274 275 Eigen::MatrixXd output_field; 276 Eigen::VectorXd output_sing; 277 278 igl::n_polyvector(V,F,HardI,HardD,output_field); 279 280 //finally update the principal directions 281 for (size_t i=0;i<mesh.face.size();i++) 282 { 283 CoordType dir1; 284 dir1[0]=output_field(i,0); 285 dir1[1]=output_field(i,1); 286 dir1[2]=output_field(i,2); 287 288 dir1.Normalize(); 289 CoordType dir2=mesh.face[i].N()^dir1; 290 dir2.Normalize(); 291 292 ScalarType Norm1=mesh.face[i].PD1().Norm(); 293 ScalarType Norm2=mesh.face[i].PD2().Norm(); 294 295 mesh.face[i].PD1()=dir1*Norm1; 296 mesh.face[i].PD2()=dir2*Norm2; 297 } 298 } 299 PickRandomDir(MeshType & mesh,int & indexF,CoordType & Dir)300 static void PickRandomDir(MeshType &mesh, 301 int &indexF, 302 CoordType &Dir) 303 { 304 indexF=rand()%mesh.fn; 305 FaceType *currF=&mesh.face[indexF]; 306 CoordType dirN=currF->N(); 307 dirN.Normalize(); 308 Dir=CoordType(1,0,0); 309 if (fabs(Dir*dirN)>0.9) 310 Dir=CoordType(0,1,0); 311 if (fabs(Dir*dirN)>0.9) 312 Dir=CoordType(0,0,1); 313 314 Dir=dirN^Dir; 315 Dir.Normalize(); 316 } 317 318 public: 319 320 struct SmoothParam 321 { 322 //the 90° rotation independence while smoothing the direction field 323 int Ndir; 324 //the weight of curvature if doing the smoothing keeping the field close to the original one 325 ScalarType alpha_curv; 326 //align the field to border or not 327 bool align_borders; 328 //threshold to consider some edge as sharp feature and to use as hard constraint (0, not use) 329 ScalarType sharp_thr; 330 //threshold to consider some edge as high curvature anisotropyand to use as hard constraint (0, not use) 331 ScalarType curv_thr; 332 //the method used to smooth MIQ or "Designing N-PolyVector Fields with Complex Polynomials" 333 SmoothMethod SmoothM; 334 //the number of faces of the ring used ot esteem the curvature 335 int curvRing; 336 //this are additional hard constraints 337 std::vector<std::pair<int,CoordType> > AddConstr; 338 SmoothParamSmoothParam339 SmoothParam() 340 { 341 Ndir=4; 342 curvRing=2; 343 alpha_curv=0.0; 344 345 align_borders=false; 346 347 SmoothM=SMMiq; 348 349 sharp_thr=0.0; 350 curv_thr=0.4; 351 } 352 353 }; 354 SelectConstraints(MeshType & mesh,SmoothParam & SParam)355 static void SelectConstraints(MeshType &mesh,SmoothParam &SParam) 356 { 357 //clear all selected faces 358 vcg::tri::UpdateFlags<MeshType>::FaceClear(mesh); 359 360 //add curvature hard constraints 361 //ScalarType Ratio=mesh.bbox.Diag()*0.01; 362 363 if (SParam.curv_thr>0) 364 AddCurvatureConstraints(mesh,SParam.curv_thr);///Ratio); 365 366 //add alignment to sharp features 367 if (SParam.sharp_thr>0) 368 AddSharpEdgesConstraints(mesh,SParam.sharp_thr); 369 370 //add border constraints 371 if (SParam.align_borders) 372 AddBorderConstraints(mesh); 373 374 //aff final constraints 375 for (int i=0;i<SParam.AddConstr.size();i++) 376 { 377 int indexF=SParam.AddConstr[i].first; 378 CoordType dir=SParam.AddConstr[i].second; 379 mesh.face[indexF].PD1()=dir; 380 mesh.face[indexF].PD2()=mesh.face[indexF].N()^dir; 381 mesh.face[indexF].PD1().Normalize(); 382 mesh.face[indexF].PD2().Normalize(); 383 mesh.face[indexF].SetS(); 384 } 385 } 386 GloballyOrient(MeshType & mesh)387 static void GloballyOrient(MeshType &mesh) 388 { 389 vcg::tri::CrossField<MeshType>::MakeDirectionFaceCoherent(mesh,true); 390 } 391 392 393 static void InitByCurvature(MeshType & mesh, 394 int Nring, 395 bool UpdateFaces=true) 396 { 397 398 tri::RequirePerVertexCurvatureDir(mesh); 399 400 Eigen::MatrixXi F; 401 Eigen::MatrixXd V; 402 403 Eigen::MatrixXd PD1,PD2,PV1,PV2; 404 MeshToMatrix<MeshType>::GetTriMeshData(mesh,F,V); 405 igl::principal_curvature(V,F,PD1,PD2,PV1,PV2,Nring); 406 407 //then copy curvature per vertex 408 for (size_t i=0;i<mesh.vert.size();i++) 409 { 410 mesh.vert[i].PD1()=CoordType(PD1(i,0),PD1(i,1),PD1(i,2)); 411 mesh.vert[i].PD2()=CoordType(PD2(i,0),PD2(i,1),PD2(i,2)); 412 mesh.vert[i].K1()=PV1(i,0); 413 mesh.vert[i].K2()=PV2(i,0); 414 } 415 if (!UpdateFaces)return; 416 vcg::tri::CrossField<MeshType>::SetFaceCrossVectorFromVert(mesh); 417 InitQualityByAnisotropyDir(mesh); 418 } 419 420 static void SmoothDirections(MeshType &mesh, 421 int Ndir, 422 SmoothMethod SMethod=SMNPoly, 423 bool HardAsS=true, 424 ScalarType alphaSoft=0) 425 { 426 427 Eigen::VectorXi HardI; //hard constraints 428 Eigen::MatrixXd HardD; //hard directions 429 Eigen::VectorXi SoftI; //soft constraints 430 Eigen::VectorXd SoftW; //weight of soft constraints 431 Eigen::MatrixXd SoftD; //soft directions 432 433 if (HardAsS) 434 CollectHardConstraints(mesh,HardI,HardD,SMethod,Ndir); 435 436 //collect soft constraints , miw only one that allows for soft constraints 437 if ((alphaSoft>0)&&(SMethod==SMMiq)) 438 CollectSoftConstraints(mesh,SoftI,SoftD,SoftW); 439 440 //add some hard constraints if are not present 441 int numC=3; 442 if ((SoftI.rows()==0)&&(HardI.rows()==0)) 443 { 444 printf("Add Forced Constraints \n"); 445 fflush(stdout); 446 HardI=Eigen::MatrixXi(numC,1); 447 448 if ((Ndir==4)&&(SMethod==SMNPoly)) 449 HardD=Eigen::MatrixXd(numC,6); 450 else 451 HardD=Eigen::MatrixXd(numC,3); 452 453 for (int i=0;i<numC;i++) 454 { 455 CoordType Dir; 456 int indexF; 457 PickRandomDir(mesh,indexF,Dir); 458 459 HardI(i,0)=indexF; 460 HardD(i,0)=Dir.X(); 461 HardD(i,1)=Dir.Y(); 462 HardD(i,2)=Dir.Z(); 463 464 if ((Ndir==4)&&(SMethod==SMNPoly)) 465 { 466 CoordType Dir1=mesh.face[indexF].N()^Dir; 467 Dir1.Normalize(); 468 HardD(i,3)=Dir1.X(); 469 HardD(i,4)=Dir1.Y(); 470 HardD(i,5)=Dir1.Z(); 471 } 472 } 473 } 474 475 //finally smooth 476 if (SMethod==SMMiq) 477 SmoothMIQ(mesh,HardI,HardD,SoftI,SoftD,SoftW,alphaSoft,Ndir); 478 else 479 { 480 assert(SMethod==SMNPoly); 481 SmoothNPoly(mesh,HardI,HardD,Ndir); 482 } 483 } 484 485 SmoothDirections(MeshType & mesh,SmoothParam SParam)486 static void SmoothDirections(MeshType &mesh,SmoothParam SParam) 487 { 488 //for the moment only cross and line field 489 490 //initialize direction by curvature if needed 491 if ((SParam.alpha_curv>0)|| 492 (SParam.sharp_thr>0)|| 493 (SParam.curv_thr>0)) 494 InitByCurvature(mesh,SParam.curvRing); 495 496 SelectConstraints(mesh,SParam); 497 //then do the actual smooth 498 SmoothDirections(mesh,SParam.Ndir,SParam.SmoothM,true,SParam.alpha_curv); 499 } 500 501 }; 502 503 } // end namespace tri 504 } // end namespace vcg 505 #endif // SMOOTHER_FIELD_H 506