1 // to be included in the library
2 
3 #ifndef __VCGLIB__TEXTCOOORD_OPTIMIZATION
4 #define __VCGLIB__TEXTCOOORD_OPTIMIZATION
5 
6 
7 #include <vcg/container/simple_temporary_data.h>
8 #ifdef _USE_OMP
9 #include <omp.h>
10 #endif
11 /*
12 
13 SINGLE PATCH TEXTURE OPTIMIZATIONS
14 
15 A set of classes to perform optimizations of disk->disk parametrizations.
16 
17 Requires texture coords to be defined per vertex (do replicate seams!).
18 
19 */
20 
21 
22 namespace vcg
23 {
24 namespace tri
25 {
26 
27 // helper function (checks that coords are inside -1..+1)
28 template <class ScalarType>
testParamCoordsPoint(const vcg::Point2<ScalarType> & p)29 bool testParamCoordsPoint(const vcg::Point2<ScalarType> &p)
30 {
31 	ScalarType eps=(ScalarType)0.00001;
32 	if (!((p.X()>=-1-eps)&&(p.X()<=1+eps)&&
33 		   (p.Y()>=-1-eps)&&(p.Y()<=1+eps)))
34 			return (false);
35 	return true;
36 }
37 
38 /* Base class for all Texture Optimizers*/
39 template<class MESH_TYPE>
40 class TexCoordOptimization{
41 protected:
42   MESH_TYPE &m;
43   SimpleTempData<typename MESH_TYPE::VertContainer, int > isFixed;
44 public:
45 
46   /* Tpyes */
47   typedef MESH_TYPE MeshType;
48   typedef typename MESH_TYPE::VertexIterator VertexIterator;
49   typedef typename MESH_TYPE::FaceIterator FaceIterator;
50   typedef typename MESH_TYPE::VertexType VertexType;
51   typedef typename MESH_TYPE::FaceType FaceType;
52   typedef typename MESH_TYPE::ScalarType ScalarType;
53 
54 
55   /* Access functions */
Mesh()56   const MeshType & Mesh() const {return m;}
Mesh()57   MeshType & Mesh() {return m;}
58 
59   /* Constructior */
TexCoordOptimization(MeshType & _m)60   TexCoordOptimization(MeshType &_m):m(_m),isFixed(_m.vert){
61    /* assert(m.HasPerVertexTexture());*/
62   }
63 
64   // initializes on current geometry
65   virtual void TargetCurrentGeometry()=0;
66 
67   // performs an interation. Returns largest movement.
68   virtual ScalarType Iterate()=0;
69 
70   // performs an iteration (faster, but it does not tell how close it is to stopping)
71   virtual void IterateBlind()=0;
72 
73   // performs <steps> iteration
IterateN(int step)74   virtual ScalarType IterateN(int step){
75     for (int i=0; i<step-1; i++) {
76       this->IterateBlind();
77     }
78     if (step>1) return this->Iterate(); else return 0;
79   }
80 
81   // performs iterations until convergence.
82   virtual int IterateUntilConvergence(ScalarType threshold=0.0001, int maxite=5000){
83     int i=0;
84     while (Iterate()>threshold) {
85 	  if (i++>maxite) return i;
86     }
87     return i;
88   }
89 
90   // desctuctor: free temporary field
~TexCoordOptimization()91   ~TexCoordOptimization(){
92     /*isFixed.Stop();*/
93   };
94 
95   // set the current border as fixed (forced to stay in position during text optimization)
SetBorderAsFixed()96   void SetBorderAsFixed(){
97     /*isFixed.Start();*/
98     for (VertexIterator v=m.vert.begin(); v!=m.vert.end(); v++) {
99 		  isFixed[v]=(v->IsB())?1:0;
100 	  }
101   }
102 
103   // everything moves, no vertex must fixed during texture optimization)
SetNothingAsFixed()104   void SetNothingAsFixed(){
105     //isFixed.Start();
106     for (VertexIterator v=m.vert.begin(); v!=m.vert.end(); v++) {
107 		  isFixed[v]=0;
108 	  }
109   }
110 
111   // fix a given vertex
112   void FixVertex(const VertexType *v, bool fix=true){
113     isFixed[v]=(fix)?1:0;
114   }
115 
IsFixed(const VertexType * v)116   bool IsFixed(const VertexType *v)
117   {
118 	return (isFixed[v]);
119   }
120 
Fixed(const FaceType * f)121   bool Fixed(const FaceType* f)
122   {
123 	return ((isFixed[f->V(0)])&&(isFixed[f->V(1)])&&(isFixed[f->V(2)]));
124   }
125 
SetSpeed(ScalarType)126   virtual void SetSpeed(ScalarType){
127     //assert(0); // by default, no speed to set
128   }
129 
GetSpeed()130   virtual ScalarType GetSpeed(){
131     //assert(0); // by default, no speed to set
132     return 0;
133   }
134 
SetTheta(int)135   virtual void SetTheta(int){
136     assert(0);
137   }
138 
GetTheta()139   virtual int GetTheta(){
140     assert(0);
141     return 0;
142   }
143 
144 };
145 
146 /*
147 AREA PRESERVING TEXTURE OPTIMIZATION
148 
149 as in: Degener, P., Meseth, J., Klein, R.
150        "An adaptable surface parameterization method."
151        Proc. of the 12th International Meshing Roundtable, 201-213 [2003].
152 
153 Features:
154 
155 :) - Balances angle and area distortions (best results!).
156 :) - Can choose how to balance area and angle preservation (see SetTheta)
157        theta=0 -> pure conformal (use MIPS instead!)
158        theta=3 -> good balance between area and angle preservation
159        theta>3 -> care more about area than about angles
160 :( - Slowest method.
161 :( - Requires a fixed boundary, else expands forever in texture space (unless theta=0).
162 :( - Diverges in presence of flipped faces (unless theta=0).
163 :( - Requires a speed parameter to be set.
164        Speed too large => when close, bounces back and forth around minimum, w/o getting any closer.
165        Lower speed => longer convercence times
166 */
167 
168 template<class MESH_TYPE>
169 class AreaPreservingTexCoordOptimization:public TexCoordOptimization<MESH_TYPE>{
170 public:
171   /* Types */
172   typedef MESH_TYPE MeshType;
173   typedef typename MESH_TYPE::VertexIterator VertexIterator;
174   typedef typename MESH_TYPE::FaceIterator FaceIterator;
175   typedef typename MESH_TYPE::VertexType VertexType;
176   typedef typename MESH_TYPE::FaceType FaceType;
177   typedef typename MESH_TYPE::ScalarType ScalarType;
178   typedef typename MESH_TYPE::CoordType CoordType;
179 
180 
181 private:
182   typedef TexCoordOptimization<MESH_TYPE> Super; // superclass (commodity)
183 
184   // extra data per face: [0..3] -> cotangents. [4] -> area*2
185   SimpleTempData<typename MESH_TYPE::FaceContainer, Point4<ScalarType> > data;
186   SimpleTempData<typename MESH_TYPE::VertContainer, Point2<ScalarType> > sum;
187 
188   std::vector<CoordType> sumX;
189   std::vector<CoordType> sumY;
190 
191   SimpleTempData<typename MESH_TYPE::VertContainer, Point2<ScalarType> > lastDir;
192   /*SimpleTempData<typename MESH_TYPE::VertContainer,omp_lock_t> lck;*/
193   SimpleTempData<typename MESH_TYPE::VertContainer, ScalarType > vSpeed;
194 
195   ScalarType totArea;
196   ScalarType speed;
197 
198   int theta;
199 
200 public:
201 
202 
203   // constructor and destructor
AreaPreservingTexCoordOptimization(MeshType & _m)204   AreaPreservingTexCoordOptimization(MeshType &_m):Super(_m),data(_m.face),sum(_m.vert),lastDir(_m.vert),vSpeed(_m.vert,1){
205     speed=(ScalarType)0.00005;
206     theta=3;
207 
208 	/*for (int i=0;i<m.vert.size();i++)
209 		omp_init_lock(&lck[i]);*/
210 
211   }
212 
~AreaPreservingTexCoordOptimization()213   ~AreaPreservingTexCoordOptimization(){
214    /* data.Stop();
215     sum.Stop();
216     Super::isFixed.Stop();*/
217   }
218 
SetSpeed(ScalarType _speed)219   void SetSpeed(ScalarType _speed){
220     speed=_speed;
221   }
222 
GetSpeed()223   ScalarType GetSpeed(){
224     return speed;
225   }
226 
227   // sets the parameter theta:
228   // good parameters are in 1..3
229   //  0 = converge to pure conformal, ignore area preservation
230   //  3 = good balance between area and conformal
231   // >3 = area more important, angle preservation less important
SetTheta(int _theta)232   void SetTheta(int _theta){
233     theta=_theta;
234   }
235 
GetTheta()236   int GetTheta(){
237     return theta;
238   }
239 
IterateBlind()240   void IterateBlind(){
241     /* todo: do as iterate, but without */
242     Iterate();
243   }
244 
Area(int i)245   ScalarType Area(int i)
246   {
247          FaceType *f=&(Super::m.face[i]);
248 	 double val=0;
249 	 if (!(Super::isFixed[f->V(0)]&& Super::isFixed[f->V(1)] && Super::isFixed[f->V(2)]))
250 		val=(f->V(1)->T().P()-f->V(0)->T().P())^(f->V(2)->T().P()-f->V(0)->T().P());
251 
252 	/* bool b0=testParamCoords(f->V(0));
253 	 bool b1=testParamCoords(f->V(1));
254 	 bool b2=testParamCoords(f->V(2));*/
255 
256 	 if(!((fabs(val)<3.14)&&(fabs(val)>=0.0)))
257 	 {
258 		 printf("v0 %lf,%lf \n",f->V(0)->T().U(),f->V(0)->T().V());
259 		 printf("v1 %lf,%lf \n",f->V(1)->T().U(),f->V(1)->T().V());
260 		 printf("v2 %lf,%lf \n",f->V(2)->T().U(),f->V(2)->T().V());
261 		 printf("Area Value %lf \n",val);
262 		 //system("pause");
263 	 }
264 
265 	 return fabs(val);
266   }
267 
InitSum()268  void InitSum()
269  {
270 	 int k;
271          auto n=Super::m.vert.size();
272          auto n1=Super::m.face.size();
273 #ifdef _USE_OMP
274 	 #pragma omp parallel for default (none) shared(n) private(k)
275 #endif
276 	 for (k=0;k<n;k++)
277 	 {
278 	   sum[k]=Point2<ScalarType>(0,0);
279 	 }
280 #ifdef _USE_OMP
281 	 #pragma omp parallel for default (none) shared(n1) private(k)
282 #endif
283 	 for (k=0;k<n1;k++)
284 	 {
285 	   sumX[k].X()=0;
286 	   sumX[k].Y()=0;
287 	   sumX[k].Z()=0;
288 	   sumY[k].X()=0;
289 	   sumY[k].Y()=0;
290 	   sumY[k].Z()=0;
291 	 }
292 #ifdef _USE_OMP
293 	 #pragma omp barrier
294 #endif
295  }
296 
getProjArea()297 ScalarType getProjArea()
298 {
299 	  int k;
300           int n=Super::m.face.size();
301 	  ScalarType tot_proj_area=0;
302 	 //# pragma omp parallel for
303 #ifdef _USE_OMP
304 	  #pragma omp parallel for default (none) shared(n) private(k) reduction(+: tot_proj_area)
305 #endif
306 	  for (k=0;k<n; k++) {
307 	      tot_proj_area+=Area(k);
308 	  }
309 #ifdef _USE_OMP
310 	  #pragma omp barrier
311 #endif
312 	  return (tot_proj_area);
313 }
314 
VertValue(const int & face,const int & vert,const double & scale)315 vcg::Point2<ScalarType> VertValue(const int &face,const int &vert,const double &scale)
316 {
317          FaceType *f=&Super::m.face[face];
318 	 /*int i=0;*/
319 	 vcg::Point2<ScalarType> t0=(f->V0(vert)->T().P());
320      vcg::Point2<ScalarType> t1=(f->V1(vert)->T().P());
321      vcg::Point2<ScalarType> t2=(f->V2(vert)->T().P());
322 	 ScalarType area2 = fabs((t1-t0) ^ (t2-t0));
323 	 ScalarType  a = (t1-t0).Norm(),
324 				 b =  ((t1-t0) * (t2-t0))/a,
325 			     c = area2 / a,
326 
327 				 m0= data[face][vert] / area2,
328 				 m1= data[face][(vert+1)%3] / area2,
329 				 m2= data[face][(vert+2)%3] / area2,
330 
331 				 mx= (b-a)/area2,
332 				 my= c/area2, // 1.0/a
333 				 mA= data[face][3]/area2* scale,
334 				 e = m0*((b-a)*(b-a)+c*c) + m1*(b*b+c*c) + m2*a*a, // as obvious
335 				 M1= mA + 1.0/mA,
336 				 M2= mA - 1.0/mA,
337 				 px= e*my,
338 				 py=-e*mx,
339 				 qx= m1*b+ m2*a,
340 				 qy= m1*c,
341 
342 
343 				 dx=pow(M1,theta-1)
344 				  	 *(px*(M1+ theta*M2) - 2.0*qx*M1),
345 
346 				 dy=pow(M1,theta-1)
347 					   *(py*(M1+ theta*M2) - 2.0*qy*M1),
348 
349 				 gy= dy/c,
350 				 gx= (dx - gy*b) / a;
351 
352 				  // 3d gradient
353 				Point2<ScalarType> val=( (t1-t0) * gx + (t2-t0) * gy ) * data[face][3];
354 
355 				return val;
356 }
357 
UpdateSum(const double & scale)358 void UpdateSum(const double &scale)
359 {
360          int n=Super::m.face.size();
361 	 int k;
362 	 FaceType *f;
363 	 ScalarType myscale=scale;
364 	 vcg::Point2<ScalarType> val0,val1,val2;
365 #ifdef _USE_OMP
366 	  #pragma omp parallel for default (none) shared(n,myscale) private(k,f,val0,val1,val2)
367 #endif
368 	  for (k=0;k<n; k++) {
369                           f=&Super::m.face[k];
370 			  val0=VertValue(k,0,myscale);
371 			  val1=VertValue(k,1,myscale);
372 			  val2=VertValue(k,2,myscale);
373 			  sumX[k].V(0)=val0.X();
374 			  sumX[k].V(1)=val1.X();
375 			  sumX[k].V(2)=val2.X();
376 			  sumY[k].V(0)=val0.Y();
377 			  sumY[k].V(1)=val1.Y();
378 			  sumY[k].V(2)=val2.Y();
379 	  }
380 #ifdef _USE_OMP
381 	  #pragma omp barrier
382 #endif
383 }
384 
385 
SumVertex()386 void SumVertex()
387 {
388 	for (unsigned int j=0; j<Super::m.face.size(); j++)
389 	{
390 		for (int i=0;i<3;i++)
391 		{
392                         VertexType *v=Super::m.face[j].V(i);
393 			sum[v].X()+=sumX[j].V(i);
394 			sum[v].Y()+=sumY[j].V(i);
395 		}
396 	}
397 }
398 
Iterate()399  ScalarType Iterate(){
400 
401 	InitSum();
402 
403 	ScalarType tot_proj_area=getProjArea();
404 
405 
406 	double scale= tot_proj_area / totArea ;
407 
408 	UpdateSum(scale);
409 
410 
411     ScalarType max=0; // max displacement
412 	// #pragma omp parallel for num_threads(4)
413  	 // for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++)
414 	/*omp_lock_t lck0;
415 	omp_init_lock(&lck0);*/
416 
417 	SumVertex();
418 
419 	//#pragma omp parallel for
420 	for (unsigned int j=0; j<Super::m.vert.size(); j++)
421 	{
422 		VertexType *v=&Super::m.vert[j];
423 
424     if (  !Super::isFixed[v] ) //if (!v->IsB())
425     {
426 		  ScalarType n=sum[v].Norm();
427 		  //printf("N %f \n",n);
428 		  if ( n > 1 ) { sum[v]/=n; n=1.0;}
429 
430 		  if (lastDir[v]*sum[v]<0) vSpeed[v]*=(ScalarType)0.85; else vSpeed[v]/=(ScalarType)0.92;
431 		  lastDir[v]= sum[v];
432 
433 		 /* if ( n*speed<=0.1 );
434 		  {*/
435 
436 			vcg::Point2f goal=v->T().P()-(sum[v] * (speed * vSpeed[v]) );
437 			bool isOK=testParamCoordsPoint<ScalarType>(goal);
438 			if (isOK)
439 				v->T().P()-=(sum[v] * (speed * vSpeed[v]) );
440 
441 
442 			n=n*speed * vSpeed[v];
443 			//#pragma omp critical
444 			max=std::max(max,n);
445      /* }*/
446 	}
447 	}
448 	return max;
449  }
450 
451 
TargetCurrentGeometry()452   void TargetCurrentGeometry(){
453 
454    /* Super::isFixed.Start();
455     data.Start();
456     sum.Start();*/
457       sumX.resize(Super::m.face.size());
458           sumY.resize(Super::m.face.size());
459 	  totArea=0;
460 	  for (FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++) {
461 		  double area2 = 	((f->V(1)->P() - f->V(0)->P() )^(f->V(2)->P() - f->V(0)->P() )).Norm();
462 		  totArea+=area2;
463 
464 		  //if (  Super::isFixed[f->V1(0)] )
465 		  for (int i=0; i<3; i++){
466 			  data[f][i]=(
467 				(f->V1(i)->P() - f->V0(i)->P() )*(f->V2(i)->P() - f->V0(i)->P() )
468 			  )/area2;
469 			  data[f][3]=area2;
470 		  }
471 	  }
472   }
473 
474 };
475 
476 
477 /*
478 MIPS TEXTURE OPTIMIZATION
479 
480 Features:
481 
482 :( - Targets angle distortions only (not ideal for texture mapping).
483 :) - Quite fast.
484 :) - Does not require fixed boundary (will auto-find a boundary -- up to a scale).
485 :) - Tends to nicely heal flipped faces
486 :( - Requires a speed parameter to be set
487     (SHOULD NOT BE LIKE THIS. BETTER IMPLEMENTATION NEEDED).
488        Speed too large => when close, bounces back and forth around minimum, w/o getting any closer.
489        Lower speed => longer convercence times
490 
491 */
492 
493 template<class MESH_TYPE>
494 class MIPSTexCoordOptimization:public TexCoordOptimization<MESH_TYPE>{
495 public:
496   /* Types */
497   typedef MESH_TYPE MeshType;
498   typedef typename MESH_TYPE::VertexIterator VertexIterator;
499   typedef typename MESH_TYPE::FaceIterator FaceIterator;
500   typedef typename MESH_TYPE::VertexType VertexType;
501   typedef typename MESH_TYPE::FaceType FaceType;
502   typedef typename MESH_TYPE::ScalarType ScalarType;
503 
504 
505 protected:
506   typedef TexCoordOptimization<MESH_TYPE> Super; // superclass (commodity)
507 
508   // extra data per face: [0..3] -> cotangents.
509   SimpleTempData<typename MESH_TYPE::FaceContainer, Point3<ScalarType> > data;
510   SimpleTempData<typename MESH_TYPE::VertContainer, Point2<ScalarType> > sum;
511 
512   ScalarType totArea;
513   ScalarType speed;
514 
515 public:
516 
517 
518   // constructor and destructor
MIPSTexCoordOptimization(MeshType & _m)519   MIPSTexCoordOptimization(MeshType &_m):Super(_m),data(_m.face),sum(_m.vert){
520     speed=(ScalarType)0.001;
521   }
522 
~MIPSTexCoordOptimization()523   ~MIPSTexCoordOptimization(){
524    /* data.Stop();
525     sum.Stop();
526     Super::isFixed.Stop();*/
527   }
528 
SetSpeed(ScalarType _speed)529   void SetSpeed(ScalarType _speed){
530     speed=_speed;
531   }
532 
GetSpeed()533   ScalarType GetSpeed(){
534     return speed;
535   }
536 
IterateBlind()537   void IterateBlind(){
538     /* todo: do as iterate, but without */
539     Iterate();
540   }
541 
Iterate()542   ScalarType Iterate(){
543 
544     #define v0 (f->V(0)->T().P())
545     #define v1 (f->V(1)->T().P())
546     #define v2 (f->V(2)->T().P())
547     #define vi (f->V(i)->T().P())
548     #define vj (f->V(j)->T().P())
549     #define vk (f->V(k)->T().P())
550 	  for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++) {
551 		  //sum[v].Zero();
552 		  sum[v]=Point2<ScalarType>(0,0);
553 	  }
554 
555     ScalarType totProjArea=0;
556 	  for (FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++) {
557       ScalarType area2 = ((v1-v0) ^ (v2-v0));
558       totProjArea+=area2;
559       ScalarType o[3] = { // (opposite edge)^2
560         (   v1-v2).SquaredNorm(),
561         (v0   -v2).SquaredNorm(),
562         (v0-v1   ).SquaredNorm(),
563       };
564       ScalarType e =( data[f][0] * o[0] +
565                       data[f][1] * o[1] +
566                       data[f][2] * o[2] ) / (area2*area2);
567 
568 		  for (int i=0; i<3; i++){
569         int j=(i+1)%3, k=(i+2)%3;
570 			  ScalarType p=(vj-vi)*(vk-vi);
571 				ScalarType gy= e*(o[k]-p) - 2*data[f][j];
572 				ScalarType gx= e*(o[j]-p) - 2*data[f][k];
573 
574 				// speed free mode: (a try!)
575 			  //sum[f->V(i)]+= ( (vj-vi) * gx + (vk-vi) * gy );// / area2;
576 
577 			  // speed mode:
578         sum[f->V(i)]+= ( (vj-vi) * gx + (vk-vi) * gy ) / area2;
579 		  }
580 	  }
581 
582     ScalarType max=0; // max displacement
583 
584  	  for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++)
585     if (  !Super::isFixed[v] )
586     {
587       // speed free mode: (a try!)
588       //v->T().P()-=speed * sum[v] *totProjArea/totArea;
589 
590       // speed mode:
591 
592       ScalarType n=sum[v].Norm(); if ( n > 1 ) { sum[v]/=n; n=1.0;}
593 		  v->T().P()-=(sum[v] ) * speed ;
594 		  if (max<n) max=n;
595 
596   	}
597   	return max;
598   	#undef v0
599     #undef v1
600     #undef v2
601   	#undef vi
602     #undef vj
603     #undef vk
604   	//printf("rejected %d\n",rejected);
605   }
606 
607 
608 
609 
TargetCurrentGeometry()610   void TargetCurrentGeometry(){
611 
612    /* Super::isFixed.Start();
613     data.Start();
614     sum.Start();*/
615 
616 	  totArea=0;
617 	  for (FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++) {
618 		  double area2 = 	((f->V(1)->P() - f->V(0)->P() )^(f->V(2)->P() - f->V(0)->P() )).Norm();
619 		  totArea+=area2;
620 		  for (int i=0; i<3; i++){
621 			  data[f][i]=(
622 				  (f->V1(i)->P() - f->V0(i)->P() )*(f->V2(i)->P() - f->V0(i)->P() )
623 			  );
624         // / area2;
625 		  }
626 	  }
627   }
628 
629 };
630 
631 
632 #if 0  // Temporarily commented out. It still have to be thoroughly tested...
633 
634 template<class MESH_TYPE>
635 class WachspressTexCoordOptimization:public TexCoordOptimization<MESH_TYPE>{
636 public:
637   /* Types */
638   typedef MESH_TYPE MeshType;
639   typedef typename MESH_TYPE::VertexIterator VertexIterator;
640   typedef typename MESH_TYPE::FaceIterator FaceIterator;
641   typedef typename MESH_TYPE::VertexType VertexType;
642   typedef typename MESH_TYPE::FaceType FaceType;
643   typedef typename MESH_TYPE::VertexType::TexCoordType::PointType PointType;
644   typedef typename MESH_TYPE::VertexType::TexCoordType::PointType::ScalarType ScalarType;
645 
646 private:
647   class Factors{
648     public:
649     ScalarType data[3][2];
650   };
651 
652   typedef TexCoordOptimization<MESH_TYPE> Super; // superclass (commodity)
653 
654   // extra data per face: [0..3] -> cotangents. [4] -> area*2
655   SimpleTempData<typename MESH_TYPE::FaceContainer, Factors > factors;
656 
657 public:
658 
659 
660   // constructor and destructor
661   WachspressTexCoordOptimization(MeshType &_m):Super(_m),factors(_m.face){
662   }
663 
664   ~WachspressTexCoordOptimization(){
665    /* factors.Stop();
666     Super::isFixed.Stop();*/
667   }
668 
669   void IterateBlind(){
670     /* todo: do as iterate, but without */
671     Iterate();
672   }
673 
674   ScalarType Iterate(){
675 
676     ScalarType max; // max displacement
677 
678 
679     for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++) {
680 		  v->div=0; v->sum.SetZero();
681 	  }
682 
683     #define vi0 (f->V(i0)->P())
684     #define vi1 (f->V(i1)->P())
685     #define vi2 (f->V(i2)->P())
686     #define EPSILON 1e-4
687 
688 	  for (FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++){
689 		  ScalarType A=(((f->V(0)->P()) - (f->V(1)->P()))^((f->V(0)->P()) - (f->V(2)->P()))).Norm();
690 		  if (A<EPSILON) continue;
691 		  for (int i0=0; i0<3; i0++) {
692 			  int i1=(i0+1)%3,i2=(i0+2)%3;
693 
694 
695   		  ScalarType fact = (vi1-vi0)*(vi2-vi0)/A;
696 				//fact=1;
697 			  if ( (!f->V(i1)->IsB()) /*|| f->V(o)->IsB()*/);{
698 				  f->V(i1)->sum += f->V(i0)->T().P() * fact;
699 				  f->V(i1)->div += fact;
700 			  }
701 			  if ( (!f->V(i2)->IsB()) /*|| f->V(o)->IsB()*/);{
702 				  f->V(i2)->sum += f->V(i0)->T().P() * fact;
703 				  f->V(i2)->div += fact;
704 			  }
705 		  }
706 	  }
707 
708 	  for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++) {
709       if (  !Super::isFixed[v] )
710 		  if (v->div>0.001) {
711 			  v->T().P() = v->sum/v->div;
712 		  }
713 	  }
714     return max;
715   }
716 
717 
718   void TargetCurrentGeometry(){
719   }
720 
721 };
722 
723 #endif
724 
725 template<class MESH_TYPE>
726 class MeanValueTexCoordOptimization:public TexCoordOptimization<MESH_TYPE>{
727 public:
728   /* Types */
729   typedef MESH_TYPE MeshType;
730   typedef typename MESH_TYPE::VertexIterator VertexIterator;
731   typedef typename MESH_TYPE::FaceIterator FaceIterator;
732   typedef typename MESH_TYPE::VertexType VertexType;
733   typedef typename MESH_TYPE::FaceType FaceType;
734   typedef typename MESH_TYPE::VertexType::TexCoordType::PointType PointType;
735   typedef typename MESH_TYPE::VertexType::TexCoordType::PointType::ScalarType ScalarType;
736 
737 private:
738   class Factors{
739     public:
740     ScalarType data[3][2];
741   };
742 
743   typedef TexCoordOptimization<MESH_TYPE> Super; // superclass (commodity)
744 
745   // extra data per face: factors
746   SimpleTempData<typename MESH_TYPE::FaceContainer, Factors > factors;
747 
748   // extra data per vertex: sums and div
749   SimpleTempData<typename MESH_TYPE::VertContainer, PointType > sum;
750   SimpleTempData<typename MESH_TYPE::VertContainer, ScalarType > div;
751 
752 public:
753 
754 
755   // constructor and destructor
MeanValueTexCoordOptimization(MeshType & _m)756   MeanValueTexCoordOptimization(MeshType &_m):Super(_m),factors(_m.face),sum(_m.vert),div(_m.vert){
757   }
758 
~MeanValueTexCoordOptimization()759   ~MeanValueTexCoordOptimization(){
760    /* factors.Stop();
761     sum.Stop();
762     div.Stop();
763     Super::isFixed.Stop();*/
764   }
765 
IterateBlind()766   void IterateBlind(){
767     /* todo: do as iterate, but without */
768     Iterate();
769   }
770 
Iterate()771   ScalarType Iterate(){
772 
773     ScalarType max=0; // max displacement
774 
775 	  for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++) {
776 		  sum[v]=PointType(0,0);
777 		  div[v]=0;
778 	  }
779 
780 	  for (FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++){
781 		  for (int i=0; i<3; i++)
782 		  for (int j=1; j<3; j++) {
783 			  int d=i, o=(i+3-j)%3;
784 			  sum[f->V(d)] += f->V(o)->T().P() * factors[f].data[i][j-1];
785 			  div[f->V(d)] += factors[f].data[i][j-1];
786 		  }
787 	  }
788 
789 	  for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++)
790     if (  !Super::isFixed[v] )
791 	  if (		div[v]>0.000001 ) {
792 		  PointType swap=v->T().P();
793 		  PointType goal=sum[v]/div[v];
794 
795 		  v->T().P() = goal*(ScalarType)0.1+swap*(ScalarType)0.9;
796 
797 		  //v->T().P()=v->RestUV*(1-v->Damp)+(sum[v]/div[v])*(v->Damp);
798 		  ScalarType temp=(swap-v->T().P()).SquaredNorm();
799 		  if (max<temp)
800 			  max=temp;
801 	  }
802     return max;
803   }
804 
TargetEquilateralGeometry()805   void TargetEquilateralGeometry(){
806 	  for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++) {
807 		  div[v]=0;
808 	  }
809 	  const ScalarType fact= 1.0 / sqrt(3.0);
810 	  for (FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++){
811 		  for (int i=0; i<3; i++)
812 			  for (int j=0; j<2; j++) {
813 				  factors[f].data[i][j] = fact;
814 				  div[f->V(i)] += fact ;
815 			  }
816 	  }
817   }
818 
TargetCurrentGeometry()819   void TargetCurrentGeometry(){
820 
821     /*Super::isFixed.Start();
822     factors.Start();
823     sum.Start();
824     div.Start();*/
825 
826     for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++) {
827 		  div[v]=0;
828 	  };
829 	  for (FaceIterator  f=Super::m.face.begin(); f!=Super::m.face.end(); f++){
830 		  for (int i=0; i<3; i++)
831 		  for (int j=1; j<3; j++) factors[f].data[i][j-1]=0;
832 	  };
833 
834     #define vs (f->V(s)->P())
835     #define vd (f->V(d)->P())
836     #define vo (f->V(o)->P())
837     #define EPSILON 1e-4
838 
839 	  for (FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++){
840 		  int s=0,d=1,o=2;
841 		  ScalarType A=((vs - vd)^(vs - vo)).Norm();
842 		  if (A<EPSILON) break;
843 		  for (int i=0; i<3; i++)
844 		  for (int j=1; j<3; j++) {
845 			  d=i; s=(i+j)%3; o=(i+3-j)%3;
846 			  {
847 
848 				  ScalarType dd=((vd-vs).Norm());
849 				  if (dd<=0.0001) continue;
850 				  ScalarType fact= ( ( vd -vo ).Norm() - ((vd-vo)*(vd-vs))/dd) /A;
851 
852 				  //if (fact<0) printf("AAAGH!");
853 				  factors[f].data[d][j-1] = fact;
854 				  //f->V(d)->sum += f->V(o)->projected * fact;
855 				  div[f->V(d)] += fact ;
856 			  } //else {
857 				  //printf(".");
858 			  //}
859 		  }
860 	  }
861 	  for (FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++){
862       for (int i=0; i<3; i++)
863 		  for (int j=1; j<3; j++)
864 		  	if (div[f->V(i)]>0.0001) {
865 		  	  //fact[i][j-1]/=div[f->V(i)];
866 		  	}
867 		  	//else f->fact[i][j-1]=0.0;
868     }
869 
870 	  /*
871     for (f=face.begin(); f!=face.end(); f++)  {
872 	  	for (int i=0; i<3; i++)
873 		  for (int j=1; j<3; j++)
874 		  	if (f->V(i)->div>0.01) {
875 		  	  f->fact[i][j-1]/=f->V(i)->div;
876 		  	}
877 		  	else f->fact[i][j-1]=0.0;
878 	  } */
879 
880 	  #undef vs
881     #undef vd
882     #undef vo
883 
884   }
885 
886 };
887 
888 
889 /* texture coords general utility functions */
890 /*++++++++++++++++++++++++++++++++++++++++++*/
891 
892 // returns false if any fold is present (faster than MarkFolds)
893 template<class MESH_TYPE>
IsTexCoordFoldFree(MESH_TYPE & m)894 bool IsTexCoordFoldFree(MESH_TYPE &m){
895 
896   assert(HasPerVertexTexCoord(m));
897 
898   typedef typename MESH_TYPE::VertexType::TexCoordType::PointType::ScalarType ScalarType;
899 
900   ScalarType lastsign=0;
901   for (typename MESH_TYPE::FaceIterator f=m.face.begin(); f!=m.face.end(); f++){
902     ScalarType sign=((f->V(1)->T().P()-f->V(0)->T().P()) ^ (f->V(2)->T().P()-f->V(0)->T().P()));
903     if (sign!=0) {
904       if (sign*lastsign<0) return false;
905       lastsign=sign;
906     }
907   }
908   return true;
909 }
910 
911 // detects and marks folded faces, by setting their quality to 0 (or 1 otherwise)
912 // returns number of folded faces
913 template<class MESH_TYPE>
MarkTexCoordFolds(MESH_TYPE & m)914 int MarkTexCoordFolds(MESH_TYPE &m){
915 
916   assert(HasPerVertexTexCoord(m));
917   assert(m.HasPerFaceQuality());
918 
919   typedef typename MESH_TYPE::VertexType::TexCoordType::PointType::ScalarType ScalarType;
920 
921   SimpleTempData<typename MESH_TYPE::FaceContainer, short> sign(m.face);
922   //sign.Start(0);
923 
924   // first pass, determine predominant sign
925   int npos=0, nneg=0;
926   ScalarType lastsign=0;
927   for (typename MESH_TYPE::FaceIterator f=m.face.begin(); f!=m.face.end(); f++){
928     ScalarType fsign=((f->V(1)->T().P()-f->V(0)->T().P()) ^ (f->V(2)->T().P()-f->V(0)->T().P()));
929     if (fsign<0) { sign[f]=-1;  nneg++; }
930     if (fsign>0) { sign[f]=+1; npos++; }
931   }
932 
933   // second pass, detect folded faces
934   int res=0;
935   short gsign= (nneg>npos)?-1:+1;
936   for (typename MESH_TYPE::FaceIterator f=m.face.begin(); f!=m.face.end(); f++){
937     if (sign[f]*gsign<0){
938       res++;
939       f->Q()=0;
940     } else f->Q()=1;
941   }
942 
943   //sign.Stop();
944 
945   return res;
946 }
947 
948 // Smooths texture coords.
949 // (can be useful to remove folds,
950 //  e.g. these created when obtaining tecture coordinates after projections)
951 template<class MESH_TYPE>
SmoothTexCoords(MESH_TYPE & m)952 void SmoothTexCoords(MESH_TYPE &m){
953 
954   assert(HasPerVertexTexCoord(m));
955 
956   typedef typename MESH_TYPE::VertexType::TexCoordType::PointType PointType;
957   typedef typename MESH_TYPE::VertexType::TexCoordType::ScalarType ScalarType;
958 
959   SimpleTempData<typename MESH_TYPE::VertContainer, int> div(m.vert);
960   SimpleTempData<typename MESH_TYPE::VertContainer, PointType > sum(m.vert);
961 
962  /* div.Start();
963   sum.Start();*/
964 
965 	for (typename MESH_TYPE::VertexIterator v=m.vert.begin(); v!=m.vert.end(); v++) {
966 		sum[v].SetZero();
967     div[v]=0;
968 	}
969 
970 	for (typename MESH_TYPE::FaceIterator f=m.face.begin(); f!=m.face.end(); f++){
971 		div[f->V(0)] +=2; sum[f->V(0)] += f->V(2)->T().P(); sum[f->V(0)] += f->V(1)->T().P();
972 		div[f->V(1)] +=2; sum[f->V(1)] += f->V(0)->T().P(); sum[f->V(1)] += f->V(2)->T().P();
973 		div[f->V(2)] +=2; sum[f->V(2)] += f->V(1)->T().P(); sum[f->V(2)] += f->V(0)->T().P();
974 	}
975 
976 	for (typename MESH_TYPE::VertexIterator v=m.vert.begin(); v!=m.vert.end(); v++)
977 	if (!v->IsB())
978   {
979 		//if (v->div>0) {
980 	    if (div[v]>0) {
981 			v->T().P() = sum[v]/((ScalarType)div[v]);
982 		}
983 	}
984 
985 	/*div.Stop();
986   sum.Stop();*/
987 
988 }
989 // MIPSTexCoordFoldHealer
990 // ----------------------
991 // Uses MIPS optimization to attempt to remove folds.
992 // Acts only in proximity of foleded faces!
993 // Use "iterateUntilConvergence" to unfold faces (returns number of performed iterations, as usual)
994 // Use "maxStarSize" (direct access) to determine size of affected patch around folded face
995 
996 // AUTO_SPEED metaparameter:
997 #define AUTO_SPEED 1
998 // if set to one, speed is reduced/increase automatically to avoid oscillating behaviour
999 // (consumes memory and CPU, but increase robustness with speed parameter and sometimes converge faster)
1000 
1001 template<class MESH_TYPE>
1002 class MIPSTexCoordFoldHealer:public MIPSTexCoordOptimization<MESH_TYPE>{
1003 public:
1004 
1005   int maxStarSize; // max star size that is affected around a folded face.. Default: 3
1006 
1007   typedef MESH_TYPE MeshType;
1008   typedef typename MESH_TYPE::VertexIterator VertexIterator;
1009   typedef typename MESH_TYPE::FaceIterator FaceIterator;
1010   typedef typename MESH_TYPE::VertexType VertexType;
1011   typedef typename MESH_TYPE::FaceType FaceType;
1012   typedef typename MESH_TYPE::ScalarType ScalarType;
1013 
1014   typedef MIPSTexCoordOptimization<MESH_TYPE> Super; // superclass (commodity)
1015 protected:
1016 
1017   SimpleTempData<typename MESH_TYPE::FaceContainer, bool > foldf;
1018   SimpleTempData<typename MESH_TYPE::VertContainer, bool > foldv;
1019 
1020 #if AUTO_SPEED
1021   SimpleTempData<typename MESH_TYPE::VertContainer, Point2<ScalarType> > lastDir;
1022   SimpleTempData<typename MESH_TYPE::VertContainer, ScalarType > lastSpeed;
1023 #endif
1024 
1025   ScalarType sign;
1026   int nfolds;
1027   FaceType* aFoldedFace;
1028 
PropagateFoldF()1029   void PropagateFoldF(){
1030     for (typename MeshType::FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++){
1031       if (foldv[ f->V(0)] || foldv[ f->V(1)] || foldv[ f->V(2) ] ) {
1032         foldf[f] = true;
1033       }
1034     }
1035   }
1036 
PropagateFoldV()1037   void PropagateFoldV(){
1038     for (typename MeshType::FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++) {
1039       if (foldf[f] ) {
1040         foldv[ f->V(0) ] = foldv[ f->V(1) ] = foldv[ f->V(2) ] = true;
1041       }
1042     }
1043   }
1044 
FindFolds()1045   bool FindFolds(){
1046 
1047     /*ScalarType lastsign=0;*/
1048     int npos=0, nneg=0;
1049     for (typename MESH_TYPE::FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++){
1050       ScalarType sign=((f->V(1)->T().P()-f->V(0)->T().P()) ^ (f->V(2)->T().P()-f->V(0)->T().P()));
1051       if (sign>0) { npos++; }
1052       if (sign<0) { nneg++; }
1053     }
1054     if (npos*nneg==0)     {sign=0; nfolds=0;} else
1055     if (npos>nneg) { sign=+1; nfolds=nneg; } else
1056     { sign=-1; nfolds=npos; };
1057 
1058     for (typename MeshType::FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++){
1059       ScalarType signf=((f->V(1)->T().P()-f->V(0)->T().P()) ^ (f->V(2)->T().P()-f->V(0)->T().P()));
1060 	 /* if ((!Super::isFixed[f->V(0)])&&(!Super::isFixed[f->V(1)])&&(!Super::isFixed[f->V(2)]))*/
1061 		foldf[f] = (signf*sign<0);
1062 	 /* else
1063 		foldf[f] =  false;*/
1064     }
1065 
1066     return true;
1067   }
1068 
1069 public:
1070 
1071  int IterateUntilConvergence(ScalarType threshold=0.0001, int maxite=50){
1072 	(void)threshold;
1073 	  for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++) foldv[v]=false;
1074     FindFolds();
1075     PropagateFoldV();
1076     PropagateFoldF();
1077    /* int i=0;*/
1078     int nite = 0, totIte=0, pass=0;
1079     while (Iterate()>0) {
1080       totIte++;
1081 	  nite++;
1082       if (nite>=maxite) {
1083         PropagateFoldV();
1084         PropagateFoldF();
1085         nite=0;
1086         if (pass++>=maxStarSize) break; // number of passes
1087       }
1088     }
1089     return totIte;
1090   }
1091 
1092 
1093   // constructor and destructor
MIPSTexCoordFoldHealer(MeshType & _m)1094   MIPSTexCoordFoldHealer(MeshType &_m):MIPSTexCoordOptimization<MeshType>(_m),foldf(_m.face),foldv(_m.vert)
1095 #if AUTO_SPEED
1096                    ,lastDir(_m.vert),lastSpeed(_m.vert,1.0)
1097 #endif
1098   {sign=0; nfolds=0;  maxStarSize=3; };
1099 
~MIPSTexCoordFoldHealer()1100   ~MIPSTexCoordFoldHealer(){
1101    /* data.Stop();
1102     sum.Stop();
1103     Super::isFixed.Stop();*/
1104   }
1105 
Iterate()1106   ScalarType Iterate(){
1107 
1108     #define v0 (f->V(0)->T().P())
1109     #define v1 (f->V(1)->T().P())
1110     #define v2 (f->V(2)->T().P())
1111     #define vi (f->V(i)->T().P())
1112     #define vj (f->V(j)->T().P())
1113     #define vk (f->V(k)->T().P())
1114 	  for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++) {
1115 		  //sum[v].Zero();
1116 		  Super::sum[v]=Point2<ScalarType>(0,0);
1117 	  }
1118 
1119     ScalarType totProjArea=0;
1120     nfolds=0;
1121 	  for (FaceIterator f=Super::m.face.begin(); f!=Super::m.face.end(); f++) {
1122       if (Super::isFixed[f->V(0)] && Super::isFixed[f->V(1)] && Super::isFixed[f->V(2)]) continue;
1123       if (!foldf[f]) continue;
1124       ScalarType area2 = ((v1-v0) ^ (v2-v0));
1125       if (area2*sign<0) nfolds++;
1126       totProjArea+=area2;
1127       ScalarType o[3] = { // (opposite edge)^2
1128         (   v1-v2).SquaredNorm(),
1129         (v0   -v2).SquaredNorm(),
1130         (v0-v1   ).SquaredNorm(),
1131       };
1132       ScalarType e =( Super::data[f][0] * o[0] +
1133                       Super::data[f][1] * o[1] +
1134                       Super::data[f][2] * o[2] ) / (area2*area2);
1135 
1136 		  for (int i=0; i<3; i++){
1137         int j=(i+1)%3, k=(i+2)%3;
1138 			  ScalarType p=(vj-vi)*(vk-vi);
1139 				ScalarType gy= e*(o[k]-p) - 2*Super::data[f][j];
1140 				ScalarType gx= e*(o[j]-p) - 2*Super::data[f][k];
1141 
1142 				// speed free mode: (a try!)
1143 			  //sum[f->V(i)]+= ( (vj-vi) * gx + (vk-vi) * gy );// / area2;
1144 
1145 			  // speed mode:
1146         Super::sum[f->V(i)]+= ( (vj-vi) * gx + (vk-vi) * gy ) / area2;
1147 		  }
1148 	  }
1149 
1150     if (nfolds==0) return 0;
1151 
1152  	  for (VertexIterator v=Super::m.vert.begin(); v!=Super::m.vert.end(); v++)
1153     if (  !Super::isFixed[v] && foldv[v] )
1154     {
1155       ScalarType n=Super::sum[v].Norm(); if ( n > 1 ) { Super::sum[v]/=n; n=1.0;}
1156 #if AUTO_SPEED
1157       if (Super::sum[v]*lastDir[v] < 0.0) lastSpeed[v]*=(ScalarType)0.8; else lastSpeed[v]*=(ScalarType)1.1;
1158       lastDir[v]=Super::sum[v];
1159 		  v->T().P()-=(Super::sum[v] ) * (Super::speed * lastSpeed[v] );
1160 #else
1161 		  v->T().P()-=(Super::sum[v] ) * Super::speed;
1162 #endif
1163   	}
1164   	return (ScalarType)nfolds;
1165   	#undef v0
1166     #undef v1
1167     #undef v2
1168   	#undef vi
1169     #undef vj
1170     #undef vk
1171   	//printf("rejected %d\n",rejected);
1172   }
1173 
1174 
1175 };
1176 
1177 
1178 }	}	// End namespace vcg::tri
1179 
1180 #endif //  __VCGLIB__TEXTCOOORD_OPTIMIZATION
1181