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