1 // This file is part of libigl, a simple c++ geometry processing library.
2 //
3 // Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.com>, Olga Diamanti <olga.diam@gmail.com>
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public License
6 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
7 // obtain one at http://mozilla.org/MPL/2.0/.
8
9 #include "comb_cross_field.h"
10
11 #include <vector>
12 #include <deque>
13 #include <Eigen/Geometry>
14 #include "per_face_normals.h"
15 #include "is_border_vertex.h"
16 #include "rotation_matrix_from_directions.h"
17
18 #include "triangle_triangle_adjacency.h"
19
20 namespace igl {
21 template <typename DerivedV, typename DerivedF>
22 class Comb
23 {
24 public:
25
26 const Eigen::PlainObjectBase<DerivedV> &V;
27 const Eigen::PlainObjectBase<DerivedF> &F;
28 const Eigen::PlainObjectBase<DerivedV> &PD1;
29 const Eigen::PlainObjectBase<DerivedV> &PD2;
30 DerivedV N;
31
32 private:
33 // internal
34 DerivedF TT;
35 DerivedF TTi;
36
37
38 private:
39
40
Sign(double a)41 static inline double Sign(double a){return (double)((a>0)?+1:-1);}
42
43
44 private:
45
46 // returns the 90 deg rotation of a (around n) most similar to target b
47 /// a and b should be in the same plane orthogonal to N
K_PI_new(const Eigen::Matrix<typename DerivedV::Scalar,3,1> & a,const Eigen::Matrix<typename DerivedV::Scalar,3,1> & b,const Eigen::Matrix<typename DerivedV::Scalar,3,1> & n)48 static inline Eigen::Matrix<typename DerivedV::Scalar, 3, 1> K_PI_new(const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& a,
49 const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& b,
50 const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& n)
51 {
52 Eigen::Matrix<typename DerivedV::Scalar, 3, 1> c = (a.cross(n)).normalized();
53 typename DerivedV::Scalar scorea = a.dot(b);
54 typename DerivedV::Scalar scorec = c.dot(b);
55 if (fabs(scorea)>=fabs(scorec))
56 return a*Sign(scorea);
57 else
58 return c*Sign(scorec);
59 }
60
61
62
63 public:
Comb(const Eigen::PlainObjectBase<DerivedV> & _V,const Eigen::PlainObjectBase<DerivedF> & _F,const Eigen::PlainObjectBase<DerivedV> & _PD1,const Eigen::PlainObjectBase<DerivedV> & _PD2)64 inline Comb(const Eigen::PlainObjectBase<DerivedV> &_V,
65 const Eigen::PlainObjectBase<DerivedF> &_F,
66 const Eigen::PlainObjectBase<DerivedV> &_PD1,
67 const Eigen::PlainObjectBase<DerivedV> &_PD2
68 ):
69 V(_V),
70 F(_F),
71 PD1(_PD1),
72 PD2(_PD2)
73 {
74 igl::per_face_normals(V,F,N);
75 igl::triangle_triangle_adjacency(F,TT,TTi);
76 }
comb(Eigen::PlainObjectBase<DerivedV> & PD1out,Eigen::PlainObjectBase<DerivedV> & PD2out)77 inline void comb(Eigen::PlainObjectBase<DerivedV> &PD1out,
78 Eigen::PlainObjectBase<DerivedV> &PD2out)
79 {
80 // PD1out = PD1;
81 // PD2out = PD2;
82 PD1out.setZero(F.rows(),3);PD1out<<PD1;
83 PD2out.setZero(F.rows(),3);PD2out<<PD2;
84
85 Eigen::VectorXi mark = Eigen::VectorXi::Constant(F.rows(),false);
86
87 std::deque<int> d;
88
89 while (!mark.all()) // Stop until all vertices are marked
90 {
91 int unmarked = 0;
92 while (mark(unmarked))
93 unmarked++;
94
95 d.push_back(unmarked);
96 mark(unmarked) = true;
97
98 while (!d.empty())
99 {
100 int f0 = d.at(0);
101 d.pop_front();
102 for (int k=0; k<3; k++)
103 {
104 int f1 = TT(f0,k);
105 if (f1==-1) continue;
106 if (mark(f1)) continue;
107
108 Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0 = PD1out.row(f0);
109 Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1 = PD1out.row(f1);
110 Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n0 = N.row(f0);
111 Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n1 = N.row(f1);
112
113
114 Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0Rot = igl::rotation_matrix_from_directions(n0, n1)*dir0;
115 dir0Rot.normalize();
116 Eigen::Matrix<typename DerivedV::Scalar, 3, 1> targD = K_PI_new(dir1,dir0Rot,n1);
117
118 PD1out.row(f1) = targD;
119 PD2out.row(f1) = n1.cross(targD).normalized();
120
121 mark(f1) = true;
122 d.push_back(f1);
123
124 }
125 }
126 }
127
128 // everything should be marked
129 for (int i=0; i<F.rows(); i++)
130 {
131 assert(mark(i));
132 }
133 }
134
135
136
137 };
138 }
139 template <typename DerivedV, typename DerivedF>
comb_cross_field(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,const Eigen::PlainObjectBase<DerivedV> & PD1,const Eigen::PlainObjectBase<DerivedV> & PD2,Eigen::PlainObjectBase<DerivedV> & PD1out,Eigen::PlainObjectBase<DerivedV> & PD2out)140 IGL_INLINE void igl::comb_cross_field(const Eigen::PlainObjectBase<DerivedV> &V,
141 const Eigen::PlainObjectBase<DerivedF> &F,
142 const Eigen::PlainObjectBase<DerivedV> &PD1,
143 const Eigen::PlainObjectBase<DerivedV> &PD2,
144 Eigen::PlainObjectBase<DerivedV> &PD1out,
145 Eigen::PlainObjectBase<DerivedV> &PD2out)
146 {
147 igl::Comb<DerivedV, DerivedF> cmb(V, F, PD1, PD2);
148 cmb.comb(PD1out, PD2out);
149 }
150
151 #ifdef IGL_STATIC_LIBRARY
152 // Explicit template instantiation
153 template void igl::comb_cross_field<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
154 template void igl::comb_cross_field<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
155 #endif
156