1 #ifndef LOCAL_PARAMETRIZATION
2 #define LOCAL_PARAMETRIZATION
3 
4 
5 
6 #include "defines.h"
7 ///fitting
8 //#include <vcg/space/fitting3.h>
9 #include <vcg/math/matrix33.h>
10 #include <vcg/space/triangle2.h>
11 #include "texcoord_optimization.h"
12 #include "mesh_operators.h"
13 
14 //#include <vcg/complex/algorithms/point_sampling.h>
15 
16 //#define samples_area 80
17 
18 template <class MeshType>
ParametrizeExternal(MeshType & to_parametrize)19 void ParametrizeExternal(MeshType &to_parametrize)
20 {
21     typedef typename MeshType::ScalarType ScalarType;
22     typedef typename MeshType::VertexType VertexType;
23 
24         std::vector<VertexType*> vertices;
25 
26     ///find first border vertex
27         VertexType* Start=NULL;
28         typename MeshType::VertexIterator Vi=to_parametrize.vert.begin();
29 
30     while ((Start==NULL)&&(Vi<to_parametrize.vert.end()))
31     {
32         if (((*Vi).IsB())&&(!(*Vi).IsD()))
33             Start=&(*Vi);
34         Vi++;
35     }
36     if (Vi==to_parametrize.vert.end())
37     {
38         assert(0);
39     }
40     ///get sorted border vertices
41 
42     FindSortedBorderVertices<MeshType>(to_parametrize,Start,vertices);
43 
44     //assert(vertices.size()>=4);
45 
46     ///find perimeter
47     ScalarType perimeter=0;
48     int size=vertices.size();
49     for (int i=0;i<size;i++)
50         perimeter+=(vertices[i]->P()-vertices[(i+1)%size]->P()).Norm();
51 
52     ///find scaling factor
53     /*ScalarType Sperimeter=(2.0*M_PI)/perimeter;*/
54 
55     ///set default texCoords
56     for (Vi=to_parametrize.vert.begin();Vi!=to_parametrize.vert.end();Vi++)
57     {
58         (*Vi).T().U()=-2;
59         (*Vi).T().V()=-2;
60     }
61 
62     ///set border vertices
63         typename std::vector<VertexType*>::iterator iteV;
64     /*ScalarType curr_perim=0;*/
65     ScalarType curr_angle=0;
66     vertices[0]->T().U()=cos(curr_angle);
67     vertices[0]->T().V()=sin(curr_angle);
68     //for (int i=1;i<vertices.size();i++)
69     //{
70     //	curr_perim+=(vertices[i]->P()-vertices[(i-1)]->P()).Norm();
71     //	//curr_perim+=perimeter/(ScalarType)size;
72     //	curr_angle=curr_perim*Sperimeter;
73     //	vertices[i]->T().U()=cos(curr_angle);
74     //	vertices[i]->T().V()=sin(curr_angle);
75     //	assert((vertices[i]->T().U()>=-1)&&(vertices[i]->T().U()<=1));
76     //	assert((vertices[i]->T().V()>=-1)&&(vertices[i]->T().V()<=1));
77     //}
78     ScalarType anglediv=(2.0*M_PI)/(ScalarType)(vertices.size());
79     curr_angle=0;
80     for (unsigned int i=1;i<vertices.size();i++)
81     {
82         curr_angle+=anglediv;
83         vertices[i]->T().U()=cos(curr_angle);
84         vertices[i]->T().V()=sin(curr_angle);
85         assert((vertices[i]->T().U()>=-1)&&(vertices[i]->T().U()<=1));
86         assert((vertices[i]->T().V()>=-1)&&(vertices[i]->T().V()<=1));
87     }
88 }
89 
90 template <class MeshType>
ParametrizeInternal(MeshType & to_parametrize)91 void ParametrizeInternal(MeshType &to_parametrize)
92 {
93     typedef typename MeshType::ScalarType ScalarType;
94     const ScalarType Eps=(ScalarType)0.0001;
95     ///set internal vertices
96         for (typename MeshType::VertexIterator Vi=to_parametrize.vert.begin();Vi!=to_parametrize.vert.end();Vi++)
97     {
98         //assert(!Vi->IsD());
99         if ((!Vi->IsB())&&(!Vi->IsD()))
100         {
101             ///find kernel
102                         std::vector<typename MeshType::VertexType*> star;
103             getVertexStar<MeshType>(&(*Vi),star);
104             ScalarType kernel=0;
105             for (unsigned int k=0;k<star.size();k++)
106                 if (star[k]->IsB())
107                 {
108                     ScalarType dist=((*Vi).P()-star[k]->P()).Norm();
109                     if (dist<Eps)
110                         dist=Eps;
111                     //ScalarType peso=1.0/(dist);
112                     ScalarType peso=(dist)/star.size();
113                     kernel+=(peso);//*dist);
114                 }
115             assert(kernel>0);
116             ///then find factor
117             kernel=1.0/kernel;
118 
119             (*Vi).T().U()=0;
120             (*Vi).T().V()=0;
121 
122             int num=0;
123             ///find weighted media
124             for (unsigned int k=0;k<star.size();k++)
125                 if (star[k]->IsB())
126                 {
127                     ScalarType dist=((*Vi).P()-star[k]->P()).Norm();
128                     if (dist<Eps)
129                         dist=Eps;
130                     //ScalarType peso=1.0/(dist);
131                     ScalarType peso=(dist)/star.size();
132                     ScalarType kval=(peso)*kernel;
133                     assert(kval>0);
134                     (*Vi).T().U()+=kval*star[k]->T().U();
135                     (*Vi).T().V()+=kval*star[k]->T().V();
136                     num++;
137                 }
138             ////on border case 2 neighbors
139             ///go next to the center
140             /*if (num==2)
141             {
142                 (*Vi).T().U()/=2.0;
143                 (*Vi).T().V()/=2.0;
144             }*/
145             /*ScalarType u=(*Vi).T().U();
146             ScalarType v=(*Vi).T().V();*/
147             assert(((*Vi).T().U()>=-1)&&((*Vi).T().U()<=1));
148             assert(((*Vi).T().V()>=-1)&&((*Vi).T().V()<=1));
149         }
150     }
151     ///smoothing of txcoords
152     InitDampRestUV(to_parametrize);
153         for (typename MeshType::VertexIterator Vi=to_parametrize.vert.begin();Vi!=to_parametrize.vert.end();Vi++)
154     {
155         if ((!Vi->IsB())&&(!Vi->IsD()))
156         {
157                         std::vector<typename MeshType::VertexType*> star;
158             getVertexStar<MeshType>(&(*Vi),star);
159             vcg::Point2<ScalarType> UV=vcg::Point2<ScalarType>(0,0);
160             for (unsigned int k=0;k<star.size();k++)
161                 UV+=star[k]->RestUV;
162             UV/=(ScalarType)star.size();
163             (*Vi).T().P()=UV;
164         }
165     }
166 }
167 
168 
169 template <class FaceType>
InterpolatePos(FaceType * f,const typename FaceType::CoordType & bary)170 typename FaceType::CoordType InterpolatePos
171 (FaceType* f, const typename FaceType::CoordType &bary)
172 {return (f->V(0)->P()*bary.X()+f->V(1)->P()*bary.Y()+f->V(2)->P()*bary.Z());}
173 
174 template <class FaceType>
InterpolateRPos(FaceType * f,const typename FaceType::CoordType & bary)175 typename FaceType::CoordType InterpolateRPos
176 (FaceType* f,const typename FaceType::CoordType &bary)
177 {
178     return (f->V(0)->RPos*bary.X()+f->V(1)->RPos*bary.Y()+f->V(2)->RPos*bary.Z());
179 }
180 
181 template <class FaceType>
InterpolateNorm(FaceType * f,const typename FaceType::CoordType & bary)182 typename FaceType::CoordType InterpolateNorm
183 (FaceType* f, const typename FaceType::CoordType &bary)
184 {
185         typedef typename FaceType::CoordType CoordType;
186     CoordType n0=f->V(0)->N();
187     CoordType n1=f->V(1)->N();
188     CoordType n2=f->V(2)->N();
189     return (n0*bary.X()+n1*bary.Y()+n2*bary.Z());
190 }
191 
192 template <class ScalarType>
Approx(const ScalarType & value)193 int Approx(const ScalarType &value)
194 {
195     ScalarType val0=floor(value);
196     ScalarType val1=ceil(value);
197     if (fabs(val0-value)<fabs(val1-value))
198         return ((int)val0);
199     else
200         return ((int)val1);
201 }
202 
203 template <class FaceType>
InterpolateColor(FaceType * f,const typename FaceType::CoordType & bary)204 vcg::Point3i InterpolateColor
205 (FaceType* f,const typename FaceType::CoordType &bary)
206 {
207         typedef typename FaceType::ScalarType ScalarType;
208     vcg::Color4b c0=f->V(0)->C();
209     vcg::Color4b c1=f->V(1)->C();
210     vcg::Color4b c2=f->V(2)->C();
211     double R=(ScalarType)c0.X()*bary.X()+(ScalarType)c1.X()*bary.Y()+(ScalarType)c2.X()*bary.Z();
212     double G=(ScalarType)c0.Y()*bary.X()+(ScalarType)c1.Y()*bary.Y()+(ScalarType)c2.Y()*bary.Z();
213     double B=(ScalarType)c0.Z()*bary.X()+(ScalarType)c1.Z()*bary.Y()+(ScalarType)c2.Z()*bary.Z();
214 
215     vcg::Point3i p=vcg::Point3i(Approx(R),Approx(G),Approx(B));
216 
217     assert(p.X()<=255);
218     assert(p.Y()<=255);
219     assert(p.Z()<=255);
220 
221     return (p);
222 }
223 
224 template <class VertexType>
ProjectPos(const VertexType & v)225 typename VertexType::CoordType ProjectPos(const VertexType &v)
226 {
227     typedef typename VertexType::FaceType FaceType;
228     typedef typename VertexType::CoordType CoordType;
229 
230     FaceType *f=v.father;
231     CoordType b=v.Bary;
232     return (InterpolatePos<FaceType>(f,b));
233 }
234 
235 template <class VertexType>
Warp(const VertexType * v)236 typename VertexType::CoordType Warp(const VertexType* v)
237 {
238         typename VertexType::FaceType * father=v->father;
239         typename VertexType::CoordType proj=father->P(0)*v->Bary.X()+father->P(1)*v->Bary.Y()+father->P(2)*v->Bary.Z();
240     return proj;
241 }
242 
243 template <class VertexType>
WarpRpos(const VertexType * v)244 typename VertexType::CoordType WarpRpos(const VertexType* v)
245 {
246         typename VertexType::FaceType * father=v->father;
247         typename VertexType::CoordType proj=father->V(0)->RPos*v->Bary.X()+father->V(1)->RPos*v->Bary.Y()+father->V(2)->RPos*v->Bary.Z();
248     return proj;
249 }
250 
251 template <class MeshType>
EstimateLengthByParam(const typename MeshType::VertexType * v0,const typename MeshType::VertexType * v1,typename MeshType::FaceType * on_edge[2])252 typename MeshType::ScalarType EstimateLengthByParam
253                                 (const typename MeshType::VertexType* v0,
254                                  const typename MeshType::VertexType* v1,
255                                  typename MeshType::FaceType* on_edge[2])
256 {
257     typedef typename MeshType::FaceType FaceType;
258     typedef typename MeshType::CoordType CoordType;
259     typedef typename MeshType::VertexType VertexType;
260     typedef typename MeshType::ScalarType ScalarType;
261 
262 //	assert((on_edge.size()==0)||(on_edge.size()==1));
263     ScalarType estimated[2]={0,0};
264     size_t num[2]={0,0};
265 
266     for (int i=0;i<2;i++)
267     {
268         FaceType *test_face=on_edge[i];
269 
270         int edge_index=EdgeIndex(test_face,v0,v1);
271         FaceType *Fopp=test_face->FFp(edge_index);
272 
273         if (test_face->vertices_bary.size()<2)
274         {
275             ScalarType dist=Distance(v0->RPos,v1->RPos);
276             //#pragma omp atomic
277             estimated[i]+=dist;
278             num[i]=0;
279             continue;
280         }
281 
282         ///collect vertices
283         std::vector<VertexType*> vertices;
284         vertices.reserve(test_face->vertices_bary.size());
285         for (unsigned int k=0;k<test_face->vertices_bary.size();k++)
286             vertices.push_back(test_face->vertices_bary[k].first);
287 
288 
289         ///collect faces
290         std::vector<FaceType*> faces;
291         getSharedFace<MeshType>(vertices,faces);
292 
293         ///get border edges
294         std::vector<std::pair<VertexType*,VertexType*> > edges;
295         for (unsigned int j=0;j<faces.size();j++)
296         {
297             FaceType*f=faces[j];
298             ///find if there is an on-edge edge
299             bool found=false;
300             int k=0;
301             while  ((k<3)&&(!found))
302             {
303                 if ((f->V0(k)->father==test_face)&&
304                     (f->V1(k)->father==test_face)&&
305                     (f->V2(k)->father==Fopp))
306                 {
307                     edges.push_back(std::pair<VertexType*,VertexType*>(f->V0(k),f->V1(k)));
308                     found=true;
309                 }
310                 k++;
311             }
312         }
313 
314         ///find if there's ant edge return initial length
315         if (edges.size()==0)
316         {
317             estimated[i]+=(Distance(v0->RPos,v1->RPos));
318             num[i]=0;
319             continue;
320         }
321         else
322         {
323             //get edge direction
324             ///store the two nearest for each vertex
325             /*VertexType *n0=edges[0].first;
326             VertexType *n1=edges[0].second;
327             ScalarType d0=(Warp(n0)-v0->P()).Norm();
328             ScalarType d1=(Warp(n1)-v1->P()).Norm();*/
329 
330             //CoordType edgedir=v0->cP()-v1->cP();
331             CoordType edgedir=v0->RPos-v1->RPos;
332             edgedir.Normalize();
333             num[i]=edges.size();
334             for (unsigned int e=0;e<edges.size();e++)
335             {
336                 VertexType *vH0=edges[e].first;
337                 VertexType *vH1=edges[e].second;
338 
339                 ///project points over the plane
340                 /*CoordType proj0=Warp(vH0);
341                 CoordType proj1=Warp(vH1);*/
342                 CoordType proj0=WarpRpos(vH0);
343                 CoordType proj1=WarpRpos(vH1);
344 
345                 CoordType dirproj=proj0-proj1;
346                 dirproj.Normalize();
347                 //estimated[i]+=fabs(dirproj*edgedir)*((vH0->P()-vH1->P()).Norm());
348                 //#pragma omp atomic
349                 estimated[i]+=fabs(dirproj*edgedir)*((vH0->RPos-vH1->RPos).Norm());
350 
351             }
352         }
353     }
354 
355     ///media of estimated values
356     ScalarType alpha0,alpha1;
357     ScalarType max_num=abstraction_num;
358 
359     if (num[0]>=max_num)
360         alpha0=1;
361     else
362         alpha0=num[0]/max_num;
363 
364     if (num[1]>=max_num)
365         alpha1=1;
366     else
367         alpha1=num[1]/max_num;
368 
369     estimated[0]= (ScalarType)(alpha0*estimated[0]+(1.0-alpha0)*(Distance(v0->RPos,v1->RPos)));
370     estimated[1]= (ScalarType)(alpha1*estimated[1]+(1.0-alpha1)*(Distance(v0->RPos,v1->RPos)));
371     return(ScalarType)((estimated[0]+estimated[1])/2.0);
372 }
373 
374 template <class MeshType>
MeanVal(const std::vector<vcg::Point2<typename MeshType::ScalarType>> & Points,std::vector<typename MeshType::ScalarType> & Lamda,typename MeshType::CoordType & p)375 void MeanVal(const std::vector<vcg::Point2<typename MeshType::ScalarType> > &Points,
376              std::vector<typename MeshType::ScalarType> &Lamda,
377              typename MeshType::CoordType &p)
378 {
379     typedef typename MeshType::ScalarType ScalarType;
380         typedef typename MeshType::CoordType CoordType;
381 
382     int size=Points.size();
383     Lamda.resize(size);
384 
385     ScalarType sum=0;
386     for (int i=0;i<size;i++)
387     {
388         int size=Points.size()-1;
389         vcg::Point2<ScalarType> Pcurr=Points[i];
390         vcg::Point2<ScalarType> Pprev=Points[(i+(size-1))%size];
391         vcg::Point2<ScalarType> Pnext=Points[(i+1)%size];
392 
393         CoordType v0=Pprev-p;
394         CoordType v1=Pcurr-p;
395         CoordType v2=Pnext-p;
396         ScalarType l=v1.Norm();
397         v0.Normalize();
398         v1.Normalize();
399         v2.Normalize();
400 
401         ScalarType Alpha0=acos(v0*v1);
402         ScalarType Alpha1=acos(v1*v2);
403 
404                 Lamda[i]=(tan(Alpha0/2.0)+tan(Alpha1/2.0))/l;
405                 sum+=Lamda[i];
406     }
407 
408     ///normalization
409     for (int i=0;i<size;i++)
410                 Lamda[i]/=sum;
411 }
412 
413 template <class FaceType>
EstimateAreaByParam(const FaceType * f)414 typename FaceType::ScalarType EstimateAreaByParam(const FaceType* f)
415 {
416     typedef typename FaceType::VertexType VertexType;
417     typedef typename FaceType::ScalarType ScalarType;
418 
419     int num=0;
420     ScalarType estimated=0;
421     for (unsigned int k=0;k<f->vertices_bary.size();k++)
422     {
423         VertexType *HresVert=f->vertices_bary[k].first;
424         estimated+=HresVert->area;
425         num++;
426     }
427 
428     ///media of estimated values
429     ScalarType alpha;
430     ScalarType max_num=abstraction_num;
431 
432     if (num>=max_num)
433         alpha=1;
434     else
435         alpha=num/max_num;
436 
437     ScalarType Rarea= (ScalarType)(((f->cV(1)->RPos-f->cV(0)->RPos)^(f->cV(2)->RPos-f->cV(0)->RPos)).Norm()/2.0);
438     estimated=(ScalarType)(alpha*estimated+(1.0-alpha)*Rarea);
439 
440     return(estimated);
441 }
442 
443 template <class MeshType>
EstimateAreaByParam(const typename MeshType::VertexType * v0,const typename MeshType::VertexType * v1,typename MeshType::FaceType * on_edge[2])444 typename MeshType::ScalarType EstimateAreaByParam
445                               (const typename MeshType::VertexType* v0,
446                                const typename MeshType::VertexType* v1,
447                                typename MeshType::FaceType* on_edge[2])
448 {
449     typedef typename MeshType::FaceType FaceType;
450     typedef typename MeshType::VertexType VertexType;
451     typedef typename MeshType::ScalarType ScalarType;
452 
453     //MeshType::PerVertexAttributeHandle<AuxiliaryVertData> handle = vcg::tri::Allocator<MeshType>::GetPerVertexAttribute<AuxiliaryVertData>(mesh,"AuxiliaryVertData");
454 
455     ScalarType estimated[2]={0,0};
456     int num[2]={0,0};
457     VertexType *v2[2];
458     for (int i=0;i<2;i++)
459     {
460         FaceType *test_face=on_edge[i];
461 
462         for (int k=0;k<3;k++)
463             if ((test_face->V(k)!=v0)&&(test_face->V(k)!=v1))
464                 v2[i]=test_face->V(k);
465 
466         for (unsigned int k=0;k<test_face->vertices_bary.size();k++)
467         {
468             VertexType *brother=test_face->vertices_bary[k].first;
469             estimated[i]+=brother->area;
470             //estimated[i]+=handle[brother].area;
471             num[i]++;
472         }
473     }
474 
475     ///media of estimated values
476     ScalarType alpha0,alpha1;
477     ScalarType max_num=abstraction_num;//20
478 
479     if (num[0]>=max_num)
480         alpha0=1;
481     else
482         alpha0=num[0]/max_num;
483 
484     if (num[1]>=max_num)
485         alpha1=1;
486     else
487         alpha1=num[1]/max_num;
488     ScalarType Rarea0= (ScalarType)(((on_edge[0]->V(1)->RPos-on_edge[0]->V(0)->RPos)^(on_edge[0]->V(2)->RPos-on_edge[0]->V(0)->RPos)).Norm()/2.0);
489     ScalarType Rarea1= (ScalarType)(((on_edge[1]->V(1)->RPos-on_edge[1]->V(0)->RPos)^(on_edge[1]->V(2)->RPos-on_edge[1]->V(0)->RPos)).Norm()/2.0);
490     estimated[0]= (ScalarType)(alpha0*estimated[0]+(1.0-alpha0)*Rarea0);
491     estimated[1]= (ScalarType)(alpha1*estimated[1]+(1.0-alpha1)*Rarea1);
492     return(ScalarType)((estimated[0]+estimated[1])/2.0);
493 }
494 
495 ///template class used to sample surface
496 template <class FaceType>
497 class VertexSampler{
498     typedef typename FaceType::CoordType CoordType;
499 public:
500     std::vector<CoordType> points;
501 
AddFace(const FaceType & f,const CoordType & bary)502     void AddFace(const FaceType &f,const CoordType & bary)
503     {points.push_back(f.P(0)*bary.X()+f.P(1)*bary.Y()+f.P(2)*bary.Z());}
504 };
505 
506 
507 ///sample 3d vertex possible's position
508 ///using area criterion
509 //template <class MeshType>
510 //void SamplingPoints(MeshType &mesh,
511 //					std::vector<typename MeshType::CoordType> &pos)
512 //{
513 //	typedef typename MeshType::CoordType CoordType;
514 //	typedef VertexSampler<MeshType::FaceType> Sampler;
515 //	pos.reserve(samples_area);
516 //	Sampler ps;
517 //	ps.points.reserve(samples_area);
518 //
519 //	vcg::tri::SurfaceSampling<MeshType,Sampler>::Montecarlo(mesh,ps,samples_area);
520 //	pos=std::vector<CoordType>(ps.points.begin(),ps.points.end());
521 //}
522 
523 template <class MeshType>
InitDampRestUV(MeshType & m)524 void InitDampRestUV(MeshType &m)
525 {
526     for (unsigned int i=0;i<m.vert.size();i++)
527         m.vert[i].RestUV=m.vert[i].T().P();
528 }
529 
530 
531 template <class MeshType>
RestoreRestUV(MeshType & m)532 void RestoreRestUV(MeshType &m)
533 {
534     for (unsigned int i=0;i<m.vert.size();i++)
535         m.vert[i].T().P()=m.vert[i].RestUV;
536 }
537 
538 #ifndef IMPLICIT
539 ///parametrize a submesh from trinagles that are incident on vertices
540 template <class MeshType>
541 void ParametrizeLocally(MeshType &parametrized,
542                         bool fix_boundary=true,
543                         bool init_border=true)
544 {
545     typedef typename MeshType::FaceType FaceType;
546     typedef typename MeshType::ScalarType ScalarType;
547     typedef typename MeshType::CoordType CoordType;
548 
549     //const ScalarType epsilon=(ScalarType)0.0001;
550 
551 
552     ///save old positions
553 
554     std::vector<CoordType> positions;
555     positions.resize(parametrized.vert.size());
556     ///set rest position
557     for (unsigned int i=0;i<parametrized.vert.size();i++)
558     {
559         positions[i]=parametrized.vert[i].P();
560         parametrized.vert[i].P()=parametrized.vert[i].RPos;
561     }
562 
563     UpdateTopologies(&parametrized);
564 
565     if (init_border)
566         ParametrizeExternal(parametrized);
567 
568     ParametrizeInternal(parametrized);
569 
570     //for (int i=0;i<parametrized.vert.size();i++)
571     //	parametrized.vert[i].RPos=parametrized.vert[i].P();
572 
573     vcg::tri::MeanValueTexCoordOptimization<MeshType> opt(parametrized);
574     vcg::tri::AreaPreservingTexCoordOptimization<MeshType> opt1(parametrized);
575     opt.SetSpeed((ScalarType)0.0005);
576     InitDampRestUV(parametrized);
577     if (fix_boundary)
578     {
579         opt.TargetEquilateralGeometry();
580         //opt.TargetCurrentGeometry();
581         opt.SetBorderAsFixed();
582         opt.IterateUntilConvergence((ScalarType)0.000001,100);
583         /*opt.Iterate();
584         opt.Iterate();*/
585     }
586     else
587     {
588         opt1.TargetCurrentGeometry();
589         //opt1.SetBorderAsFixed();
590         opt1.IterateUntilConvergence((ScalarType)0.000001,200);
591     }
592 
593     ///assert parametrization
594     for (unsigned int i=0;i<parametrized.face.size();i++)
595     {
596         FaceType *f=&parametrized.face[i];
597         vcg::Point2<ScalarType> tex0=vcg::Point2<ScalarType>(f->V(0)->T().U(),f->V(0)->T().V());
598         vcg::Point2<ScalarType> tex1=vcg::Point2<ScalarType>(f->V(1)->T().U(),f->V(1)->T().V());
599         vcg::Point2<ScalarType> tex2=vcg::Point2<ScalarType>(f->V(2)->T().U(),f->V(2)->T().V());
600 #ifndef NDEBUG
601         ScalarType area=(tex1-tex0)^(tex2-tex0);
602         assert(area>0);
603 #endif
604     }
605     ///restore position
606     for (unsigned int i=0;i<parametrized.vert.size();i++)
607         parametrized.vert[i].P()=positions[i];
608 }
609 #else
610 
611 ///parametrize a submesh keeping fixed the borders
612 template <class MeshType>
613 void ParametrizeLocally(MeshType &parametrized,
614                         bool fix_boundary=true,
615                         bool init_border=true)
616 {
617     typedef typename MeshType::FaceType FaceType;
618     typedef typename MeshType::ScalarType ScalarType;
619     typedef typename MeshType::CoordType CoordType;
620 
621     ///save old positions
622     std::vector<CoordType> positions;
623     positions.resize(parametrized.vert.size());
624     ///set rest position
625     for (unsigned int i=0;i<parametrized.vert.size();i++)
626     {
627         positions[i]=parametrized.vert[i].P();
628         parametrized.vert[i].P()=parametrized.vert[i].RPos;
629     }
630 
631     UpdateTopologies(&parametrized);
632 
633 
634     if (init_border)
635         ParametrizeExternal(parametrized);
636 
637     ParametrizeInternal(parametrized);
638 
639     //for (int i=0;i<parametrized.vert.size();i++)
640     //	parametrized.vert[i].RPos=parametrized.vert[i].P();
641 
642     vcg::tri::MeanValueTexCoordOptimization<MeshType> opt(parametrized);
643     vcg::tri::AreaPreservingTexCoordOptimization<MeshType> opt1(parametrized);
644     opt.SetSpeed((ScalarType)0.0005);
645     InitDampRestUV(parametrized);
646     if (fix_boundary)
647     {
648         opt.TargetEquilateralGeometry();
649         //opt.TargetCurrentGeometry();
650         opt.SetBorderAsFixed();
651         opt.IterateUntilConvergence((ScalarType)0.000001,100);
652         /*opt.Iterate();
653         opt.Iterate();*/
654     }
655     else
656     {
657         opt1.TargetCurrentGeometry();
658         //opt1.SetBorderAsFixed();
659         opt1.IterateUntilConvergence((ScalarType)0.000001,200);
660     }
661 
662     ///assert parametrization
663     for (unsigned int i=0;i<parametrized.face.size();i++)
664     {
665         FaceType *f=&parametrized.face[i];
666         vcg::Point2<ScalarType> tex0=vcg::Point2<ScalarType>(f->V(0)->T().U(),f->V(0)->T().V());
667         vcg::Point2<ScalarType> tex1=vcg::Point2<ScalarType>(f->V(1)->T().U(),f->V(1)->T().V());
668         vcg::Point2<ScalarType> tex2=vcg::Point2<ScalarType>(f->V(2)->T().U(),f->V(2)->T().V());
669         vcg::Triangle2<typename MeshType::ScalarType> t2d=vcg::Triangle2<typename MeshType::ScalarType>(tex0,tex1,tex2);
670 #ifndef NDEBUG
671         ScalarType area=(tex1-tex0)^(tex2-tex0);
672         assert(area>0);
673 #endif
674     }
675     ///restore position
676     for (unsigned int i=0;i<parametrized.vert.size();i++)
677         parametrized.vert[i].P()=positions[i];
678 }
679 #endif
680 template <class MeshType>
ForceInParam(vcg::Point2<typename MeshType::ScalarType> & UV,MeshType & domain)681 void ForceInParam(vcg::Point2<typename MeshType::ScalarType> &UV,MeshType &domain)
682 {
683     typedef typename MeshType::FaceType FaceType;
684     typedef typename MeshType::ScalarType ScalarType;
685 
686     ScalarType minDist=(ScalarType)1000.0;
687     vcg::Point2<ScalarType> closest;
688     vcg::Point2<ScalarType> center=vcg::Point2<ScalarType>(0,0);
689     for (unsigned int i=0;i<domain.face.size();i++)
690     {
691         FaceType *f=&domain.face[i];
692         vcg::Point2<ScalarType> tex0=vcg::Point2<ScalarType>(f->V(0)->T().U(),f->V(0)->T().V());
693         vcg::Point2<ScalarType> tex1=vcg::Point2<ScalarType>(f->V(1)->T().U(),f->V(1)->T().V());
694         vcg::Point2<ScalarType> tex2=vcg::Point2<ScalarType>(f->V(2)->T().U(),f->V(2)->T().V());
695         center+=tex0;
696         center+=tex1;
697         center+=tex2;
698         vcg::Triangle2<ScalarType> t2d=vcg::Triangle2<ScalarType>(tex0,tex1,tex2);
699         ScalarType dist;
700         vcg::Point2<ScalarType> temp;
701         t2d.PointDistance(UV,dist,temp);
702         if (dist<minDist)
703         {
704             minDist=dist;
705             closest=temp;
706         }
707     }
708     center/=(ScalarType)(domain.face.size()*3);
709     UV=closest*(ScalarType)0.95+center*(ScalarType)0.05;
710 }
711 
712 
713 template <class VertexType>
testParamCoords(VertexType * v)714 bool testParamCoords(VertexType *v)
715 {
716     typedef typename VertexType::ScalarType ScalarType;
717     ScalarType eps=(ScalarType)0.00001;
718     if (!(((v->T().P().X()>=-1-eps)&&(v->T().P().X()<=1+eps)&&
719            (v->T().P().Y()>=-1-eps)&&(v->T().P().Y()<=1+eps))))
720             return (false);
721 
722     return true;
723 }
724 
725 template <class MeshType>
testParamCoords(MeshType & domain)726 bool testParamCoords(MeshType &domain)
727 {
728     for (unsigned int i=0;i<domain.vert.size();i++)
729     {
730                 typename MeshType::VertexType *v=&domain.vert[i];
731                 bool b=testParamCoords<typename MeshType::VertexType>(v);
732         if (!b)
733         {
734             #ifndef _MESHLAB
735             printf("\n position %lf,%lf \n",v->T().U(),v->T().V());
736             #endif
737             return false;
738         }
739     }
740     return true;
741 }
742 
743 template <class CoordType>
testBaryCoords(CoordType & bary)744 bool testBaryCoords(CoordType &bary)
745 {
746     typedef typename CoordType::ScalarType ScalarType;
747     ///test
748     float eps=(ScalarType)0.0001;
749     if(!(fabs(bary.X()+bary.Y()+bary.Z()-1.0)<eps))
750         return false;
751     if(!((bary.X()<=1.0)&&(bary.X()>=-eps)&&(bary.Y()<=1.0)&&(bary.Y()>=-eps)&&(bary.Z()<=1.0)&&(bary.Z()>=-eps)))
752         return false;
753     return true;
754 }
755 
756 template <class CoordType>
NormalizeBaryCoords(CoordType & bary)757 bool NormalizeBaryCoords(CoordType &bary)
758 {
759   typedef typename CoordType::ScalarType ScalarType;
760 
761   ScalarType EPS=(ScalarType)0.00000001;
762     bool isOK=testBaryCoords(bary);
763     if (!isOK)
764         return false;
765 
766     typedef typename CoordType::ScalarType ScalarType;
767 
768     ///test <0
769   if (bary.X()<0) bary.X()=EPS;
770   if (bary.Y()<0) bary.Y()=EPS;
771   if (bary.Z()<0) bary.Z()=EPS;
772 
773     ///test >1
774   if (bary.X()>1.0) bary.X()=1.0-EPS;
775   if (bary.Y()>1.0) bary.Y()=1.0-EPS;
776   if (bary.Z()>1.0) bary.Z()=1.0-EPS;
777 
778     ///test sum
779     ScalarType diff=bary.X()+bary.Y()+bary.Z()-1.0;
780     bary.X()-=(diff+EPS);
781 
782     if (bary.X()<0)
783         bary.X()=EPS;
784     return true;
785 }
786 
787 template <class MeshType>
AssingFather(typename MeshType::VertexType & v,typename MeshType::FaceType * father,typename MeshType::CoordType & bary,MeshType & domain)788 void AssingFather(typename MeshType::VertexType &v,
789                                     typename MeshType::FaceType *father,
790                                     typename MeshType::CoordType &bary,
791                   MeshType & domain)
792 {
793 #ifdef _DEBUG
794     const typename MeshType::ScalarType eps=(typename MeshType::ScalarType)0.00001;
795     assert(vcg::tri::IsValidPointer(domain,father));
796     assert(!(father->IsD()));
797     assert(!(father==NULL));
798     assert((bary.X()>=0)&&(bary.X()<=1)&&
799                  (bary.Y()>=0)&&(bary.Y()<=1)&&
800                  (bary.Z()>=0)&&(bary.Z()<=1)&&
801                  ((bary.X()+bary.Y()+bary.Z())<=1+eps)&&
802                  ((bary.X()+bary.Y()+bary.Z())>=1-eps));
803 #endif
804     v.father=father;
805     v.Bary=bary;
806 }
807 
808 template <class MeshType>
testParametrization(MeshType & domain,MeshType & Hlev)809 bool testParametrization(MeshType &domain,
810              MeshType &Hlev)
811 {
812     typedef typename MeshType::FaceType FaceType;
813     typedef typename MeshType::VertexType VertexType;
814     bool is_good=true;
815     int num_del=0;
816     int num_null=0;
817     int fath_son=0;
818     int wrong_address=0;
819     for (unsigned int i=0;i<Hlev.vert.size();i++)
820     {
821         VertexType *v=&Hlev.vert[i];
822         bool isGoodAddr=true;
823         if ((v->father-&(*domain.face.begin()))>=(int)domain.face.size())
824         {
825 //			printf("\n ADDRESS EXCEEDS OF %d \n",v->father-&(*domain.face.begin()));
826             wrong_address++;
827             is_good=false;
828             isGoodAddr=false;
829         }
830         if ((isGoodAddr)&&(v->father==NULL))
831         {
832             //printf("\n PAR ERROR : father NULL\n");
833             num_null++;
834             is_good=false;
835         }
836         if ((isGoodAddr)&&(v->father->IsD()))
837         {
838             //printf("\n PAR ERROR : father DELETED \n");
839             num_del++;
840             is_good=false;
841         }
842         if ((isGoodAddr)&&(!(((v->Bary.X()>=0)&&(v->Bary.X()<=1))&&
843         ((v->Bary.Y()>=0)&&(v->Bary.Y()<=1))&&
844         ((v->Bary.Z()>=0)&&(v->Bary.Z()<=1)))))
845         {
846             printf("\n PAR ERROR 0: bary coords exceeds: %f,%f,%f \n",v->Bary.X(),v->Bary.Y(),v->Bary.Z());
847             /*system("pause");*/
848             NormalizeBaryCoords(v->Bary);
849             is_good=false;
850         }
851     }
852     for (unsigned int i=0;i<domain.face.size();i++)
853     {
854         FaceType *face=&domain.face[i];
855         if (!face->IsD())
856         {
857             for (unsigned int j=0;j<face->vertices_bary.size();j++)
858             {
859                 VertexType *v=face->vertices_bary[j].first;
860                 if (v->father!=face)
861                 {
862                     //printf("\n PAR ERROR : Father<->son \n");
863                     fath_son++;
864                     v->father=face;
865                     is_good=false;
866                 }
867             }
868         }
869     }
870 
871     if (num_del>0)
872         printf("\n PAR ERROR %d Father isDel  \n",num_del);
873     if (num_null>0)
874         printf("\n PAR ERROR %d Father isNull \n",num_null);
875     if (fath_son>0)
876         printf("\n PAR ERROR %d Father<->son  \n",fath_son);
877     if (wrong_address>0)
878     {
879         printf("\n PAR ERROR %d Wrong Address Num Faces %d\n",wrong_address,domain.fn);
880         /*system("pause");*/
881     }
882     return (is_good);
883 }
884 
885 template <class MeshType>
NonFolded(MeshType & parametrized)886 bool NonFolded(MeshType &parametrized)
887 {
888     //const ScalarType epsilon=0.00001;
889     typedef typename MeshType::FaceType FaceType;
890     typedef typename MeshType::ScalarType ScalarType;
891     ///assert parametrization
892     for (unsigned int i=0;i<parametrized.face.size();i++)
893     {
894         FaceType *f=&parametrized.face[i];
895         if (!((f->V(0)->IsB())&&(f->V(1)->IsB())&&(f->V(2)->IsB())))
896         {
897 
898             vcg::Point2<ScalarType> tex0=vcg::Point2<ScalarType>(f->V(0)->T().U(),f->V(0)->T().V());
899             vcg::Point2<ScalarType> tex1=vcg::Point2<ScalarType>(f->V(1)->T().U(),f->V(1)->T().V());
900             vcg::Point2<ScalarType> tex2=vcg::Point2<ScalarType>(f->V(2)->T().U(),f->V(2)->T().V());
901             ScalarType area=(tex1-tex0)^(tex2-tex0);
902             if (area<=0)
903                 return false;
904         }
905     }
906     return true;
907 }
908 
909 template <class MeshType>
NonFolded(MeshType & parametrized,std::vector<typename MeshType::FaceType * > & folded)910 bool NonFolded(MeshType &parametrized,std::vector<typename MeshType::FaceType*> &folded)
911 {
912 
913     typedef typename MeshType::FaceType FaceType;
914     typedef typename MeshType::ScalarType ScalarType;
915 
916     const ScalarType epsilon=(ScalarType)0.00001;
917     folded.resize(0);
918     ///assert parametrization
919     for (unsigned int i=0;i<parametrized.face.size();i++)
920     {
921         FaceType *f=&parametrized.face[i];
922         if (!((f->V(0)->IsB())&&(f->V(1)->IsB())&&(f->V(2)->IsB())))
923         {
924 
925             vcg::Point2<ScalarType> tex0=vcg::Point2<ScalarType>(f->V(0)->T().U(),f->V(0)->T().V());
926             vcg::Point2<ScalarType> tex1=vcg::Point2<ScalarType>(f->V(1)->T().U(),f->V(1)->T().V());
927             vcg::Point2<ScalarType> tex2=vcg::Point2<ScalarType>(f->V(2)->T().U(),f->V(2)->T().V());
928             ScalarType area=(tex1-tex0)^(tex2-tex0);
929             if (area<=epsilon)
930                 folded.push_back(f);
931         }
932     }
933     return (folded.size()==0);
934 }
935 //getFoldedFaces(std::vector)
936 ///parametrize a submesh from trinagles that are incident on vertices with equi-area subdivision
937 template <class MeshType>
938 void ParametrizeStarEquilateral(MeshType &parametrized,
939                                 const typename MeshType::ScalarType &radius=1)
940 {
941     typedef typename MeshType::ScalarType ScalarType;
942     typedef typename MeshType::VertexType VertexType;
943 
944     UpdateTopologies(&parametrized);
945 
946     //set borders
947 
948     ///find first border & non border vertex
949     std::vector<VertexType*> non_border;
950     VertexType* Start=NULL;
951     for (unsigned int i=0;i<parametrized.vert.size();i++)
952     {
953         VertexType* vert=&parametrized.vert[i];
954         if ((Start==NULL)&&(vert->IsB()))
955             Start=vert;
956 
957         if (!vert->IsB())
958             non_border.push_back(vert);
959     }
960     assert(non_border.size()!=0);
961 
962     ///get sorted border vertices
963   std::vector<VertexType*> vertices;
964     FindSortedBorderVertices<MeshType>(parametrized,Start,vertices);
965 
966     ///set border vertices
967     int num=vertices.size();
968         typename std::vector<VertexType*>::iterator iteV;
969     ScalarType curr_angle=0;
970     vertices[0]->T().U()=cos(curr_angle)*radius;
971     vertices[0]->T().V()=sin(curr_angle)*radius;
972     ScalarType division=(2*M_PI)/(ScalarType)num;
973     ///set border
974     for (unsigned int i=1;i<vertices.size();i++)
975     {
976         curr_angle+=division;
977         vertices[i]->T().U()=radius*cos(curr_angle);
978         vertices[i]->T().V()=radius*sin(curr_angle);
979     }
980 
981     if (non_border.size()==1)
982     {
983         ///if non-border vertex is one then set it to zero otherwise
984         ///set it to the average of neighbors
985         non_border[0]->T().P()=vcg::Point2<ScalarType>(0,0);
986     }
987     else
988     {
989         ///set media of star vertices
990         assert(non_border.size()==2);
991         for (unsigned int i=0;i<non_border.size();i++)
992         {
993             VertexType *v=non_border[i];
994             v->T().P()=vcg::Point2<ScalarType>(0,0);
995             int ariety_vert=0;
996             std::vector<VertexType*> star;
997             getVertexStar<MeshType>(v,star);
998             for (unsigned int k=0;k<star.size();k++)
999             {
1000                 if ((!star[k]->IsD())&&(star[k]->IsB()))
1001                 {
1002                     v->T().P()+=star[k]->T().P();
1003                     ariety_vert++;
1004                 }
1005             }
1006             v->T().P()/=(ScalarType)ariety_vert;
1007         }
1008         ///test particular cases
1009         if (!NonFolded(parametrized))
1010         {
1011             std::vector<VertexType*> shared;
1012             getSharedVertexStar<MeshType>(non_border[0],non_border[1],shared);
1013 
1014             assert(shared.size()==2);
1015             assert(shared[0]->IsB());
1016             assert(shared[1]->IsB());
1017             assert(shared[0]!=shared[1]);
1018 
1019             //ScalarType epsilon=(ScalarType)0.001;
1020             ///then get the media of two shared vertices
1021             vcg::Point2<ScalarType> uvAve=shared[0]->T().P()+shared[1]->T().P();
1022             assert(uvAve.Norm()>(ScalarType)0.001);
1023             uvAve.Normalize();
1024             vcg::Point2<ScalarType> p0=uvAve*(ScalarType)0.3;
1025             vcg::Point2<ScalarType> p1=uvAve*(ScalarType)(-0.3);
1026             ///then test and set right assignment
1027             non_border[0]->T().P()=p0;
1028             non_border[1]->T().P()=p1;
1029             if (!NonFolded(parametrized)){
1030                 non_border[0]->T().P()=p1;
1031                 non_border[1]->T().P()=p0;
1032             }
1033         }
1034     }
1035 
1036     ///final assert parametrization
1037     assert(NonFolded(parametrized));
1038 
1039 }
1040 
1041 ///given the mesh and the two edges (same) seen from face[0] and face[1] of the mesh construct
1042 ///a diamond parametrization using equilateral triangles of edge edge_len
1043 template <class MeshType>
1044 void ParametrizeDiamondEquilateral(MeshType &parametrized,
1045                                    const int &edge0,const int &edge1,
1046                                    const typename MeshType::ScalarType &edge_len=1)
1047 {
1048 
1049     typedef typename MeshType::FaceType FaceType;
1050     typedef typename FaceType::ScalarType ScalarType;
1051     typedef typename FaceType::VertexType VertexType;
1052 
1053     ScalarType h=(sqrt(3.0)/2.0)*edge_len;
1054 
1055     FaceType *fd0=&parametrized.face[0];
1056 #ifndef NDEBUG
1057     FaceType *fd1=&parametrized.face[1];
1058 #endif
1059     assert(fd0->FFp(edge0)==fd1);
1060     assert(fd1->FFp(edge1)==fd0);
1061 
1062     ///get 2 vertex on the edge
1063     VertexType *v0=fd0->V(edge0);
1064     VertexType *v1=fd0->V((edge0+1)%3);
1065 
1066 #ifndef NDEBUG
1067     VertexType *vtest0=fd1->V(edge1);
1068     VertexType *vtest1=fd1->V((edge1+1)%3);
1069 
1070     assert(v0!=v1);
1071     assert(vtest0!=vtest1);
1072     assert(((v0==vtest0)&&(v1==vtest1))||((v1==vtest0)&&(v0==vtest1)));
1073 #endif
1074 
1075     ///other 2 vertex
1076     VertexType *v2=parametrized.face[0].V((edge0+2)%3);
1077     VertexType *v3=parametrized.face[1].V((edge1+2)%3);
1078     assert((v2!=v3)&&(v0!=v2)&&(v0!=v3)&&(v1!=v2)&&(v1!=v3));
1079 
1080     ///assign texcoords
1081     v0->T().P()=vcg::Point2<ScalarType>(0,-edge_len/2.0);
1082     v1->T().P()=vcg::Point2<ScalarType>(0,edge_len/2.0);
1083     v2->T().P()=vcg::Point2<ScalarType>(-h,0);
1084     v3->T().P()=vcg::Point2<ScalarType>(h,0);
1085 
1086     ///test
1087     assert(NonFolded(parametrized));
1088 }
1089 
1090 ///given the mesh and the two edges (same) seen from face[0] and face[1] of the mesh construct
1091 ///a diamond parametrization using equilateral triangles of edge edge_len
1092 template <class MeshType>
1093 void ParametrizeFaceEquilateral(MeshType &parametrized,
1094                                    const typename MeshType::ScalarType &edge_len=1)
1095 {
1096     typedef typename MeshType::FaceType FaceType;
1097     typedef typename FaceType::ScalarType ScalarType;
1098 
1099     ScalarType h=(sqrt(3.0)/2.0)*edge_len;
1100 
1101     FaceType *f_param=&(parametrized.face[0]);
1102 
1103     f_param->V(0)->T().P()=vcg::Point2<ScalarType>(edge_len/2.0,0);
1104     f_param->V(1)->T().P()=vcg::Point2<ScalarType>(0,h);
1105     f_param->V(2)->T().P()=vcg::Point2<ScalarType>(-edge_len/2.0,0);
1106 }
1107 
1108 
1109 ///parametrize and create a submesh from trinagles that are incident on
1110 /// vertices .... seturn a vetor of original faces
1111 template <class MeshType>
ParametrizeLocally(MeshType & parametrized,const std::vector<typename MeshType::VertexType * > & subset,std::vector<typename MeshType::FaceType * > & orderedFaces,std::vector<typename MeshType::VertexType * > & orderedVertex)1112 void ParametrizeLocally(MeshType &parametrized,
1113                         const std::vector<typename MeshType::VertexType*> &subset,
1114                         std::vector<typename MeshType::FaceType*> &orderedFaces,
1115                         std::vector<typename MeshType::VertexType*> &orderedVertex)
1116 {
1117     typedef typename MeshType::FaceType FaceType;
1118     typedef typename FaceType::VertexType VertexType;
1119 
1120     orderedFaces.clear();
1121         std::vector<VertexType*> vertices;
1122 
1123 
1124     ///get faces referenced by vertices
1125     getSharedFace<MeshType>(subset,orderedFaces);
1126 
1127     ///do a first copy of the mesh
1128     ///and parametrize it
1129     ///NB: order of faces is maintained
1130     CopyMeshFromFaces<MeshType>(orderedFaces,orderedVertex,parametrized);
1131 
1132     //CreateMeshVertexStar(subset,orderedFaces,parametrized);
1133     ParametrizeLocally(parametrized);
1134 
1135 }
1136 
1137 
1138 
1139 template <class MeshType>
InterpolateUV(const typename MeshType::FaceType * f,const typename MeshType::CoordType & bary,typename MeshType::ScalarType & U,typename MeshType::ScalarType & V)1140 void InterpolateUV(const typename MeshType::FaceType* f,
1141            const typename MeshType::CoordType &bary,
1142            typename MeshType::ScalarType &U,
1143            typename MeshType::ScalarType &V)
1144 {
1145     U=bary.X()*f->cV(0)->T().U()+bary.Y()*f->cV(1)->T().U()+bary.Z()*f->cV(2)->T().U();
1146     V=bary.X()*f->cV(0)->T().V()+bary.Y()*f->cV(1)->T().V()+bary.Z()*f->cV(2)->T().V();
1147     /*if ((!((U>=-1)&&(U<=1)))||(!((V>=-1)&&(V<=1))))
1148     {
1149         printf("Bary:%f,%f,%f \n",bary.X(),bary.Y(),bary.Z());
1150         printf("texCoord:%f,%f \n",U,V);
1151         assert(0);
1152     }*/
1153     //assert ((U>=-1)&&(U<=1));
1154     //assert ((V>=-1)&&(V<=1));
1155 }
1156 
1157 template <class MeshType>
GetBaryFaceFromUV(const MeshType & m,const typename MeshType::ScalarType & U,const typename MeshType::ScalarType & V,typename MeshType::CoordType & bary,int & index)1158 bool GetBaryFaceFromUV(const MeshType &m,
1159                        const typename MeshType::ScalarType &U,
1160                        const typename MeshType::ScalarType &V,
1161                        typename MeshType::CoordType &bary,
1162                        int &index)
1163 {
1164     typedef typename MeshType::ScalarType ScalarType;
1165     const ScalarType _EPSILON = ScalarType(0.0000001);
1166     /*assert ((U>=-1)&&(U<=1));
1167     assert ((V>=-1)&&(V<=1));*/
1168     for (unsigned int i=0;i<m.face.size();i++)
1169     {
1170                 const typename  MeshType::FaceType *f=&m.face[i];
1171         vcg::Point2<ScalarType> tex0=vcg::Point2<ScalarType>(f->cV(0)->T().U(),f->cV(0)->T().V());
1172         vcg::Point2<ScalarType> tex1=vcg::Point2<ScalarType>(f->cV(1)->T().U(),f->cV(1)->T().V());
1173         vcg::Point2<ScalarType> tex2=vcg::Point2<ScalarType>(f->cV(2)->T().U(),f->cV(2)->T().V());
1174 
1175         vcg::Point2<ScalarType> test=vcg::Point2<ScalarType>(U,V);
1176         vcg::Triangle2<ScalarType> t2d=vcg::Triangle2<ScalarType>(tex0,tex1,tex2);
1177         ScalarType area=(tex1-tex0)^(tex2-tex0);
1178         //assert(area>-_EPSILON);
1179         ///then find if the point 2d falls inside
1180         if ((area>_EPSILON)&&(t2d.InterpolationParameters(test,bary.X(),bary.Y(),bary.Z())))
1181         {
1182             index=i;
1183             ///approximation errors
1184             ScalarType sum=0;
1185             for (int x=0;x<3;x++)
1186             {
1187                 if (((bary[x])<=0)&&(bary[x]>=-_EPSILON))
1188                     bary[x]=0;
1189                 else
1190                 if (((bary[x])>=1)&&(bary[x]<=1+_EPSILON))
1191                     bary[x]=1;
1192                 sum+=bary[x];
1193             }
1194             if (sum==0)
1195                 printf("error SUM %f \n",sum);
1196 
1197             bary/=sum;
1198             return true;
1199         }
1200     }
1201 
1202     return (false);
1203 }
1204 
1205 template <class FaceType>
GetBaryFaceFromUV(std::vector<FaceType * > faces,const typename FaceType::ScalarType & U,const typename FaceType::ScalarType & V,typename FaceType::CoordType & bary,int & index)1206 bool GetBaryFaceFromUV(std::vector<FaceType*> faces,
1207                        const typename FaceType::ScalarType &U,
1208                        const typename FaceType::ScalarType &V,
1209                        typename FaceType::CoordType &bary,
1210                        int &index)
1211 {
1212     typedef typename FaceType::ScalarType ScalarType;
1213     typedef typename FaceType::ScalarType ScalarType;
1214     const ScalarType _EPSILON = ScalarType(0.0000001);
1215     /*assert ((U>=-1)&&(U<=1));
1216     assert ((V>=-1)&&(V<=1));*/
1217     for (unsigned int i=0;i<faces.size();i++)
1218     {
1219         FaceType *f=faces[i];
1220         vcg::Point2<ScalarType> tex0=vcg::Point2<ScalarType>(f->V(0)->T().U(),f->V(0)->T().V());
1221         vcg::Point2<ScalarType> tex1=vcg::Point2<ScalarType>(f->V(1)->T().U(),f->V(1)->T().V());
1222         vcg::Point2<ScalarType> tex2=vcg::Point2<ScalarType>(f->V(2)->T().U(),f->V(2)->T().V());
1223 
1224         vcg::Point2<ScalarType> test=vcg::Point2<ScalarType>(U,V);
1225         vcg::Triangle2<ScalarType> t2d=vcg::Triangle2<ScalarType>(tex0,tex1,tex2);
1226         ScalarType area=fabs((tex1-tex0)^(tex2-tex0));
1227         //assert(area>-_EPSILON);
1228         ///then find if the point 2d falls inside
1229         if ((area>_EPSILON)&&(t2d.InterpolationParameters(test,bary.X(),bary.Y(),bary.Z())))
1230         {
1231             index=i;
1232 
1233             ///approximation errors
1234             ScalarType sum=0;
1235             for (int x=0;x<3;x++)
1236             {
1237                 if (((bary[x])<=0)&&(bary[x]>=-_EPSILON))
1238                     bary[x]=0;
1239                 else
1240                 if (((bary[x])>=1)&&(bary[x]<=1+_EPSILON))
1241                     bary[x]=1;
1242                 sum+=fabs(bary[x]);
1243             }
1244             if (sum==0)
1245                 printf("error SUM %f \n",sum);
1246 
1247             bary/=sum;
1248             /*if (!((bary.X()>=0)&& (bary.X()<=1)))
1249                 printf("error %f \n",bary.X());*/
1250             /*ScalarType diff=(1.0-bary.X()-bary.Y()-bary.Z());
1251             bary.X()+=diff;*/
1252             return true;
1253         }
1254     }
1255 
1256     return (false);
1257 }
1258 
1259 template <class MeshType>
GetBaryFaceFromUV(const MeshType & m,const typename MeshType::ScalarType & U,const typename MeshType::ScalarType & V,const std::vector<typename MeshType::FaceType * > & orderedFaces,typename MeshType::CoordType & bary,typename MeshType::FaceType * & chosen)1260 bool GetBaryFaceFromUV(const MeshType &m,
1261                        const typename MeshType::ScalarType &U,
1262                        const typename MeshType::ScalarType &V,
1263                        const std::vector<typename MeshType::FaceType*> &orderedFaces,
1264                        typename MeshType::CoordType &bary,
1265                        typename MeshType::FaceType* &chosen)
1266 {
1267     int index;
1268     bool found=GetBaryFaceFromUV(m,U,V,bary,index);
1269     if(!found) {
1270             chosen=0;
1271             return false;
1272         }
1273     chosen=orderedFaces[index];
1274     return true;
1275 }
1276 
1277 template <class MeshType>
1278 bool GetCoordFromUV(const MeshType &m,
1279                     const typename MeshType::ScalarType &U,
1280                     const typename MeshType::ScalarType &V,
1281                     typename MeshType::CoordType &val,
1282                     bool rpos=false)
1283 {
1284     typedef typename MeshType::ScalarType ScalarType;
1285     const ScalarType _EPSILON = (ScalarType)0.00001;
1286     for (unsigned int i=0;i<m.face.size();i++)
1287     {
1288                 const typename MeshType::FaceType *f=&m.face[i];
1289         vcg::Point2<ScalarType> tex0=vcg::Point2<ScalarType>(f->cV(0)->T().U(),f->cV(0)->T().V());
1290         vcg::Point2<ScalarType> tex1=vcg::Point2<ScalarType>(f->cV(1)->T().U(),f->cV(1)->T().V());
1291         vcg::Point2<ScalarType> tex2=vcg::Point2<ScalarType>(f->cV(2)->T().U(),f->cV(2)->T().V());
1292 
1293         vcg::Point2<ScalarType> test=vcg::Point2<ScalarType>(U,V);
1294         vcg::Triangle2<ScalarType> t2d=vcg::Triangle2<ScalarType>(tex0,tex1,tex2);
1295         ScalarType area=(tex1-tex0)^(tex2-tex0);
1296         ///then find if the point 2d falls inside
1297                 typename MeshType::CoordType bary;
1298         if ((area>_EPSILON)&&(t2d.InterpolationParameters(test,bary.X(),bary.Y(),bary.Z())))
1299         {
1300             ///approximation errors
1301             for (int x=0;x<3;x++)
1302             {
1303                 if (((bary[x])<=0)&&(bary[x]>=-_EPSILON))
1304                     bary[x]=0;
1305                 else
1306                 if (((bary[x])>=1)&&(bary[x]<=1+_EPSILON))
1307                     bary[x]=1;
1308             }
1309             ScalarType diff= (ScalarType)(1.0-bary.X()-bary.Y()-bary.Z());
1310             bary.X()+=diff;
1311             if (!rpos)
1312                 val=f->cP(0)*bary.X()+f->cP(1)*bary.Y()+f->cP(0)*bary.Z();
1313             else
1314                 val=f->cV(0)->RPos*bary.X()+f->cV(1)->RPos*bary.Y()+f->cV(2)->RPos*bary.Z();
1315 
1316             return true;
1317         }
1318     }
1319     return false;
1320 }
1321 
1322 template <class MeshType>
GetSmallestUVEdgeSize(const MeshType & m)1323 typename MeshType::ScalarType GetSmallestUVEdgeSize(const MeshType &m)
1324 {
1325     typedef typename MeshType::ScalarType ScalarType;
1326 
1327     ScalarType smallest=100.f;
1328     assert(m.fn>0);
1329     for (int i=0;i<m.face.size();i++)
1330     {
1331 
1332         ///approximation errors
1333             for (int j=0;j<3;j++)
1334             {
1335 
1336                 vcg::Point2<ScalarType> uv0=m.face[i].V(j)->T().P();
1337                 vcg::Point2<ScalarType> uv1=m.face[i].V((j+1)%3)->T().P();
1338                 ScalarType test=(uv0-uv1).Norm();
1339                 if (test<smallest)
1340                     smallest=test;
1341             }
1342     }
1343     return smallest;
1344 }
1345 
1346 template <class MeshType>
GetSmallestUVHeight(const MeshType & m)1347 typename MeshType::ScalarType GetSmallestUVHeight(const MeshType &m)
1348 {
1349     typedef typename MeshType::ScalarType ScalarType;
1350     ScalarType smallest=(ScalarType)100.0;
1351     ScalarType eps=(ScalarType)0.0001;
1352     assert(m.fn>0);
1353     for (unsigned int i=0;i<m.face.size();i++)
1354     {
1355                 const typename MeshType::FaceType *f=&m.face[i];
1356         ///approximation errors
1357         for (int j=0;j<3;j++)
1358         {
1359             vcg::Point2<ScalarType> uv0=f->cV(j)->cT().P();
1360             vcg::Point2<ScalarType> uv1=f->cV1(j)->cT().P();
1361             vcg::Point2<ScalarType> uv2=f->cV2(j)->cT().P();
1362             ScalarType area=fabs((uv1-uv0)^(uv2-uv0));
1363             ScalarType base=(uv1-uv2).Norm();
1364             ScalarType h_test=area/base;
1365             if (h_test<smallest)
1366                     smallest=h_test;
1367             }
1368     }
1369     if (smallest<eps)
1370         smallest=(ScalarType)eps;
1371     if (smallest>(ScalarType)0.05)
1372         smallest=(ScalarType)0.05;
1373     return smallest;
1374 }
1375 
1376 template <class MeshType>
ParametrizeStarEquilateral(typename MeshType::VertexType * center,bool)1377 void ParametrizeStarEquilateral(typename MeshType::VertexType *center,
1378                                 bool /*subvertices=true*/)
1379 {
1380     ///initialize domain
1381     typedef typename MeshType::VertexType VertexType;
1382     typedef typename MeshType::FaceType FaceType;
1383     typedef typename MeshType::CoordType CoordType;
1384 
1385     MeshType parametrized;
1386     std::vector<VertexType*> vertices,ordVert;
1387     std::vector<VertexType*> HresVert;
1388     std::vector<FaceType*> faces;
1389     vertices.push_back(center);
1390     getSharedFace<MeshType>(vertices,faces);
1391     CopyMeshFromFaces<MeshType>(faces,ordVert,parametrized);
1392 
1393     ///parametrize and then copy back
1394     ParametrizeStarEquilateral<MeshType>(parametrized);
1395     for (unsigned int i=0;i<ordVert.size();i++)
1396         ordVert[i]->T().P()=parametrized.vert[i].T().P();
1397 
1398     ///initialize sub-vertices
1399     getHresVertex<FaceType>(faces,HresVert);
1400     for (unsigned int i=0;i<HresVert.size();i++)
1401     {
1402         FaceType *father=HresVert[i]->father;
1403         CoordType Bary=HresVert[i]->Bary;
1404     InterpolateUV<MeshType>(father,Bary,HresVert[i]->T().U(),HresVert[i]->T().V());
1405     }
1406 }
1407 
1408 #endif
1409