1 /****************************************************************************
2 * VCGLib o o *
3 * Visual and Computer Graphics Library o o *
4 * _ O _ *
5 * Copyright(C) 2004-2016 \/)\/ *
6 * Visual Computing Lab /\/| *
7 * ISTI - Italian National Research Council | *
8 * \ *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
20 * for more details. *
21 * *
22 ****************************************************************************/
23 #include <vcg/math/histogram.h>
24 #include <vcg/complex/algorithms/update/curvature.h>
25 #include <vcg/complex/algorithms/update/flag.h>
26 #include <algorithm>
27
28 #ifndef VCG_TANGENT_FIELD_OPERATORS
29 #define VCG_TANGENT_FIELD_OPERATORS
30
31 namespace vcg {
32 namespace tri{
33
34
35
36
37 template < typename ScalarType >
InterpolateNRosy2D(const std::vector<vcg::Point2<ScalarType>> & V,const std::vector<ScalarType> & W,const int N)38 vcg::Point2<ScalarType> InterpolateNRosy2D(const std::vector<vcg::Point2<ScalarType> > &V,
39 const std::vector<ScalarType> &W,
40 const int N)
41 {
42 // check parameter
43 assert(V.size() == W.size());
44 assert( N > 0);
45
46 // create a vector of angles
47 std::vector<ScalarType> angles(V.size(), 0);
48
49 // make angle mod 2pi/N by multiplying times N
50 for (size_t i = 0; i < V.size(); i++)
51 angles[i] = std::atan2(V[i].Y(), V[i].X() ) * N;
52
53 // create vector of directions
54 std::vector<vcg::Point2<ScalarType> > VV(V.size(), vcg::Point2<ScalarType>(0,0));
55
56 // compute directions
57 for (size_t i = 0; i < V.size(); i++) {
58 VV[i].X() = std::cos(angles[i]);
59 VV[i].Y() = std::sin(angles[i]);
60 }
61
62 // average vector
63 vcg::Point2<ScalarType> Res(0,0);
64
65 // compute average of the unit vectors
66 ScalarType Sum=0;
67 for (size_t i = 0; i < VV.size(); i++)
68 {
69 Res += VV[i] * W[i];
70 Sum+=W[i];
71 }
72 assert(Sum>0);
73 Res /=Sum;
74
75 //R /= VV.rows();
76
77 // scale them back
78 ScalarType a = std::atan2(Res.Y(), Res.X()) / N;
79 Res.X() = std::cos(a);
80 Res.Y() = std::sin(a);
81
82 return Res;
83 }
84
85 template < typename ScalarType >
InterpolateNRosy3D(const std::vector<vcg::Point3<ScalarType>> & V,const std::vector<vcg::Point3<ScalarType>> & Norm,const std::vector<ScalarType> & W,const int N,const vcg::Point3<ScalarType> & TargetN)86 vcg::Point3<ScalarType> InterpolateNRosy3D(const std::vector<vcg::Point3<ScalarType> > &V,
87 const std::vector<vcg::Point3<ScalarType> > &Norm,
88 const std::vector<ScalarType> &W,
89 const int N,
90 const vcg::Point3<ScalarType> &TargetN)
91 {
92 typedef typename vcg::Point3<ScalarType> CoordType;
93 ///create a reference frame along TargetN
94 CoordType TargetZ=TargetN;
95 TargetZ.Normalize();
96 CoordType U=CoordType(1,0,0);
97 if (fabs(TargetZ*U)>0.999)
98 U=CoordType(0,1,0);
99
100 CoordType TargetX=TargetZ^U;
101 CoordType TargetY=TargetX^TargetZ;
102 TargetX.Normalize();
103 TargetY.Normalize();
104 vcg::Matrix33<ScalarType> RotFrame=vcg::TransformationMatrix(TargetX,TargetY,TargetZ);
105 vcg::Matrix33<ScalarType> RotFrameInv=vcg::Inverse(RotFrame);
106 std::vector<vcg::Point2<ScalarType> > Cross2D;
107 ///rotate each vector to transform to 2D
108 for (size_t i=0;i<V.size();i++)
109 {
110 CoordType NF=Norm[i];
111 NF.Normalize();
112 CoordType Vect=V[i];
113 Vect.Normalize();
114 //ScalarType Dot=fabs(Vect*NF);
115 //std::cout << "V[i] " << V[i].X() << " " << V[i].Y() << std::endl << std::flush;
116
117 ///rotate the vector to become tangent to the reference plane
118 vcg::Matrix33<ScalarType> RotNorm=vcg::RotationMatrix(Norm[i],TargetN);
119 //std::cout << "Norm[i] " << Norm[i].X() << " " << Norm[i].Y() << " " << Norm[i].Z()<< std::endl;
120 //std::cout << "TargetN " << TargetN.X() << " " << TargetN.Y() << " " << TargetN.Z()<< std::endl<< std::flush;
121
122 CoordType rotV=RotNorm*V[i];
123 //assert(fabs(rotV*TargetN)<0.000001);
124 rotV.Normalize();
125 //std::cout << "rotV " << rotV.X() << " " << rotV.Y() << " " << rotV.Z()<< std::endl<< std::flush;
126
127 ///trassform to the reference frame
128 rotV=RotFrame*rotV;
129 assert(!isnan(rotV.X()));
130 assert(!isnan(rotV.Y()));
131
132 //it's 2D from now on
133 Cross2D.push_back(vcg::Point2<ScalarType>(rotV.X(),rotV.Y()));
134
135 }
136
137 vcg::Point2<ScalarType> AvDir2D=InterpolateNRosy2D(Cross2D,W,N);
138 assert(!isnan(AvDir2D.X()));
139 assert(!isnan(AvDir2D.Y()));
140 CoordType AvDir3D=CoordType(AvDir2D.X(),AvDir2D.Y(),0);
141 //transform back to 3D
142 AvDir3D=RotFrameInv*AvDir3D;
143 return AvDir3D;
144 }
145
146 template <class MeshType>
147 class CrossField
148 {
149 typedef typename MeshType::FaceType FaceType;
150 typedef typename MeshType::VertexType VertexType;
151 typedef typename MeshType::CoordType CoordType;
152 typedef typename MeshType::ScalarType ScalarType;
153 typedef typename vcg::face::Pos<FaceType> PosType;
154 typedef typename vcg::Triangle3<ScalarType> TriangleType;
155
156 private:
157
SubDivideDir(const CoordType & Edge0,const CoordType & Edge1,std::vector<CoordType> & SubDEdges,int Nsub)158 static void SubDivideDir(const CoordType &Edge0,
159 const CoordType &Edge1,
160 std::vector<CoordType> &SubDEdges,
161 int Nsub)
162 {
163 CoordType Dir0=Edge0;
164 CoordType Dir1=Edge1;
165 ScalarType a=Edge0.Norm();
166 Dir0.Normalize();
167 Dir1.Normalize();
168 /*
169 *
170 *
171 * A
172 * / |
173 * Dir1 / |
174 * ^ / |
175 * | c / |
176 * / / |b
177 * | / |
178 * / / |
179 * / |
180 * / |
181 * B____________________C ->Dir0
182 * a
183 *
184 * */
185 ScalarType BTotal=vcg::Angle(Dir0,Dir1);
186 //get angle step
187 ScalarType StepAngle=BTotal/(ScalarType)Nsub;
188 //get the other edge
189 CoordType Edge2=Edge1-Edge0;
190
191 //and its direction
192 CoordType Dir2=Edge2;
193 Dir2.Normalize();
194 //find angle C
195 ScalarType C=vcg::Angle(Dir2,-Dir0);
196 //safety checks
197 assert(BTotal<=(M_PI));
198 assert(BTotal>=0);
199 assert(C<=(M_PI));
200 assert(C>=0);
201
202 //push the first one
203 SubDEdges.push_back(Edge0);
204 for (size_t i=1;i<Nsub;i++)
205 {
206 //find angle interval
207 ScalarType B=StepAngle*(ScalarType)i;
208 ScalarType A=(M_PI-C-B);
209 assert(A<=(M_PI));
210 assert(A>=0);
211 //find current lenght
212 ScalarType b=a*sin(B)/sin(A);
213 //then move along the direction of edge b
214 CoordType intervB=Dir2*b;
215 //finally find absolute position summing edge 0
216 intervB+=Edge0;
217 SubDEdges.push_back(intervB);
218 }
219 //push the last one
220 SubDEdges.push_back(Edge1);
221 }
222
FindSubDir(vcg::Triangle3<ScalarType> T3,size_t Nvert,std::vector<CoordType> & SubDEdges,int Nsub)223 static void FindSubDir(vcg::Triangle3<ScalarType> T3,
224 size_t Nvert,
225 std::vector<CoordType> &SubDEdges,
226 int Nsub)
227 {
228 CoordType P0=T3.P0(Nvert);
229 CoordType P1=T3.P1(Nvert);
230 CoordType P2=T3.P2(Nvert);
231 P1-=P0;
232 P2-=P0;
233 SubDivideDir(P1,P2,SubDEdges,Nsub);
234 for (size_t j=0;j<SubDEdges.size();j++)
235 SubDEdges[j]+=P0;
236 }
237
SubdivideTris(vcg::Triangle3<ScalarType> T3,size_t Nvert,std::vector<vcg::Triangle3<ScalarType>> & SubTris,int Nsub)238 static void SubdivideTris(vcg::Triangle3<ScalarType> T3,
239 size_t Nvert,
240 std::vector<vcg::Triangle3<ScalarType> > &SubTris,
241 int Nsub)
242 {
243 std::vector<CoordType> SplittedPos;
244 FindSubDir(T3,Nvert,SplittedPos,Nsub);
245 CoordType P0=T3.P(Nvert);
246 //then create the triangles
247 for (size_t j=0;j<SplittedPos.size()-1;j++)
248 {
249 TriangleType T(P0,SplittedPos[j+1],SplittedPos[j]);
250 SubTris.push_back(T);
251 }
252 }
253
Sign(ScalarType a)254 static ScalarType Sign(ScalarType a){return (ScalarType)((a>0)?+1:-1);}
255
256 static void FindSubTriangles(const typename vcg::face::Pos<FaceType> &vPos,
257 std::vector<TriangleType> &SubFaces,
258 std::vector<FaceType*> &OriginalFace,
259 int numSub=3)
260 {
261 SubFaces.clear();
262 OriginalFace.clear();
263
264 //not ct on border
265 assert(!vPos.IsBorder());
266 //the vertex should be the one on the pos
267 assert(vPos.F()->V(vPos.E())==vPos.V());
268
269 // collect all faces around the star of the vertex
270 std::vector<PosType> StarPos;
271 vcg::face::VFOrderedStarFF(vPos, StarPos);
272
273 //all the infos for the strip of triangles
274
275 VertexType* v0=vPos.V();
276 CoordType P0=v0->P();
277 //create the strip of triangles
278 for (size_t i = 0; i < StarPos.size(); ++i)
279 {
280 PosType currPos=StarPos[i];
281 FaceType *CurrF=currPos.F();
282 OriginalFace.push_back(CurrF);
283
284 VertexType *v1=CurrF->V2(currPos.E());
285 VertexType *v2=CurrF->V1(currPos.E());
286 CoordType P1=v1->P();
287 CoordType P2=v2->P();
288 assert(v1!=v2);
289 assert(v0!=v1);
290 assert(v0!=v2);
291
292 SubdivideTris(vcg::Triangle3<ScalarType>(P0,P1,P2),0,SubFaces,numSub);
293 }
294 }
295
InterpolateCrossSubTriangles(const std::vector<TriangleType> & SubFaces,const std::vector<FaceType * > & OriginalFace,std::vector<CoordType> & Dir,std::vector<CoordType> & NormSubF)296 static void InterpolateCrossSubTriangles(const std::vector<TriangleType> &SubFaces,
297 const std::vector<FaceType*> &OriginalFace,
298 std::vector<CoordType> &Dir,
299 std::vector<CoordType> &NormSubF)
300 {
301 Dir.clear();
302 //check
303 assert(SubFaces.size()>OriginalFace.size());
304
305 int SubDivisionSize=SubFaces.size()/OriginalFace.size();
306 assert(SubDivisionSize>=3);
307 assert(SubDivisionSize%2==1);//should be odd
308 int MiddlePos=SubDivisionSize/2+1;//the one in the middle that should preserve face's cross field
309
310 //calculate the interpolation weights and intervals
311 std::vector<std::pair<FaceType*,FaceType*> > InterpFaces;
312 std::vector<ScalarType> InterpWeights;
313
314 //the one in the middle should spond to one of the
315 //two faces the rest is interpolated
316 for (size_t i=0;i<SubFaces.size();i++)
317 {
318 //find the interval and get the index in the sub_interval
319 int interval=i/SubDivisionSize;
320 int sub_int=i%SubDivisionSize;
321
322 int IndexF0=-1,IndexF1=-1;
323 ScalarType alpha;
324
325 //get the index of the two faces upon which to interpolate
326 if (sub_int<MiddlePos)
327 {
328 IndexF0=(interval+(OriginalFace.size()-1)) % OriginalFace.size();
329 IndexF1=interval;
330 alpha=1-(ScalarType)(sub_int+MiddlePos)/(ScalarType)SubDivisionSize;
331 }
332 else
333 if (sub_int>MiddlePos)
334 {
335 IndexF0=interval;
336 IndexF1=(interval+1) % OriginalFace.size();
337 alpha=1-(sub_int-MiddlePos)/(ScalarType)SubDivisionSize;
338 }
339 else //sub_int==MiddlePos
340 {
341 IndexF0=interval;
342 IndexF1=(interval+1) % OriginalFace.size();
343 alpha=1;
344 }
345 assert((IndexF0>=0)&&(IndexF0<OriginalFace.size()));
346 assert((IndexF1>=0)&&(IndexF1<OriginalFace.size()));
347
348 FaceType* F0=OriginalFace[IndexF0];
349 FaceType* F1=OriginalFace[IndexF1];
350
351 InterpFaces.push_back(std::pair<FaceType*,FaceType*>(F0,F1));
352 InterpWeights.push_back(alpha);
353 }
354 assert(InterpFaces.size()==InterpWeights.size());
355 //then calculate the interpolated cross field
356 for (size_t i=0;i<InterpFaces.size();i++)
357 {
358 std::vector<vcg::Point3<ScalarType> > TangVect;
359 std::vector<vcg::Point3<ScalarType> > Norms;
360 std::vector<ScalarType> W;
361 Norms.push_back(InterpFaces[i].first->N());
362 Norms.push_back(InterpFaces[i].second->N());
363 CoordType Dir0=InterpFaces[i].first->PD1();
364 CoordType Dir1=InterpFaces[i].second->PD1();
365 TangVect.push_back(Dir0);
366 TangVect.push_back(Dir1);
367 ScalarType CurrW=InterpWeights[i];
368 W.push_back(CurrW);
369 W.push_back(1-CurrW);
370
371 CoordType TargetN=InterpFaces[i].first->N();
372 if (CurrW<0.5)
373 TargetN=InterpFaces[i].second->N();
374
375 NormSubF.push_back(TargetN);
376 //CoordType TargetN=vcg::Normal(SubFaces[i].P(0),SubFaces[i].P(1),SubFaces[i].P(2));
377 TargetN.Normalize();
378 CoordType InterpD0=InterpolateNRosy3D(TangVect,Norms,W,4,TargetN);
379 //CoordType InterpD0=Dir1;
380 //if (CurrW>0.5)InterpD0=Dir0;
381
382 Dir.push_back(InterpD0);
383 }
384 }
385
turn(const CoordType & norm,const CoordType & d0,const CoordType & d1)386 static ScalarType turn (const CoordType &norm, const CoordType &d0, const CoordType &d1)
387 {
388 //first check if they are coplanar up to a certain delta
389 return ((d0 ^ d1).normalized() * norm);
390 }
391
InterpolateDir(const CoordType & Dir0,const CoordType & Dir1,const CoordType & Sep0,const CoordType & Sep1,const TriangleType & t0,const TriangleType & t1,CoordType & Interpolated,size_t & Face)392 static void InterpolateDir(const CoordType &Dir0,
393 const CoordType &Dir1,
394 const CoordType &Sep0,
395 const CoordType &Sep1,
396 const TriangleType &t0,
397 const TriangleType &t1,
398 CoordType &Interpolated,
399 size_t &Face)
400 {
401 //find smallest edge
402 ScalarType smallestE=std::numeric_limits<ScalarType>::max();
403 for (int j=0;j<3;j++)
404 {
405 ScalarType L0=(t0.P0(j)-t0.P1(j)).Norm();
406 ScalarType L1=(t1.P0(j)-t1.P1(j)).Norm();
407 if (L0<smallestE) smallestE=L0;
408 if (L1<smallestE) smallestE=L1;
409 }
410
411 //safety check
412 assert(t0.P(0)==t1.P(0));
413
414 CoordType Origin=t0.P(0);
415 TriangleType T0Rot(CoordType(0,0,0),t0.P(1)-Origin,t0.P(2)-Origin);
416 TriangleType T1Rot(CoordType(0,0,0),t1.P(1)-Origin,t1.P(2)-Origin);
417
418 //then rotate normal of T0 to match with normal of T1
419 CoordType N0=vcg::Normal(T0Rot.cP(0),T0Rot.cP(1),T0Rot.P(2));
420 CoordType N1=vcg::Normal(T1Rot.cP(0),T1Rot.cP(1),T1Rot.cP(2));
421 N0.Normalize();
422 N1.Normalize();
423 vcg::Matrix33<ScalarType> rotation=vcg::RotationMatrix(N0,N1);
424
425 //transform the first triangle
426 T0Rot.P(1)=rotation*T0Rot.P(1);
427 T0Rot.P(2)=rotation*T0Rot.P(2);
428
429 //also rotate the directions
430 CoordType Dir0Rot=rotation*Dir0;
431 CoordType Dir1Rot=Dir1;
432 CoordType Sep0Rot=rotation*Sep0;
433 CoordType Sep1Rot=Sep1;
434
435 //find the interpolation angles
436 ScalarType Angle0=vcg::Angle(Dir0Rot,Sep0Rot);
437 ScalarType Angle1=vcg::Angle(Dir1Rot,Sep1Rot);
438 assert(Angle0>=0);
439 assert(Angle1>=0);
440 assert(Angle0<=(M_PI/2));
441 assert(Angle1<=(M_PI/2));
442 ScalarType alpha=0.5;//Angle0/(Angle0+Angle1);
443
444 //then interpolate the direction
445 //CoordType Interp=Dir0Rot*(1-alpha)+Dir1Rot*alpha;
446 Interpolated=Sep0Rot*(1-alpha)+Sep1Rot*alpha;
447 Interpolated.Normalize();
448 Interpolated*=smallestE;
449
450 //then find the triangle which falls into
451 CoordType bary0,bary1;
452 bool Inside0=vcg::InterpolationParameters(T0Rot,Interpolated,bary0);
453 bool Inside1=vcg::InterpolationParameters(T1Rot,Interpolated,bary1);
454 assert(Inside0 || Inside1);
455 // if (!(Inside0 || Inside1))
456 // {
457 // std::cout << "Not Inside " << Interpolated.X() << "," << Interpolated.Y() << "," << Interpolated.Z() << std::endl;
458 // std::cout << "bary0 " << bary0.X() << "," << bary0.Y() << "," << bary0.Z() << std::endl;
459 // std::cout << "bary1 " << bary1.X() << "," << bary1.Y() << "," << bary1.Z() << std::endl;
460 // std::cout << "Diff0 " << fabs(bary0.Norm() - 1) << std::endl;
461 // std::cout << "Diff1 " << fabs(bary1.Norm() - 1) << std::endl;
462 // }
463
464 if (Inside0)
465 {
466 Interpolated=t0.P(0)*bary0.X()+t0.P(1)*bary0.Y()+t0.P(2)*bary0.Z();
467 Interpolated-=Origin;
468 Face=0;
469 }
470 else
471 Face=1;
472
473 //otherwise is already in the tangent space of t0
474 Interpolated.Normalize();
475 }
476
ReduceOneDirectionField(std::vector<CoordType> & directions,std::vector<FaceType * > & faces)477 static void ReduceOneDirectionField(std::vector<CoordType> &directions,
478 std::vector<FaceType*> &faces)
479 {
480 //compute distribution and find closest pait
481 std::vector<ScalarType> DirProd;
482 ScalarType MaxProd=-1;
483 int MaxInd0=-1,MaxInd1=-1;
484 for (size_t i=0;i<directions.size();i++)
485 {
486 ScalarType prod=0;
487 for (size_t j=1;j<directions.size();j++)
488 {
489 int Index0=i;
490 int Index1=(i+j)%directions.size();
491 CoordType Dir0=directions[Index0];
492 CoordType Dir1=directions[Index1];
493 ScalarType currP=(Dir0*Dir1);
494 if (currP>MaxProd)
495 {
496 MaxProd=currP;
497 MaxInd0=Index0;
498 MaxInd1=Index1;
499 }
500 prod+=currP;
501 }
502 DirProd.push_back(prod);
503 }
504 assert(MaxInd0!=MaxInd1);
505
506 //then find the one to be deleted
507 int IndexDel=MaxInd0;
508 if (DirProd[MaxInd1]>DirProd[MaxInd0])IndexDel=MaxInd1;
509
510 std::vector<CoordType> SwapV(directions.begin(),directions.end());
511 std::vector<FaceType*> SwapF(faces.begin(),faces.end());
512
513 directions.clear();
514 faces.clear();
515
516 for (size_t i=0;i<SwapV.size();i++)
517 {
518 if (i==IndexDel)continue;
519 directions.push_back(SwapV[i]);
520 faces.push_back(SwapF[i]);
521 }
522 }
523
524 public:
525
526 static size_t FindSeparatrices(const typename vcg::face::Pos<FaceType> &vPos,
527 std::vector<CoordType> &directions,
528 std::vector<FaceType*> &faces,
529 std::vector<TriangleType> &WrongTris,
530 int expVal=-1,
531 int numSub=3)
532 {
533
534 directions.clear();
535
536 //not ct on border
537 assert(!vPos.IsBorder());
538 //the vertex should be the one on the pos
539 assert(vPos.F()->V(vPos.E())==vPos.V());
540
541 std::vector<TriangleType> SubFaces;
542 std::vector<CoordType> Dir1,Dir2;
543 std::vector<CoordType> Norms;
544 std::vector<FaceType*> OriginalFaces;
545
546 //find subfaces
547 FindSubTriangles(vPos,SubFaces,OriginalFaces,numSub);
548
549 //interpolate the cross field
550 InterpolateCrossSubTriangles(SubFaces,OriginalFaces,Dir1,Norms);
551
552 //then calculate the orthogonal
553 for (size_t i=0;i<Dir1.size();i++)
554 {
555 CoordType TargetN=Norms[i];//vcg::Normal(SubFaces[i].P(0),SubFaces[i].P(1),SubFaces[i].P(2));
556 CoordType CrossDir2=Dir1[i]^TargetN;
557 CrossDir2.Normalize();
558 Dir2.push_back(CrossDir2);
559 }
560
561 //then check the triangles
562 CoordType CenterStar=vPos.V()->P();
563 for (size_t i = 0; i < SubFaces.size(); ++i)
564 {
565 TriangleType t0=SubFaces[i];
566 TriangleType t1=SubFaces[(i+1)%SubFaces.size()];
567 CoordType N0=Norms[i];//vcg::Normal(t0.P(0),t0.P(1),t0.P(2));
568 CoordType N1=Norms[(i+1)%Norms.size()];//vcg::Normal(t1.P(0),t1.P(1),t1.P(2));
569 N0.Normalize();
570 N1.Normalize();
571
572 std::vector<CoordType> SubDEdges0;
573 std::vector<CoordType> SubDEdges1;
574 FindSubDir(t0,0,SubDEdges0,2);
575 FindSubDir(t1,0,SubDEdges1,2);
576 CoordType Bary0=SubDEdges0[1];
577 CoordType Bary1=SubDEdges1[1];
578 //then get the directions to the barycenter
579 Bary0-=CenterStar;
580 Bary1-=CenterStar;
581 Bary0.Normalize();
582 Bary1.Normalize();
583
584 //then get the cross field of the 2 faces
585 CoordType Dir1F0=Dir1[i];
586 CoordType Dir2F0=Dir2[i];
587
588 assert(fabs(Dir1F0*N0)<0.001);
589 assert(fabs(Dir1F0*Dir2F0)<0.001);
590
591 if ((Dir1F0*Bary0)<0)Dir1F0=-Dir1F0;
592 if ((Dir2F0*Bary0)<0)Dir2F0=-Dir2F0;
593
594 CoordType Dir1F1=Dir1[(i+1)%Dir1.size()];
595 CoordType Dir2F1=Dir2[(i+1)%Dir2.size()];
596
597 assert(fabs(Dir1F1*N1)<0.001);
598 assert(fabs(Dir1F1*Dir2F1)<0.01);
599
600 //find the most similar rotation of the cross field
601 vcg::Matrix33<ScalarType> rotation=vcg::RotationMatrix(N0,N1);
602 CoordType Dir1F0R=rotation*Dir1F0;
603 CoordType Dir2F0R=rotation*Dir2F0;
604
605 //then get the closest upf to K*PI/2 rotations
606 Dir1F1=vcg::tri::CrossField<MeshType>::K_PI(Dir1F1,Dir1F0R,N1);
607 Dir2F1=vcg::tri::CrossField<MeshType>::K_PI(Dir2F1,Dir2F0R,N1);
608
609 //then check if cross the direction of the barycenter
610 ScalarType versef0D1 = turn(N0, Bary0, Dir1F0);
611 ScalarType versef1D1 = turn(N1, Bary1, Dir1F1);
612 ScalarType versef0D2 = turn(N0, Bary0, Dir2F0);
613 ScalarType versef1D2 = turn(N1, Bary1, Dir2F1);
614
615 if ((Bary0*Bary1)<0)
616 {
617
618 WrongTris.push_back(t0);
619 WrongTris.push_back(t1);
620 }
621
622
623 CoordType InterpDir;
624 size_t tri_Index=-1;
625 if ((versef0D1 * versef1D1) < ScalarType(0))
626 {
627 InterpolateDir(Dir1F0,Dir1F1,Bary0,Bary1,t0,t1,InterpDir,tri_Index);
628 }
629 else
630 {
631 if ((versef0D2 * versef1D2 )< ScalarType(0))
632 InterpolateDir(Dir2F0,Dir2F1,Bary0,Bary1,t0,t1,InterpDir,tri_Index);
633 }
634 //no separatrix found continue
635 if (tri_Index==-1)continue;
636
637 //retrieve original face
638 assert((tri_Index==0)||(tri_Index==1));
639 int OrigFIndex=((i+tri_Index)%SubFaces.size())/numSub;
640 assert(OrigFIndex>=0);
641 assert(OrigFIndex<OriginalFaces.size());
642
643 FaceType* currF=OriginalFaces[OrigFIndex];
644 //add the data
645 directions.push_back(InterpDir);
646 faces.push_back(currF);
647 }
648 if (expVal==-1)return directions.size();
649 if (directions.size()<=expVal)return directions.size();
650
651 size_t sampledDir=directions.size();
652 int to_erase=directions.size()-expVal;
653 do
654 {
655 ReduceOneDirectionField(directions,faces);
656 to_erase--;
657 }while (to_erase!=0);
658 return sampledDir;
659 }
660
FollowDirection(const FaceType & f0,const FaceType & f1,const CoordType & dir0)661 static CoordType FollowDirection(const FaceType &f0,
662 const FaceType &f1,
663 const CoordType &dir0)
664 {
665 ///first it rotate dir to match with f1
666 CoordType dirR=vcg::tri::CrossField<MeshType>::Rotate(f0,f1,dir0);
667 ///then get the closest upf to K*PI/2 rotations
668 CoordType dir1=f1.cPD1();
669 CoordType ret=vcg::tri::CrossField<MeshType>::K_PI(dir1,dirR,f1.cN());
670 return ret;
671 }
672
FollowDirectionI(const FaceType & f0,const FaceType & f1,const CoordType & dir0)673 static int FollowDirectionI(const FaceType &f0,
674 const FaceType &f1,
675 const CoordType &dir0)
676 {
677 ///first it rotate dir to match with f1
678 CoordType dirTarget=FollowDirection(f0,f1,dir0);
679 CoordType dir[4];
680 CrossVector(f1,dir);
681 ScalarType best=-1;
682 int ret=-1;
683 for (int i=0;i<4;i++)
684 {
685 ScalarType dot=dir[i]*dirTarget;
686 if (dot>best)
687 {
688 best=dot;
689 ret=i;
690 }
691 }
692 assert(ret!=-1);
693
694 return ret;
695 }
696
FollowDirection(const FaceType & f0,const FaceType & f1,int dir0)697 static int FollowDirection(const FaceType &f0,
698 const FaceType &f1,
699 int dir0)
700 {
701 ///first it rotate dir to match with f1
702 CoordType dirS=CrossVector(f0,dir0);
703 CoordType dirR=vcg::tri::CrossField<MeshType>::Rotate(f0,f1,dirS);
704 ///then get the closest upf to K*PI/2 rotations
705 CoordType dir1=f1.cPD1();
706 //int ret=I_K_PI(dir1,dirR,f1.cN());
707 CoordType dir[4];
708 CrossVector(f1,dir);
709 ScalarType best=-1;
710 int ret=-1;
711 for (int i=0;i<4;i++)
712 {
713 ScalarType dot=dir[i]*dirR;
714 if (dot>best)
715 {
716 best=dot;
717 ret=i;
718 }
719 }
720 assert(ret!=-1);
721
722 return ret;
723 }
724
FollowLineDirection(const FaceType & f0,const FaceType & f1,int dir)725 static int FollowLineDirection(const FaceType &f0,
726 const FaceType &f1,
727 int dir)
728 {
729 ///first it rotate dir to match with f1
730 CoordType dir0=CrossVector(f0,dir);
731 CoordType dir0R=vcg::tri::CrossField<MeshType>::Rotate(f0,f1,dir0);
732 ///then get the closest upf to K*PI/2 rotations
733 CoordType dir1_test=CrossVector(f1,dir);
734 CoordType dir2_test=-dir1_test;
735 if ((dir1_test*dir0R)>(dir2_test*dir0R))
736 return dir;
737 return ((dir+2)%4);
738
739 }
740
741 ///fird a tranformation matrix to transform
742 ///the 3D space to 2D tangent space specified
743 ///by the cross field (where Z=0)
TransformationMatrix(const FaceType & f)744 static vcg::Matrix33<ScalarType> TransformationMatrix(const FaceType &f)
745 {
746 typedef typename FaceType::CoordType CoordType;
747 typedef typename FaceType::ScalarType ScalarType;
748
749 ///transform to 3d
750 CoordType axis0=f.cPD1();
751 CoordType axis1=f.cPD2();//axis0^f.cN();
752 CoordType axis2=f.cN();
753
754 return (vcg::TransformationMatrix(axis0,axis1,axis2));
755 }
756
757
758 ///transform a given angle in tangent space wrt X axis of
759 ///tangest space will return the sponding 3D vector
TangentAngleToVect(const FaceType & f,const ScalarType & angle)760 static CoordType TangentAngleToVect(const FaceType &f,const ScalarType &angle)
761 {
762 ///find 2D vector
763 vcg::Point2<ScalarType> axis2D=vcg::Point2<ScalarType>(cos(angle),sin(angle));
764 CoordType axis3D=CoordType(axis2D.X(),axis2D.Y(),0);
765 vcg::Matrix33<ScalarType> Trans=TransformationMatrix(f);
766 vcg::Matrix33<ScalarType> InvTrans=Inverse(Trans);
767 ///then transform
768 return (InvTrans*axis3D);
769 }
770
771 ///find an angle with respect to dirX on the plane perpendiculr to DirZ
772 ///dirX and dirZ should be perpendicular
TangentVectToAngle(const CoordType dirX,const CoordType dirZ,const CoordType & vect3D)773 static ScalarType TangentVectToAngle(const CoordType dirX,
774 const CoordType dirZ,
775 const CoordType &vect3D)
776 {
777 const CoordType dirY=dirX^dirZ;
778 dirX.Normalize();
779 dirY.Normalize();
780 dirZ.Normalize();
781 vcg::Matrix33<ScalarType> Trans=TransformationMatrix(dirX,dirY,dirZ);
782 ///trensform the vector to the reference frame by rotating it
783 CoordType vect_transf=Trans*vect3D;
784
785 ///then put to zero to the Z coordinate
786 vcg::Point2<ScalarType> axis2D=vcg::Point2<ScalarType>(vect_transf.X(),vect_transf.Y());
787 axis2D.Normalize();
788
789 ///then find the angle with respact to axis 0
790 ScalarType alpha=atan2(axis2D.Y(),axis2D.X()); ////to sum up M_PI?
791 if (alpha<0)
792 alpha=(2*M_PI+alpha);
793 if (alpha<0)
794 alpha=0;
795 return alpha;
796 }
797
798 ///find an angle with respect to the tangent frame of given face
VectToAngle(const FaceType & f,const CoordType & vect3D)799 static ScalarType VectToAngle(const FaceType &f,const CoordType &vect3D)
800 {
801 vcg::Matrix33<ScalarType> Trans=TransformationMatrix(f);
802
803 ///trensform the vector to the reference frame by rotating it
804 CoordType vect_transf=Trans*vect3D;
805
806 ///then put to zero to the Z coordinate
807 vcg::Point2<ScalarType> axis2D=vcg::Point2<ScalarType>(vect_transf.X(),vect_transf.Y());
808 axis2D.Normalize();
809
810 ///then find the angle with respact to axis 0
811 ScalarType alpha=atan2(axis2D.Y(),axis2D.X()); ////to sum up M_PI?
812 if (alpha<0)
813 alpha=(2*M_PI+alpha);
814 if (alpha<0)
815 alpha=0;
816 return alpha;
817 }
818
819 ///transform a cross field into a couple of angles
820 static void CrossFieldToAngles(const FaceType &f,
821 ScalarType &alpha1,
822 ScalarType &alpha2,
823 int RefEdge=1)
824 {
825 CoordType axis0=f.cP1(RefEdge)-f.cP0(RefEdge);
826 axis0.Normalize();
827 CoordType axis2=f.cN();
828 axis2.Normalize();
829 CoordType axis1=axis2^axis0;
830 axis1.Normalize();
831
832
833 vcg::Matrix33<ScalarType> Trans=vcg::TransformationMatrix(axis0,axis1,axis2);
834
835 //trensform the vector to the reference frame by rotating it
836 CoordType trasfPD1=Trans*f.cPD1();
837 CoordType trasfPD2=Trans*f.cPD2();
838
839 //then find the angle with respact to axis 0
840 alpha1=atan2(trasfPD1.Y(),trasfPD1.X());
841 alpha2=atan2(trasfPD2.Y(),trasfPD2.X());
842 }
843
844 ///transform a cross field into a couple of angles
845 static void AnglesToCrossField(FaceType &f,
846 const ScalarType &alpha1,
847 const ScalarType &alpha2,
848 int RefEdge=1)
849 {
850 CoordType axis0=f.cP1(RefEdge)-f.cP0(RefEdge);
851 axis0.Normalize();
852 CoordType axis2=f.cN();
853 axis2.Normalize();
854 CoordType axis1=axis2^axis0;
855 axis1.Normalize();
856
857 vcg::Matrix33<ScalarType> Trans=vcg::TransformationMatrix(axis0,axis1,axis2);
858 vcg::Matrix33<ScalarType> InvTrans=Inverse(Trans);
859
860 CoordType PD1=CoordType(cos(alpha1),sin(alpha1),0);
861 CoordType PD2=CoordType(cos(alpha2),sin(alpha2),0);
862
863 //then transform and store in the face
864 f.PD1()=(InvTrans*PD1);
865 f.PD2()=(InvTrans*PD2);
866 }
867
868 ///return the 4 directiona of the cross field in 3D
869 ///given a first direction as input
CrossVector(const CoordType & dir0,const CoordType & norm,CoordType axis[4])870 static void CrossVector(const CoordType &dir0,
871 const CoordType &norm,
872 CoordType axis[4])
873 {
874 axis[0]=dir0;
875 axis[1]=norm^axis[0];
876 axis[2]=-axis[0];
877 axis[3]=-axis[1];
878 }
879
880 ///return the 4 direction in 3D of
881 ///the cross field of a given face
CrossVector(const FaceType & f,CoordType axis[4])882 static void CrossVector(const FaceType &f,
883 CoordType axis[4])
884 {
885 CoordType dir0=f.cPD1();
886 CoordType dir1=f.cPD2();
887 axis[0]=dir0;
888 axis[1]=dir1;
889 axis[2]=-dir0;
890 axis[3]=-dir1;
891 }
892
893 ///return the 4 direction in 3D of
894 ///the cross field of a given face
CrossVector(const VertexType & v,CoordType axis[4])895 static void CrossVector(const VertexType &v,
896 CoordType axis[4])
897 {
898 CoordType dir0=v.cPD1();
899 CoordType dir1=v.cPD2();
900 axis[0]=dir0;
901 axis[1]=dir1;
902 axis[2]=-dir0;
903 axis[3]=-dir1;
904 }
905
906 ///return a specific direction given an integer 0..3
907 ///considering the reference direction of the cross field
CrossVector(const FaceType & f,const int & index)908 static CoordType CrossVector(const FaceType &f,
909 const int &index)
910 {
911 assert((index>=0)&&(index<4));
912 CoordType axis[4];
913 CrossVector(f,axis);
914 return axis[index];
915 }
916
917 ///return a specific direction given an integer 0..3
918 ///considering the reference direction of the cross field
CrossVector(const VertexType & v,const int & index)919 static CoordType CrossVector(const VertexType &v,
920 const int &index)
921 {
922 assert((index>=0)&&(index<4));
923 CoordType axis[4];
924 CrossVector(v,axis);
925 return axis[index];
926 }
927
928 ///set the cross field of a given face
SetCrossVector(FaceType & f,CoordType dir0,CoordType dir1)929 static void SetCrossVector(FaceType &f,
930 CoordType dir0,
931 CoordType dir1)
932 {
933 f.PD1()=dir0;
934 f.PD2()=dir1;
935 }
936
937 ///set the face cross vector from vertex one
SetFaceCrossVectorFromVert(FaceType & f)938 static void SetFaceCrossVectorFromVert(FaceType &f)
939 {
940 const CoordType &t0=f.V(0)->PD1();
941 const CoordType &t1=f.V(1)->PD1();
942 const CoordType &t2=f.V(2)->PD1();
943 const CoordType &N0=f.V(0)->N();
944 const CoordType &N1=f.V(0)->N();
945 const CoordType &N2=f.V(0)->N();
946 const CoordType &NF=f.N();
947 const CoordType bary=CoordType(0.33333,0.33333,0.33333);
948 CoordType tF0,tF1;
949 tF0=InterpolateCrossField(t0,t1,t2,N0,N1,N2,NF,bary);
950 tF1=NF^tF0;
951 tF0.Normalize();
952 tF1.Normalize();
953 SetCrossVector(f,tF0,tF1);
954
955 //then set the magnitudo
956 ScalarType mag1,mag2;
957 for (int i=0;i<3;i++)
958 {
959 vcg::Matrix33<ScalarType> rotN=vcg::RotationMatrix(f.V(i)->N(),f.N());
960 CoordType rotatedDir=rotN*f.V(i)->PD1();
961
962 if (fabs(rotatedDir*tF0)>fabs(rotatedDir*tF1))
963 {
964 mag1+=fabs(f.V(i)->K1());
965 mag2+=fabs(f.V(i)->K2());
966 }
967 else
968 {
969 mag1+=fabs(f.V(i)->K2());
970 mag2+=fabs(f.V(i)->K1());
971 }
972 }
973
974 f.K1()=mag1/(ScalarType)3;
975 f.K2()=mag2/(ScalarType)3;
976
977 }
978
SetFaceCrossVectorFromVert(MeshType & mesh)979 static void SetFaceCrossVectorFromVert(MeshType &mesh)
980 {
981 for (unsigned int i=0;i<mesh.face.size();i++)
982 {
983 FaceType *f=&mesh.face[i];
984 if (f->IsD())continue;
985 SetFaceCrossVectorFromVert(*f);
986 }
987 }
988
989 ///set the face cross vector from vertex one
SetVertCrossVectorFromFace(VertexType & v)990 static void SetVertCrossVectorFromFace(VertexType &v)
991 {
992 std::vector<FaceType *> faceVec;
993 std::vector<int> index;
994 vcg::face::VFStarVF(&v,faceVec,index);
995 std::vector<CoordType> TangVect;
996 std::vector<CoordType> Norms;
997 for (unsigned int i=0;i<faceVec.size();i++)
998 {
999 TangVect.push_back(faceVec[i]->PD1());
1000 Norms.push_back(faceVec[i]->N());
1001 }
1002 std::vector<ScalarType> Weights(TangVect.size(),1.0/(ScalarType)TangVect.size());
1003 CoordType NRef=v.N();
1004 CoordType N0=faceVec[0]->N();
1005 CoordType DirRef=faceVec[0]->PD1();
1006 ///find the rotation matrix that maps between normals
1007 vcg::Matrix33<ScalarType> rotation=vcg::RotationMatrix(N0,NRef);
1008 DirRef=rotation*DirRef;
1009
1010
1011 CoordType tF1=vcg::tri::CrossField<MeshType>::InterpolateCrossField(TangVect,Weights,Norms,NRef);
1012 tF1.Normalize();
1013 CoordType tF2=NRef^tF1;
1014 tF2.Normalize();
1015 v.PD1()=tF1;
1016 v.PD2()=tF2;
1017 }
1018
SetVertCrossVectorFromFace(MeshType & mesh)1019 static void SetVertCrossVectorFromFace(MeshType &mesh)
1020 {
1021 for (unsigned int i=0;i<mesh.vert.size();i++)
1022 {
1023 VertexType *v=&mesh.vert[i];
1024 if (v->IsD())continue;
1025 SetVertCrossVectorFromFace(*v);
1026 }
1027 }
1028
1029 ///rotate a given vector from the tangent space
1030 ///of f0 to the tangent space of f1 by considering the difference of normals
Rotate(const FaceType & f0,const FaceType & f1,const CoordType & dir3D)1031 static CoordType Rotate(const FaceType &f0,const FaceType &f1,const CoordType &dir3D)
1032 {
1033 CoordType N0=f0.cN();
1034 CoordType N1=f1.cN();
1035
1036 ///find the rotation matrix that maps between normals
1037 vcg::Matrix33<ScalarType> rotation=vcg::RotationMatrix(N0,N1);
1038 CoordType rotated=rotation*dir3D;
1039 return rotated;
1040 }
1041
1042 // returns the 90 deg rotation of a (around n) most similar to target b
1043 /// a and b should be in the same plane orthogonal to N
K_PI(const CoordType & a,const CoordType & b,const CoordType & n)1044 static CoordType K_PI(const CoordType &a, const CoordType &b, const CoordType &n)
1045 {
1046 CoordType c = (a^n).normalized();///POSSIBLE SOURCE OF BUG CHECK CROSS PRODUCT
1047 ScalarType scorea = a*b;
1048 ScalarType scorec = c*b;
1049 if (fabs(scorea)>=fabs(scorec)) return a*Sign(scorea); else return c*Sign(scorec);
1050 }
1051
1052 // returns the 90 deg rotation of a (around n) most similar to target b
1053 /// a and b should be in the same plane orthogonal to N
I_K_PI(const CoordType & a,const CoordType & b,const CoordType & n)1054 static int I_K_PI(const CoordType &a, const CoordType &b, const CoordType &n)
1055 {
1056 CoordType c = (n^a).normalized();
1057 ScalarType scorea = a*b;
1058 ScalarType scorec = c*b;
1059 if (fabs(scorea)>=fabs(scorec))///0 or 2
1060 {
1061 if (scorea>0)return 0;
1062 return 2;
1063 }else ///1 or 3
1064 {
1065 if (scorec>0)return 1;
1066 return 3;
1067 }
1068 }
1069
1070 ///interpolate cross field with barycentric coordinates
InterpolateCrossField(const CoordType & t0,const CoordType & t1,const CoordType & t2,const CoordType & n0,const CoordType & n1,const CoordType & n2,const CoordType & target_n,const CoordType & bary)1071 static CoordType InterpolateCrossField(const CoordType &t0,
1072 const CoordType &t1,
1073 const CoordType &t2,
1074 const CoordType &n0,
1075 const CoordType &n1,
1076 const CoordType &n2,
1077 const CoordType &target_n,
1078 const CoordType &bary)
1079 {
1080 std::vector<CoordType > V,Norm;
1081 std::vector<ScalarType > W;
1082 V.push_back(t0);
1083 V.push_back(t1);
1084 V.push_back(t2);
1085 Norm.push_back(n0);
1086 Norm.push_back(n1);
1087 Norm.push_back(n2);
1088 W.push_back(bary.X());
1089 W.push_back(bary.Y());
1090 W.push_back(bary.Z());
1091
1092 CoordType sum=vcg::tri::InterpolateNRosy3D(V,Norm,W,4,target_n);
1093 return sum;
1094 }
1095
1096 ///interpolate cross field with barycentric coordinates using normalized weights
InterpolateCrossField(const std::vector<CoordType> & TangVect,const std::vector<ScalarType> & Weight,const std::vector<CoordType> & Norms,const CoordType & BaseNorm)1097 static CoordType InterpolateCrossField(const std::vector<CoordType> &TangVect,
1098 const std::vector<ScalarType> &Weight,
1099 const std::vector<CoordType> &Norms,
1100 const CoordType &BaseNorm)
1101 {
1102
1103 CoordType sum=InterpolateNRosy3D(TangVect,Norms,Weight,4,BaseNorm);
1104 return sum;
1105 }
1106
1107 ///interpolate cross field with scalar weight
InterpolateCrossFieldLine(const typename FaceType::CoordType & t0,const typename FaceType::CoordType & t1,const typename FaceType::CoordType & n0,const typename FaceType::CoordType & n1,const typename FaceType::CoordType & target_n,const typename FaceType::ScalarType & weight)1108 static typename FaceType::CoordType InterpolateCrossFieldLine(const typename FaceType::CoordType &t0,
1109 const typename FaceType::CoordType &t1,
1110 const typename FaceType::CoordType &n0,
1111 const typename FaceType::CoordType &n1,
1112 const typename FaceType::CoordType &target_n,
1113 const typename FaceType::ScalarType &weight)
1114 {
1115 std::vector<CoordType > V,Norm;
1116 std::vector<ScalarType > W;
1117 V.push_back(t0);
1118 V.push_back(t1);
1119 Norm.push_back(n0);
1120 Norm.push_back(n1);
1121 W.push_back(weight);
1122 W.push_back(1-weight);
1123 InterpolateNRosy3D(V,Norm,&W,4,target_n);
1124 }
1125
1126
1127 ///return the difference of two cross field, values between [0,1]
DifferenceCrossField(const typename FaceType::CoordType & t0,const typename FaceType::CoordType & t1,const typename FaceType::CoordType & n)1128 static typename FaceType::ScalarType DifferenceCrossField(const typename FaceType::CoordType &t0,
1129 const typename FaceType::CoordType &t1,
1130 const typename FaceType::CoordType &n)
1131 {
1132 CoordType trans0=t0;
1133 CoordType trans1=K_PI(t1,t0,n);
1134 ScalarType diff = vcg::AngleN(trans0,trans1)/(M_PI_4);
1135 return diff;
1136 }
1137
1138 ///return the difference of two cross field, values between [0,1]
DifferenceLineField(const typename FaceType::CoordType & t0,const typename FaceType::CoordType & t1,const typename FaceType::CoordType & n)1139 static typename FaceType::ScalarType DifferenceLineField(const typename FaceType::CoordType &t0,
1140 const typename FaceType::CoordType &t1,
1141 const typename FaceType::CoordType &n)
1142 {
1143 CoordType trans0=t0;
1144 CoordType trans1=t1;
1145 if ((trans0*trans1)<0)trans1=-trans1;
1146 ScalarType angleD=vcg::Angle(trans0,trans1);
1147 assert(angleD>=0);
1148 assert(angleD<=M_PI_2);
1149 return (angleD/M_PI_2);
1150 }
1151
1152 ///return the difference of two cross field, values between [0,1]
DifferenceCrossField(const typename vcg::Point2<ScalarType> & t0,const typename vcg::Point2<ScalarType> & t1)1153 static typename FaceType::ScalarType DifferenceCrossField(const typename vcg::Point2<ScalarType> &t0,
1154 const typename vcg::Point2<ScalarType> &t1)
1155 {
1156 CoordType t03D=CoordType(t0.X(),t0.Y(),0);
1157 CoordType t13D=CoordType(t1.X(),t1.Y(),0);
1158 CoordType Norm=CoordType(0,0,1);
1159 // CoordType n=CoordType(0,0,1);
1160 // CoordType trans1=K_PI(t13D,t03D,n);
1161 // ScalarType diff=vcg::AngleN(trans0,trans1)/(M_PI_4);
1162 //ScalarType diff = 1-fabs(trans0*trans1);
1163 return DifferenceCrossField(t03D,t13D,Norm);
1164 }
1165
1166 ///return the difference of two cross field, values between [0,1]
DifferenceLineField(const typename vcg::Point2<ScalarType> & t0,const typename vcg::Point2<ScalarType> & t1)1167 static typename FaceType::ScalarType DifferenceLineField(const typename vcg::Point2<ScalarType> &t0,
1168 const typename vcg::Point2<ScalarType> &t1)
1169 {
1170 CoordType t03D=CoordType(t0.X(),t0.Y(),0);
1171 CoordType t13D=CoordType(t1.X(),t1.Y(),0);
1172 CoordType Norm=CoordType(0,0,1);
1173 // CoordType n=CoordType(0,0,1);
1174 // CoordType trans1=K_PI(t13D,t03D,n);
1175 // ScalarType diff=vcg::AngleN(trans0,trans1)/(M_PI_4);
1176 //ScalarType diff = 1-fabs(trans0*trans1);
1177 return DifferenceLineField(t03D,t13D,Norm);
1178 }
1179
1180 ///compute the mismatch between 2 directions
1181 ///each one si perpendicular to its own normal
MissMatchByCross(const CoordType & dir0,const CoordType & dir1,const CoordType & N0,const CoordType & N1)1182 static int MissMatchByCross(const CoordType &dir0,
1183 const CoordType &dir1,
1184 const CoordType &N0,
1185 const CoordType &N1)
1186 {
1187 CoordType dir0Rot=Rotate(dir0,N0,N1);
1188 CoordType dir1Rot=dir1;
1189
1190 dir0Rot.Normalize();
1191 dir1Rot.Normalize();
1192
1193 ScalarType angle_diff=VectToAngle(dir0Rot,N0,dir1Rot);
1194
1195 ScalarType step=M_PI/2.0;
1196 int i=(int)floor((angle_diff/step)+0.5);
1197 int k=0;
1198 if (i>=0)
1199 k=i%4;
1200 else
1201 k=(-(3*i))%4;
1202 return k;
1203 }
1204
1205 ///compute the mismatch between 2 faces
MissMatchByCross(const FaceType & f0,const FaceType & f1)1206 static int MissMatchByCross(const FaceType &f0,
1207 const FaceType &f1)
1208 {
1209 //CoordType dir0=CrossVector(f0,0);
1210 CoordType dir1=CrossVector(f1,0);
1211
1212 CoordType dir1Rot=Rotate(f1,f0,dir1);
1213 dir1Rot.Normalize();
1214
1215 ScalarType angle_diff=VectToAngle(f0,dir1Rot);
1216
1217 ScalarType step=M_PI/2.0;
1218 int i=(int)floor((angle_diff/step)+0.5);
1219 int k=0;
1220 if (i>=0)
1221 k=i%4;
1222 else
1223 k=(-(3*i))%4;
1224 return k;
1225 }
1226
1227
1228 ///return true if a given vertex is singular,
1229 ///return also the missmatch
IsSingularByCross(const VertexType & v,int & missmatch)1230 static bool IsSingularByCross(const VertexType &v,int &missmatch)
1231 {
1232 typedef typename VertexType::FaceType FaceType;
1233 ///check that is on border..
1234 if (v.IsB())return false;
1235
1236 std::vector<face::Pos<FaceType> > posVec;
1237 //SortedFaces(v,faces);
1238 face::Pos<FaceType> pos(v.cVFp(), v.cVFi());
1239 vcg::face::VFOrderedStarFF(pos, posVec);
1240
1241 int curr_dir=0;
1242 for (unsigned int i=0;i<posVec.size();i++)
1243 {
1244 FaceType *curr_f=posVec[i].F();
1245 FaceType *next_f=posVec[(i+1)%posVec.size()].F();
1246
1247 //find the current missmatch
1248 curr_dir=FollowDirection(*curr_f,*next_f,curr_dir);
1249 }
1250 missmatch=curr_dir;
1251 return(curr_dir!=0);
1252 }
1253
1254 ///select singular vertices
UpdateSingularByCross(MeshType & mesh)1255 static void UpdateSingularByCross(MeshType &mesh)
1256 {
1257 bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular"));
1258 bool hasSingularIndex = vcg::tri::HasPerVertexAttribute(mesh,std::string("SingularIndex"));
1259
1260 typename MeshType::template PerVertexAttributeHandle<bool> Handle_Singular;
1261 typename MeshType::template PerVertexAttributeHandle<int> Handle_SingularIndex;
1262
1263 if (hasSingular)
1264 Handle_Singular=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<bool>(mesh,std::string("Singular"));
1265 else
1266 Handle_Singular=vcg::tri::Allocator<MeshType>::template AddPerVertexAttribute<bool>(mesh,std::string("Singular"));
1267
1268 if (hasSingularIndex)
1269 Handle_SingularIndex=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<int>(mesh,std::string("SingularIndex"));
1270 else
1271 Handle_SingularIndex=vcg::tri::Allocator<MeshType>::template AddPerVertexAttribute<int>(mesh,std::string("SingularIndex"));
1272
1273 for (size_t i=0;i<mesh.vert.size();i++)
1274 {
1275 if (mesh.vert[i].IsD())continue;
1276 //if (mesh.vert[i].IsB())continue;
1277
1278 int missmatch;
1279 if (IsSingularByCross(mesh.vert[i],missmatch))
1280 {
1281 Handle_Singular[i]=true;
1282 Handle_SingularIndex[i]=missmatch;
1283 }
1284 else
1285 {
1286 Handle_Singular[i]=false;
1287 Handle_SingularIndex[i]=0;
1288 }
1289 }
1290 }
1291
1292
GradientToCross(const FaceType & f,const vcg::Point2<ScalarType> & UV0,const vcg::Point2<ScalarType> & UV1,const vcg::Point2<ScalarType> & UV2,CoordType & dirU,CoordType & dirV)1293 static void GradientToCross(const FaceType &f,
1294 const vcg::Point2<ScalarType> &UV0,
1295 const vcg::Point2<ScalarType> &UV1,
1296 const vcg::Point2<ScalarType> &UV2,
1297 CoordType &dirU,
1298 CoordType &dirV)
1299 {
1300 vcg::Point2<ScalarType> Origin2D=(UV0+UV1+UV2)/3;
1301 CoordType Origin3D=(f.cP(0)+f.cP(1)+f.cP(2))/3;
1302
1303 vcg::Point2<ScalarType> UvT0=UV0-Origin2D;
1304 vcg::Point2<ScalarType> UvT1=UV1-Origin2D;
1305 vcg::Point2<ScalarType> UvT2=UV2-Origin2D;
1306
1307 CoordType PosT0=f.cP(0)-Origin3D;
1308 CoordType PosT1=f.cP(1)-Origin3D;
1309 CoordType PosT2=f.cP(2)-Origin3D;
1310
1311 CoordType Bary0,Bary1;
1312 vcg::InterpolationParameters2(UvT0,UvT1,UvT2,vcg::Point2<ScalarType>(1,0),Bary0);
1313 vcg::InterpolationParameters2(UvT0,UvT1,UvT2,vcg::Point2<ScalarType>(0,1),Bary1);
1314
1315 //then transport to 3D
1316 dirU=PosT0*Bary0.X()+PosT1*Bary0.Y()+PosT2*Bary0.Z();
1317 dirV=PosT0*Bary1.X()+PosT1*Bary1.Y()+PosT2*Bary1.Z();
1318
1319 // dirU-=Origin3D;
1320 // dirV-=Origin3D;
1321 dirU.Normalize();
1322 dirV.Normalize();
1323 //orient coherently
1324 CoordType Ntest=dirU^dirV;
1325 CoordType NTarget=vcg::Normal(f.cP(0),f.cP(1),f.cP(2));
1326 if ((Ntest*NTarget)<0)dirV=-dirV;
1327
1328 // //then make them orthogonal
1329 // CoordType dirAvg=dirU^dirV;
1330 CoordType dirVTarget=NTarget^dirU;
1331 CoordType dirUTarget=NTarget^dirV;
1332
1333 dirUTarget.Normalize();
1334 dirVTarget.Normalize();
1335 if ((dirUTarget*dirU)<0)dirUTarget=-dirUTarget;
1336 if ((dirVTarget*dirV)<0)dirVTarget=-dirVTarget;
1337
1338 dirU=(dirU+dirUTarget)/2;
1339 dirV=(dirV+dirVTarget)/2;
1340
1341 dirU.Normalize();
1342 dirV.Normalize();
1343
1344 // ///compute non normalized normal
1345 // CoordType n = f.cN();
1346
1347 // CoordType p0 =f.cP(1) - f.cP(0);
1348 // CoordType p1 =f.cP(2) - f.cP(1);
1349 // CoordType p2 =f.cP(0) - f.cP(2);
1350
1351 // CoordType t[3];
1352 // t[0] = -(p0 ^ n);
1353 // t[1] = -(p1 ^ n);
1354 // t[2] = -(p2 ^ n);
1355
1356 // dirU = t[1]*UV0.X() + t[2]*UV1.X() + t[0]*UV2.X();
1357 // dirV = t[1]*UV0.Y() + t[2]*UV1.Y() + t[0]*UV2.Y();
1358 }
1359
MakeDirectionFaceCoherent(FaceType * f0,FaceType * f1)1360 static void MakeDirectionFaceCoherent(FaceType *f0,
1361 FaceType *f1)
1362 {
1363 CoordType dir0=f0->PD1();
1364 CoordType dir1=f1->PD1();
1365
1366 CoordType dir0Rot=Rotate(*f0,*f1,dir0);
1367 dir0Rot.Normalize();
1368
1369 CoordType targD=K_PI(dir1,dir0Rot,f1->N());
1370
1371 f1->PD1()=targD;
1372 f1->PD2()=f1->N()^targD;
1373 f1->PD2().Normalize();
1374 }
1375
AdjustDirectionsOnTangentspace(MeshType & mesh)1376 static void AdjustDirectionsOnTangentspace(MeshType &mesh)
1377 {
1378 for (size_t i=0;i<mesh.face.size();i++)
1379 {
1380 FaceType *f=&mesh.face[i];
1381 if (f->IsD())continue;
1382 CoordType Ntest=mesh.face[i].PD1()^mesh.face[i].PD2();
1383 Ntest.Normalize();
1384 CoordType Ntarget=mesh.face[i].N();
1385 if ((Ntest*Ntarget)>0.999)continue;
1386
1387 //find the rotation matrix that maps between normals
1388 vcg::Matrix33<ScalarType> rotation=vcg::RotationMatrix(Ntest,Ntarget);
1389 mesh.face[i].PD1()=rotation*mesh.face[i].PD1();
1390 mesh.face[i].PD2()=rotation*mesh.face[i].PD2();
1391 }
1392 }
1393
OrientDirectionFaceCoherently(MeshType & mesh)1394 static void OrientDirectionFaceCoherently(MeshType &mesh)
1395 {
1396 for (size_t i=0;i<mesh.face.size();i++)
1397 {
1398 FaceType *f=&mesh.face[i];
1399 if (f->IsD())continue;
1400 CoordType Ntest= CoordType::Construct( mesh.face[i].PD1()^mesh.face[i].PD2() );
1401 if ((Ntest*vcg::Normal(f->P(0),f->P(1),f->P(2)))<0)mesh.face[i].PD2()=-mesh.face[i].PD2();
1402 }
1403 }
1404
1405 static void MakeDirectionFaceCoherent(MeshType &mesh,
1406 bool normal_diff=true)
1407 {
1408 vcg::tri::UpdateFlags<MeshType>::FaceClearV(mesh);
1409 vcg::tri::UpdateTopology<MeshType>::FaceFace(mesh);
1410
1411 typedef typename std::pair<FaceType*,FaceType*> FacePair;
1412 std::vector<std::pair<ScalarType,FacePair> > heap;
1413
1414 while (true)
1415 {
1416 bool found=false;
1417 for (int i=0; i<(int)mesh.face.size(); i++)
1418 {
1419 FaceType *f=&(mesh.face[i]);
1420 if (f->IsD())continue;
1421 if (f->IsV())continue;
1422 f->SetV();
1423 found=true;
1424
1425 for (int i=0;i<f->VN();i++)
1426 {
1427 FaceType *Fopp=f->FFp(i);
1428 if (Fopp==f)continue;
1429
1430 FacePair entry(f,Fopp);
1431
1432 ScalarType val=0;
1433 if (normal_diff)val=-(f->N()-Fopp->N()).Norm();
1434
1435 heap.push_back(std::pair<ScalarType,FacePair>(val,entry));
1436 }
1437 break;
1438 }
1439
1440 if (!found)
1441 {
1442
1443 vcg::tri::UpdateFlags<MeshType>::FaceClearV(mesh);
1444 return;///all faces has been visited
1445 }
1446
1447 std::make_heap (heap.begin(),heap.end());
1448
1449 while (!heap.empty())
1450 {
1451 std::pop_heap(heap.begin(), heap.end());
1452 FaceType *f0=heap.back().second.first;
1453 FaceType *f1=heap.back().second.second;
1454 assert(f0->IsV());
1455 heap.pop_back();
1456
1457 MakeDirectionFaceCoherent(f0,f1);
1458 f1->SetV();
1459 for (int k=0; k<f1->VN(); k++)
1460 {
1461 FaceType* f2 = f1->FFp(k);
1462 if (f2->IsV())continue;
1463 if (f2->IsD())continue;
1464 if (f2==f1)continue;
1465
1466 FacePair entry(f1,f2);
1467
1468 ScalarType val=0;
1469 if (normal_diff)val=-(f1->N()-f2->N()).Norm();
1470
1471 heap.push_back(std::pair<ScalarType,FacePair>(val,entry));
1472 std::push_heap (heap.begin(),heap.end());
1473 }
1474 }
1475 }
1476
1477 }
1478
1479 ///transform curvature to UV space
1480 static vcg::Point2<ScalarType> CrossToUV(FaceType &f,int numD=0)
1481 {
1482 typedef typename FaceType::ScalarType ScalarType;
1483 typedef typename FaceType::CoordType CoordType;
1484
1485 CoordType Curv=CrossVector(f,numD);
1486 Curv.Normalize();
1487
1488 CoordType bary3d=(f.P(0)+f.P(1)+f.P(2))/3.0;
1489 vcg::Point2<ScalarType> Uv0=f.V(0)->T().P();
1490 vcg::Point2<ScalarType> Uv1=f.V(1)->T().P();
1491 vcg::Point2<ScalarType> Uv2=f.V(2)->T().P();
1492 vcg::Point2<ScalarType> baryUV=(Uv0+Uv1+Uv2)/3.0;
1493 CoordType direct3d=bary3d+Curv;
1494 CoordType baryCoordsUV;
1495 vcg::InterpolationParameters<FaceType,ScalarType>(f,direct3d,baryCoordsUV);
1496 vcg::Point2<ScalarType> curvUV=baryCoordsUV.X()*Uv0+
1497 baryCoordsUV.Y()*Uv1+
1498 baryCoordsUV.Z()*Uv2-baryUV;
1499 curvUV.Normalize();
1500 return curvUV;
1501 }
1502
InitDirFromWEdgeUV(MeshType & mesh)1503 static void InitDirFromWEdgeUV(MeshType &mesh)
1504 {
1505 for (size_t i=0;i<mesh.face.size();i++)
1506 {
1507 vcg::Point2<ScalarType> UV0=mesh.face[i].WT(0).P();
1508 vcg::Point2<ScalarType> UV1=mesh.face[i].WT(1).P();
1509 vcg::Point2<ScalarType> UV2=mesh.face[i].WT(2).P();
1510 GradientToCross(mesh.face[i],UV0,UV1,UV2,
1511 CoordType::Construct(mesh.face[i].PD1()),
1512 CoordType::Construct(mesh.face[i].PD2()) );
1513 }
1514 OrientDirectionFaceCoherently(mesh);
1515 }
1516
1517 };///end class
1518 } //End Namespace Tri
1519 } // End Namespace vcg
1520 #endif
1521