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