1 #ifndef DUAL_OPTIMIZER 2 #define DUAL_OPTIMIZER 3 #include <wrap/callback.h> 4 #ifndef IMPLICIT 5 #include "texcoord_optimization.h" 6 #else 7 #include <poisson_solver.h> 8 #endif 9 template <class MeshType> 10 class BaryOptimizatorDual 11 { 12 typedef typename MeshType::VertexType VertexType; 13 typedef typename MeshType::FaceType FaceType; 14 typedef typename MeshType::CoordType CoordType; 15 typedef typename MeshType::ScalarType ScalarType; 16 #ifndef IMPLICIT 17 typedef typename vcg::tri::AreaPreservingTexCoordOptimization<MeshType> OptType; 18 typedef typename vcg::tri::MeanValueTexCoordOptimization<MeshType> OptType1; 19 #else 20 typedef typename PoissonSolver<MeshType> OptType; 21 /*typedef typename PoissonSolver<MeshType> OptType1;*/ 22 #endif 23 //typedef typename vcg::tri::MeanValueTexCoordOptimization<MeshType> OptType; 24 25 EnergyType EType; 26 27 public: 28 struct param_domain{ 29 MeshType *domain; 30 std::vector<FaceType*> ordered_faces; 31 }; 32 33 //OptType *optimizer; 34 35 ///set of star meshes to optimize barycentryc coords locally 36 std::vector<param_domain> star_meshes; 37 std::vector<param_domain> diamond_meshes; 38 std::vector<param_domain> face_meshes; 39 ///structures for the optimization 40 std::vector<MeshType*> HRES_meshes; 41 std::vector<std::vector<VertexType*> > Ord_HVert; 42 43 ///high resolution mesh and domain mesh 44 MeshType *domain; 45 MeshType *h_res_mesh; 46 47 ///initialize star parametrization InitStarEquilateral()48 void InitStarEquilateral()//const ScalarType &average_area=1) 49 { 50 ///for each vertex 51 int index=0; 52 for (unsigned int i=0;i<domain->vert.size();i++) 53 { 54 if (!(domain->vert[i].IsD())) 55 { 56 std::vector<VertexType*> starCenter; 57 starCenter.push_back(&domain->vert[i]); 58 59 star_meshes[index].domain=new MeshType(); 60 61 ///create star 62 CreateMeshVertexStar(starCenter,star_meshes[index].ordered_faces,*star_meshes[index].domain); 63 64 ///and parametrize it 65 ParametrizeStarEquilateral<MeshType>(*star_meshes[index].domain,1.0); 66 67 index++; 68 } 69 } 70 } 71 72 void InitDiamondEquilateral(const ScalarType &edge_len=1.0) 73 { 74 75 ///for each vertex 76 int index=0; 77 for (unsigned int i=0;i<domain->face.size();i++) 78 { 79 if (!(domain->face[i].IsD())) 80 { 81 FaceType *f0=&domain->face[i]; 82 //for each edge 83 for (int j=0;j<3;j++) 84 { 85 FaceType * f1=f0->FFp(j); 86 if (f1<f0) 87 { 88 ///add to domain map 89 ///end domain mapping 90 91 int num0=j; 92 int num1=f0->FFi(j); 93 94 ///copy the mesh 95 std::vector<FaceType*> faces; 96 faces.push_back(f0); 97 faces.push_back(f1); 98 99 diamond_meshes[index].domain=new MeshType(); 100 101 ///create a copy of the mesh 102 std::vector<VertexType*> orderedVertex; 103 CopyMeshFromFaces<MeshType>(faces,orderedVertex,*diamond_meshes[index].domain); 104 UpdateTopologies<MeshType>(diamond_meshes[index].domain); 105 106 ///set other components 107 diamond_meshes[index].ordered_faces.resize(2); 108 diamond_meshes[index].ordered_faces[0]=f0; 109 diamond_meshes[index].ordered_faces[1]=f1; 110 ///parametrize locally 111 ParametrizeDiamondEquilateral<MeshType>(*diamond_meshes[index].domain,num0,num1,edge_len); 112 113 index++; 114 } 115 } 116 } 117 } 118 } 119 120 121 void InitFaceEquilateral(const ScalarType &edge_len=1) 122 { 123 ///for each vertex 124 int index=0; 125 for (unsigned int i=0;i<domain->face.size();i++) 126 { 127 if (!(domain->face[i].IsD())) 128 { 129 FaceType *f0=&domain->face[i]; 130 131 std::vector<FaceType*> faces; 132 faces.push_back(f0); 133 134 ///create the mesh 135 face_meshes[index].domain=new MeshType(); 136 std::vector<VertexType*> orderedVertex; 137 CopyMeshFromFaces<MeshType>(faces,orderedVertex,*face_meshes[index].domain); 138 139 assert(face_meshes[index].domain->vn==3); 140 assert(face_meshes[index].domain->fn==1); 141 142 ///initialize auxiliary structures 143 face_meshes[index].ordered_faces.resize(1); 144 face_meshes[index].ordered_faces[0]=f0; 145 146 ///parametrize it 147 ParametrizeFaceEquilateral<MeshType>(*face_meshes[index].domain,edge_len); 148 ///add to search structures 149 /*faceMap.insert(std::pair<FaceType*,param_domain*>(f0,&face_meshes[index]));*/ 150 index++; 151 } 152 } 153 } 154 155 ///given a point and a face return the half-star in witch it falls getVertexStar(const CoordType & point,FaceType * f)156 int getVertexStar(const CoordType &point,FaceType *f) 157 { 158 CoordType edge0=(f->P(0)+f->P(1))/2.0; 159 CoordType edge1=(f->P(1)+f->P(2))/2.0; 160 CoordType edge2=(f->P(2)+f->P(0))/2.0; 161 CoordType Center=(f->P(0)+f->P(1)+f->P(2))/3.0; 162 CoordType vect0=edge0-point; 163 CoordType vect1=Center-point; 164 CoordType vect2=edge2-point; 165 CoordType Norm=f->N(); 166 ScalarType in0=(vect0^vect1)*Norm; 167 ScalarType in1=(vect1^vect2)*Norm; 168 if ((in0>=0)&&(in1>=0)) 169 return 0; 170 171 vect0=edge0-point; 172 vect1=Center-point; 173 vect2=edge1-point; 174 in0=(vect1^vect0)*Norm; 175 in1=(vect2^vect1)*Norm; 176 if ((in0>=0)&&(in1>=0)) 177 return 1; 178 179 vect0=edge1-point; 180 vect1=Center-point; 181 vect2=edge2-point; 182 in0=(vect1^vect0)*Norm; 183 in1=(vect2^vect1)*Norm; 184 assert((in0>=0)&&(in1>=0)); 185 return 2; 186 } 187 188 ///given a point and a face return the half-diamond edge index in witch it falls getEdgeDiamond(const CoordType & point,FaceType * f)189 int getEdgeDiamond(const CoordType &point,FaceType *f) 190 { 191 CoordType Center=(f->P(0)+f->P(1)+f->P(2))/3.0; 192 CoordType vect0=f->P(1)-point; 193 CoordType vect1=Center-point; 194 CoordType vect2=f->P(0)-point; 195 CoordType Norm=f->N(); 196 ScalarType in0=(vect0^vect1)*Norm; 197 ScalarType in1=(vect1^vect2)*Norm; 198 if ((in0>=0)&&(in1>=0)) 199 return 0; 200 201 vect0=f->P(2)-point; 202 vect1=Center-point; 203 vect2=f->P(1)-point; 204 in0=(vect0^vect1)*Norm; 205 in1=(vect1^vect2)*Norm; 206 if ((in0>=0)&&(in1>=0)) 207 return 1; 208 209 vect0=f->P(0)-point; 210 vect1=Center-point; 211 vect2=f->P(2)-point; 212 in0=(vect0^vect1)*Norm; 213 in1=(vect1^vect2)*Norm; 214 assert((in0>=0)&&(in1>=0)); 215 return 2; 216 } 217 218 219 ///initialize Star Submeshes InitStarSubdivision()220 void InitStarSubdivision() 221 { 222 223 int index=0; 224 HRES_meshes.clear(); 225 Ord_HVert.clear(); 226 ///initialilze vector of meshes 227 HRES_meshes.resize(star_meshes.size()); 228 Ord_HVert.resize(star_meshes.size()); 229 /*HVert.resize(star_meshes.size());*/ 230 for (unsigned int i=0;i<HRES_meshes.size();i++) 231 HRES_meshes[i]=new MeshType(); 232 233 ///for each vertex of base domain 234 for (unsigned int i=0;i<domain->vert.size();i++) 235 { 236 VertexType *center=&domain->vert[i]; 237 if (!center->IsD()) 238 { 239 ///copy current parametrization of star 240 for (unsigned int k=0;k<star_meshes[index].ordered_faces.size();k++) 241 { 242 FaceType *param=&star_meshes[index].domain->face[k]; 243 FaceType *original=star_meshes[index].ordered_faces[k]; 244 for (int v=0;v<3;v++) 245 original->V(v)->T().P()=param->V(v)->T().P(); 246 } 247 248 ///get h res vertex on faces composing the star 249 std::vector<VertexType*> Hres,inDomain; 250 getHresVertex<FaceType>(star_meshes[index].ordered_faces,Hres); 251 252 ///find out the vertices falling in the substar 253 /*HVert[index].reserve(Hres.size()/2);*/ 254 for (unsigned int k=0;k<Hres.size();k++) 255 { 256 VertexType* chosen; 257 VertexType* test=Hres[k]; 258 CoordType proj=Warp(test); 259 FaceType * father=test->father; 260 CoordType bary=test->Bary; 261 ///get index of half-star 262 int index=getVertexStar(proj,father); 263 chosen=father->V(index); 264 ///if is part of current half star 265 if (chosen==center) 266 { 267 inDomain.push_back(test); 268 ///parametrize it 269 InterpolateUV<MeshType>(father,bary,test->T().U(),test->T().V()); 270 } 271 } 272 ///create Hres mesh already parametrized 273 std::vector<FaceType*> OrderedFaces; 274 CopyMeshFromVertices<MeshType>(inDomain,Ord_HVert[index],OrderedFaces,*HRES_meshes[index]); 275 index++; 276 } 277 } 278 } 279 280 ///initialize Star Submeshes InitDiamondSubdivision()281 void InitDiamondSubdivision() 282 { 283 284 int index=0; 285 HRES_meshes.clear(); 286 Ord_HVert.clear(); 287 ///initialilze vector of meshes 288 HRES_meshes.resize(diamond_meshes.size()); 289 Ord_HVert.resize(diamond_meshes.size()); 290 /*HVert.resize(star_meshes.size());*/ 291 for (unsigned int i=0;i<HRES_meshes.size();i++) 292 HRES_meshes[i]=new MeshType(); 293 294 ///for each edge of base domain 295 for (unsigned int i=0;i<domain->face.size();i++) 296 { 297 FaceType *f0=&domain->face[i]; 298 if (f0->IsD()) 299 break; 300 //for each edge 301 for (int eNum=0;eNum<3;eNum++) 302 { 303 FaceType * f1=f0->FFp(eNum); 304 if (f1<f0) 305 { 306 ///copy current parametrization of diamond 307 for (unsigned int k=0;k<diamond_meshes[index].ordered_faces.size();k++) 308 { 309 FaceType *param=&diamond_meshes[index].domain->face[k]; 310 FaceType *original=diamond_meshes[index].ordered_faces[k]; 311 for (int v=0;v<3;v++) 312 original->V(v)->T().P()=param->V(v)->T().P(); 313 } 314 315 ///get h res vertex on faces composing the diamond 316 std::vector<VertexType*> Hres,inDomain; 317 getHresVertex<FaceType>(diamond_meshes[index].ordered_faces,Hres); 318 319 ///find out the vertices falling in the half-diamond 320 /*HVert[index].reserve(Hres.size()/2);*/ 321 for (unsigned int k=0;k<Hres.size();k++) 322 { 323 //VertexType* chosen; 324 VertexType* test=Hres[k]; 325 CoordType proj=Warp(test); 326 FaceType * father=test->father; 327 CoordType bary=test->Bary; 328 ///get index of half-star 329 int index=getEdgeDiamond(proj,father); 330 ///if is part of current half star 331 if (index==eNum) 332 { 333 inDomain.push_back(test); 334 ///parametrize it 335 InterpolateUV<MeshType>(father,bary,test->T().U(),test->T().V()); 336 } 337 } 338 ///create Hres mesh already parametrized 339 std::vector<FaceType*> OrderedFaces; 340 CopyMeshFromVertices<MeshType>(inDomain,Ord_HVert[index],OrderedFaces,*HRES_meshes[index]); 341 index++; 342 } 343 } 344 } 345 } 346 347 ///initialize Star Submeshes InitFaceSubdivision()348 void InitFaceSubdivision() 349 { 350 351 int index=0; 352 HRES_meshes.clear(); 353 Ord_HVert.clear(); 354 ///initialilze vector of meshes 355 HRES_meshes.resize(face_meshes.size()); 356 Ord_HVert.resize(face_meshes.size()); 357 358 for (unsigned int i=0;i<HRES_meshes.size();i++) 359 HRES_meshes[i]=new MeshType(); 360 361 ///for each face of base domain 362 for (unsigned int i=0;i<domain->face.size();i++) 363 { 364 FaceType *f0=&domain->face[i]; 365 if (f0->IsD()) 366 break; 367 ///copy current parametrization of face 368 FaceType *param=&face_meshes[index].domain->face[0]; 369 FaceType *original=face_meshes[index].ordered_faces[0]; 370 371 assert(face_meshes[index].domain->vn==3); 372 assert(face_meshes[index].domain->fn==1); 373 assert(face_meshes[index].ordered_faces.size()==1); 374 assert(original==f0); 375 376 for (int v=0;v<3;v++) 377 original->V(v)->T().P()=param->V(v)->T().P(); 378 379 ///get h res vertex on faces composing the diamond 380 std::vector<VertexType*> inDomain; 381 getHresVertex<FaceType>(face_meshes[index].ordered_faces,inDomain); 382 383 ///transform in UV 384 for (unsigned int k=0;k<inDomain.size();k++) 385 { 386 VertexType* test=inDomain[k]; 387 FaceType * father=test->father; 388 assert(father==f0); 389 CoordType bary=test->Bary; 390 InterpolateUV<MeshType>(father,bary,test->T().U(),test->T().V()); 391 } 392 ///create Hres mesh already parametrized 393 std::vector<FaceType*> OrderedFaces; 394 CopyMeshFromVertices<MeshType>(inDomain,Ord_HVert[index],OrderedFaces,*HRES_meshes[index]); 395 index++; 396 } 397 } 398 399 MinimizeStep(const int & phaseNum)400 void MinimizeStep(const int &phaseNum) 401 { 402 403 404 //Ord_HVert[index] 405 for (unsigned int i=0;i<HRES_meshes.size();i++) 406 { 407 408 MeshType *currMesh=HRES_meshes[i]; 409 if (currMesh->fn>0) 410 { 411 UpdateTopologies<MeshType>(currMesh); 412 413 ///on star 414 int numDom=1; 415 switch (phaseNum) 416 { 417 case 0:numDom=6;break;//star 418 case 1:numDom=2;break;//diam 419 case 2:numDom=1;break;//face 420 } 421 ///save previous values 422 #ifndef IMPLICIT 423 InitDampRestUV(*currMesh); 424 bool b=UnFold<MeshType>(*currMesh,numDom); 425 bool isOK=testParamCoords<MeshType>(*currMesh); 426 427 if ((!b)||(!isOK)) 428 RestoreRestUV<MeshType>(*currMesh); 429 ///save previous values 430 InitDampRestUV(*currMesh); 431 432 433 434 ///NEW SETTING SPEED 435 436 ScalarType edge_esteem=GetSmallestUVHeight(*currMesh); 437 438 439 ScalarType speed0=edge_esteem*0.2; 440 ScalarType conv=edge_esteem*0.01; 441 442 if (accuracy>1) 443 conv*=1.0/(ScalarType)((accuracy-1)*10.0); 444 #endif 445 #ifndef IMPLICIT 446 if (EType==EN_EXTMips) 447 { 448 OptType opt(*currMesh); 449 opt.TargetCurrentGeometry(); 450 opt.SetBorderAsFixed(); 451 opt.SetSpeed(speed0); 452 opt.IterateUntilConvergence(conv); 453 } 454 else 455 if (EType==EN_MeanVal) 456 { 457 OptType1 opt(*currMesh); 458 opt.TargetCurrentGeometry(); 459 opt.SetBorderAsFixed(); 460 opt.SetSpeed(speed0); 461 opt.IterateUntilConvergence(conv); 462 } 463 #else 464 OptType opt(*currMesh); 465 opt.SetBorderAsFixed(); 466 opt.SolvePoisson(); 467 #endif 468 //opt.IterateUntilConvergence(); 469 470 ///test for uv errors 471 //bool IsOK=true; 472 for (unsigned int j=0;j<currMesh->vert.size();j++) 473 { 474 VertexType *ParamVert=&currMesh->vert[j]; 475 ScalarType u=ParamVert->T().U(); 476 ScalarType v=ParamVert->T().V(); 477 if ((!((u<=1.001)&&(u>=-1.001)))|| 478 (!(v<=1.001)&&(v>=-1.001))) 479 { 480 //IsOK=false; 481 482 for (unsigned int k=0;k<currMesh->vert.size();k++) 483 currMesh->vert[k].T().P()=currMesh->vert[k].RestUV; 484 break; 485 } 486 } 487 //reassing fathers and bary coordinates 488 for (unsigned int j=0;j<currMesh->vert.size();j++) 489 { 490 VertexType *ParamVert=&currMesh->vert[j]; 491 VertexType *OrigVert=Ord_HVert[i][j]; 492 ScalarType u=ParamVert->T().U(); 493 ScalarType v=ParamVert->T().V(); 494 ///then get face falling into and estimate (alpha,beta,gamma) 495 CoordType bary; 496 BaseFace* chosen; 497 param_domain *currDom; 498 switch (phaseNum) 499 { 500 case 0:currDom=&star_meshes[i];break;//star 501 case 1:currDom=&diamond_meshes[i];break;//diam 502 case 2:currDom=&face_meshes[i];break;//face 503 } 504 /*assert(currDom->domain->vn==3); 505 assert(currDom->domain->fn==1);*/ 506 bool inside=GetBaryFaceFromUV(*currDom->domain,u,v,currDom->ordered_faces,bary,chosen); 507 if (!inside) 508 { 509 /*#ifndef _MESHLAB*/ 510 printf("\n OUTSIDE %f,%f \n",u,v); 511 /*#endif*/ 512 vcg::Point2<ScalarType> UV=vcg::Point2<ScalarType>(u,v); 513 ForceInParam<MeshType>(UV,*currDom->domain); 514 u=UV.X(); 515 v=UV.Y(); 516 inside=GetBaryFaceFromUV(*currDom->domain,u,v,currDom->ordered_faces,bary,chosen); 517 //assert(0); 518 } 519 assert(inside); 520 //OrigVert->father=chosen; 521 //OrigVert->Bary=bary; 522 AssingFather(*OrigVert,chosen,bary,*domain); 523 } 524 } 525 ///delete current mesh 526 delete(HRES_meshes[i]); 527 } 528 529 ///clear father and bary 530 for (unsigned int i=0;i<domain->face.size();i++) 531 domain->face[i].vertices_bary.clear(); 532 533 ///set face-vertex link 534 for (unsigned int i=0;i<h_res_mesh->vert.size();i++) 535 { 536 BaseVertex *v=&h_res_mesh->vert[i]; 537 if (!v->IsD()) 538 { 539 BaseFace *f=v->father; 540 CoordType bary=v->Bary; 541 f->vertices_bary.push_back(std::pair<VertexType*,CoordType>(v,bary)); 542 } 543 } 544 } 545 546 547 int accuracy; 548 vcg::CallBackPos *cb; 549 int step; 550 551 public: 552 553 554 void Init(MeshType &_domain, 555 MeshType &_h_res_mesh, 556 vcg::CallBackPos *_cb, 557 int _accuracy=1, 558 EnergyType _EType=EN_EXTMips) 559 { 560 EType=_EType; 561 562 step=0; 563 cb=_cb; 564 accuracy=_accuracy; 565 566 vcg::tri::UpdateNormal<MeshType>::PerFaceNormalized(_domain); 567 568 domain=&_domain; 569 h_res_mesh=&_h_res_mesh; 570 571 ///initialize STARS 572 star_meshes.resize(domain->vn); 573 InitStarEquilateral(); 574 /*InitStarSubdivision();*/ 575 576 ///initialize DIAMONDS 577 int num_edges=0; 578 for (unsigned int i=0;i<domain->face.size();i++) 579 { 580 if (!(domain->face[i].IsD())) 581 { 582 FaceType *f0=&domain->face[i]; 583 //for each edge 584 for (int j=0;j<3;j++) 585 { 586 FaceType * f1=f0->FFp(j); 587 if (f1<f0) 588 num_edges++; 589 } 590 } 591 } 592 diamond_meshes.resize(num_edges); 593 InitDiamondEquilateral(); 594 595 ///initialize FACES 596 face_meshes.resize(domain->fn); 597 InitFaceEquilateral(); 598 599 ///init minimizer 600 for (unsigned int i=0;i<h_res_mesh->vert.size();i++) 601 h_res_mesh->vert[i].P()=h_res_mesh->vert[i].RPos; 602 603 //InitDampRestUV(*h_res_mesh); 604 } 605 606 PrintAttributes()607 void PrintAttributes() 608 { 609 610 int done=step; 611 int total=6; 612 ScalarType ratio=(ScalarType)done/total; 613 int percent=(int)(ratio*(ScalarType)100); 614 ScalarType distArea=ApproxAreaDistortion<BaseMesh>(*h_res_mesh,domain->fn); 615 ScalarType distAngle=ApproxAngleDistortion<BaseMesh>(*h_res_mesh); 616 char ret[200]; 617 sprintf(ret," PERFORM GLOBAL OPTIMIZATION Area distorsion:%4f ; ANGLE distorsion:%4f ",distArea,distAngle); 618 (*cb)(percent,ret); 619 } 620 621 void Optimize(ScalarType gap=0.5,int max_step=10) 622 { 623 int k=0; 624 ScalarType distArea=ApproxAreaDistortion<BaseMesh>(*h_res_mesh,domain->fn); 625 ScalarType distAngle=ApproxAngleDistortion<BaseMesh>(*h_res_mesh); 626 ScalarType distAggregate0=geomAverage<ScalarType>(distArea+1.0,distAngle+1.0,3,1)-1; 627 bool ContinueOpt=true; 628 PatchesOptimizer<BaseMesh> DomOpt(*domain,*h_res_mesh); 629 step++; 630 631 #ifndef IMPLICIT 632 DomOpt.OptimizePatches(); 633 PrintAttributes(); 634 #endif 635 while (ContinueOpt) 636 { 637 ///domain Optimization 638 k++; 639 640 InitStarSubdivision(); 641 MinimizeStep(0); 642 643 InitDiamondSubdivision(); 644 MinimizeStep(1); 645 646 InitFaceSubdivision(); 647 MinimizeStep(2); 648 step++; 649 PrintAttributes(); 650 651 distArea=ApproxAreaDistortion<BaseMesh>(*h_res_mesh,domain->fn); 652 distAngle=ApproxAngleDistortion<BaseMesh>(*h_res_mesh); 653 ScalarType distAggregate1=geomAverage<ScalarType>(distArea+1.0,distAngle+1.0,3,1)-1; 654 ScalarType NewGap=((distAggregate0-distAggregate1)*ScalarType(100.0))/distAggregate0; 655 if ((NewGap<gap)||(k>max_step)) 656 ContinueOpt=false; 657 distAggregate0=distAggregate1; 658 } 659 } 660 }; 661 #endif 662