MorphoGraphX
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Curvature.hpp
1 #ifndef CURVATURE_H
2 #define CURVATURE_H
3 
4 #include <Config.hpp>
5 
6 #include <Geometry.hpp>
7 #include <Information.hpp>
8 
9 #include <CImg.h>
10 #include <math.h>
11 #include <string>
12 #include <vector>
13 
14 #define CIMat cimg_library::CImg<float>
15 #define CIMatL cimg_library::CImgList<float>
16 
17 namespace mgx {
18 namespace util {
19 #ifndef isnan
20 using std::isnan;
21 #endif
22 
24 public:
25  virtual ~CurvatureMeasure() {
26  }
27  virtual float measure(float ev1, float ev2) = 0;
28  virtual std::string get_type() = 0;
29  virtual CurvatureMeasure* copy() = 0;
30 };
31 
33 public:
34  float measure(float ev1, float ev2) {
35  return ev1 * ev2;
36  }
37  std::string get_type() {
38  return "gaussian";
39  }
40  virtual CurvatureMeasure* copy() {
41  return new GaussianCurvature();
42  }
43 };
44 
46 public:
47  float measure(float ev1, float ev2) {
48  return (ev1 * ev1) + (ev2 * ev2);
49  }
50  std::string get_type() {
51  return "sumsquare";
52  }
53  virtual CurvatureMeasure* copy() {
54  return new SumSquareCurvature();
55  }
56 };
57 
59 public:
60  float measure(float ev1, float ev2) {
61  return sqrt((ev1 * ev1) + (ev2 * ev2));
62  }
63  std::string get_type() {
64  return "rootsumsquare";
65  }
66  virtual CurvatureMeasure* copy() {
67  return new RootSumSquareCurvature();
68  }
69 };
70 
72 public:
73  float measure(float ev1, float ev2) {
74  return (ev1 + ev2) / 2;
75  }
76  std::string get_type() {
77  return "average";
78  }
79  virtual CurvatureMeasure* copy() {
80  return new AverageCurvature();
81  }
82 };
83 
85 public:
86  float measure(float ev1, float ev2) {
87  return std::min(fabs(ev1), fabs(ev2));
88  }
89  std::string get_type() {
90  return "minimal";
91  }
92  virtual CurvatureMeasure* copy() {
93  return new MinimalCurvature();
94  }
95 };
96 
98 public:
99  float measure(float ev1, float ev2) {
100  return std::max(fabs(ev1), fabs(ev2));
101  }
102  std::string get_type() {
103  return "maximal";
104  }
105  virtual CurvatureMeasure* copy() {
106  return new MaximalCurvature();
107  }
108 };
109 
111 public:
112  float measure(float ev1, float ev2) {
113  return ev1 / ev2;
114  }
115  std::string get_type() {
116  return "anisotropy";
117  }
118  virtual CurvatureMeasure* copy() {
119  return new AnisotropyCurvature();
120  }
121 };
122 
124 public:
125  float measure(float ev1, float ev2) {
126  return (ev1 * ev2 < 0 ? -1.0f : 1.0f) * (fabs(ev1) + fabs(ev2)) / 2.0f;
127  }
128  std::string get_type() {
129  return "signedaverageabs";
130  }
131  virtual CurvatureMeasure* copy() {
132  return new SignedAverageAbsCurvature();
133  }
134 };
135 
136 class Curvature {
137 public:
138  // default constructor
139  Curvature()
140  : W(2, 2)
141  , ev1(0)
142  , ev2(0)
143  {
144  cM = new AverageCurvature();
145  }
146 
147  // convenience constructor
148  Curvature(std::string s)
149  : W(2, 2)
150  , ev1(0)
151  , ev2(0)
152  {
153  if(s == "gaussian")
154  cM = new GaussianCurvature();
155  else if(s == "average")
156  cM = new AverageCurvature();
157  else if(s == "minimal")
158  cM = new MinimalCurvature();
159  else if(s == "maximal")
160  cM = new MaximalCurvature();
161  else if(s == "sumsquare")
162  cM = new SumSquareCurvature();
163  else if(s == "rootsumsquare")
164  cM = new RootSumSquareCurvature();
165  else if(s == "signedaverageabs")
166  cM = new SignedAverageAbsCurvature();
167  else if(s == "anisotropy")
168  cM = new AnisotropyCurvature();
169  else
170  cM = new AverageCurvature();
171  }
172 
173  Curvature(const Curvature& copy)
174  : W(copy.W)
175  , b1(copy.b1)
176  , b2(copy.b2)
177  , b3(copy.b3)
178  , ed1(copy.ed1)
179  , ed2(copy.ed2)
180  , ev1(copy.ev1)
181  , ev2(copy.ev2)
182  , cM(0)
183  {
184  if(copy.cM)
185  cM = copy.cM->copy();
186  }
187 
188  ~Curvature() {
189  delete cM;
190  }
191 
192  // update the curvature object, first entry of pos and norm must be the position, surface normal of the point of
193  // interest
194  void update(const std::vector<Point3f>& pos, const std::vector<Point3f>& nrml)
195  {
196 
197  try {
198  if(nrml.size() > 0 && nrml.size() == pos.size()) {
199 
200  // Calculate local Coordinate system
201  update_local_Basis(nrml[0]);
202 
203  // Set up Least Square Problem
204  int cnt = pos.size();
205  CIMat X = CIMat(7, 3 * (cnt - 1));
206  CIMat b = CIMat(1, 3 * (cnt - 1));
207  Point3f vpos = pos[0];
208  int id = 0;
209 
210  for(int i = 1; i < cnt; i++) {
211  // Transform into local coordinates
212  Point3f vw = pos[i] - vpos;
213  float w1 = vw * b1;
214  float w2 = vw * b2;
215  float w3 = vw * b3;
216 
217  // Transform into local coordinates
218  float n1, n2, n3;
219  n1 = nrml[i] * b1;
220  n2 = nrml[i] * b2;
221  n3 = nrml[i] * b3;
222 
223  X(0, id) = 0.5 * w1 * w1;
224  X(1, id) = w1 * w2;
225  X(2, id) = 0.5 * w2 * w2;
226  X(3, id) = w1 * w1 * w1;
227  X(4, id) = w1 * w1 * w2;
228  X(5, id) = w1 * w2 * w2;
229  X(6, id) = w2 * w2 * w2;
230  b(id) = w3;
231  id++;
232 
233  X(0, id) = w1;
234  X(1, id) = w2;
235  X(2, id) = 0;
236  X(3, id) = 3 * w1 * w1;
237  X(4, id) = 2 * w1 * w2;
238  X(5, id) = w2 * w2;
239  X(6, id) = 0;
240  b(id) = -n1 / n3;
241  id++;
242 
243  X(0, id) = 0;
244  X(1, id) = w1;
245  X(2, id) = w2;
246  X(3, id) = 0;
247  X(4, id) = w1 * w1;
248  X(5, id) = 2 * w1 * w2;
249  X(6, id) = 3 * w2 * w2;
250  b(id) = -n2 / n3;
251  id++;
252  }
253  CIMat Xt = X.get_transpose();
254 
255  CIMat Q = Xt * X;
256 
257  b = Xt * b;
258 
259  float det = Q.det();
260  if(!isnan(det) && !(det == 0)) {
261  // solve least-squares problem
262  b.solve(Q);
263 
264  // update Weingarten-matrix
265  W(0, 0) = b[0];
266  W(1, 0) = W(0, 1) = b[1];
267  W(1, 1) = b[2];
268 
269  // calculate eigenvalues and eigenvectors
270  CIMatL bL = W.get_symmetric_eigen();
271  ev1 = bL(0)(0, 0);
272  ev2 = bL(0)(1, 0);
273  ed1 = -bL(1)(0, 0) * b1 - bL(1)(0, 1) * b2;
274  ed2 = -bL(1)(1, 0) * b1 - bL(1)(1, 1) * b2;
275  if(fabs(ev1) < fabs(ev2)) {
276  using std::swap;
277  swap(ev1, ev2);
278  swap(ed1, ed2);
279  }
280  }
281  }
282  return;
283  }
284  catch(...) {
285  Information::out << "Error calculating curvature" << endl;
286  }
287  }
288 
289  // get principal directions of curvature
290  void get_eigenVectors(Point3f& evec1, Point3f& evec2)
291  {
292  evec1 = ed1;
293  evec2 = ed2;
294  return;
295  }
296 
297  // get principal directions of curvature
298  void get_eigenValues(float& eval1, float& eval2)
299  {
300  eval1 = ev1;
301  eval2 = ev2;
302  return;
303  }
304 
305  // calculate a single-valued measure of curvature (for example "gaussian")
306  float get_curvature() {
307  return cM->measure(ev1, ev2);
308  }
309 
310  // get the name of the curvature measure (for example "gaussian")
311  std::string get_type() {
312  return cM->get_type();
313  }
314 
315  // set the curvature measure
316  void set_Measure(CurvatureMeasure* _cM)
317  {
318  if(cM)
319  delete cM;
320  cM = _cM;
321  return;
322  }
323 
324 private:
325  // calculate local basis with b1,b2, in-plane and b3 out-of-plane basis vectors
326  void update_local_Basis(Point3f nrml)
327  {
328  nrml = nrml / nrml.norm();
329  do {
330  b1 = Point3f(cimg_library::cimg::rand(), cimg_library::cimg::rand(), cimg_library::cimg::rand());
331  b1 = b1 / b1.norm();
332  } while(fabs(b1 * nrml) == 1);
333  b1 = b1 - (b1 * nrml) * nrml;
334  b1 = b1 / b1.norm();
335  b2 = b1.cross(nrml);
336  b3 = nrml;
337 
338  return;
339  }
340 
341  // Weingarten Matrix
342  CIMat W;
343 
344  // Local Basis Vectors
345  Point3f b1;
346  Point3f b2;
347  Point3f b3;
348 
349  // Principal directions of curvature
350  Point3f ed1;
351  Point3f ed2;
352 
353  // Principal values of curvature
354  float ev1;
355  float ev2;
356 
357  CurvatureMeasure* cM;
358 };
359 } // namespace util
360 } // namespace mgx
361 #endif
Definition: Curvature.hpp:84
CU_HOST_DEVICE Vector cross(const Vector &other) const
Compute the cross product as this x other.
Definition: Vector.hpp:641
Definition: Curvature.hpp:97
Definition: Curvature.hpp:45
CU_HOST_DEVICE T norm() const
Euclidean norm of the vector.
Definition: Vector.hpp:522
Definition: Curvature.hpp:32
Definition: Curvature.hpp:136
Common definitions and utilities for all geometry algorithms This file is shared by cuda...
Definition: Curvature.hpp:58
Definition: Curvature.hpp:110
Definition: Curvature.hpp:123
Definition: Curvature.hpp:23
Definition: Curvature.hpp:71