1 // This file is part of libigl, a simple c++ geometry processing library.
2 //
3 // Copyright (C) 2013 Daniele Panozzo <daniele.panozzo@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 #include "principal_curvature.h"
9 #include <iostream>
10 #include <fstream>
11 #include <iomanip>
12 #include <iostream>
13 #include <queue>
14 #include <list>
15 #include <cmath>
16 #include <limits>
17
18 #include <Eigen/SparseCholesky>
19
20 // Lib IGL includes
21 #include <igl/adjacency_list.h>
22 #include <igl/per_face_normals.h>
23 #include <igl/per_vertex_normals.h>
24 #include <igl/avg_edge_length.h>
25 #include <igl/vertex_triangle_adjacency.h>
26
27 typedef enum
28 {
29 SPHERE_SEARCH,
30 K_RING_SEARCH
31 } searchType;
32
33 typedef enum
34 {
35 AVERAGE,
36 PROJ_PLANE
37 } normalType;
38
39 class CurvatureCalculator
40 {
41 public:
42 /* Row number i represents the i-th vertex, whose columns are:
43 curv[i][0] : K2
44 curv[i][1] : K1
45 curvDir[i][0] : PD1
46 curvDir[i][1] : PD2
47 */
48 std::vector< std::vector<double> > curv;
49 std::vector< std::vector<Eigen::Vector3d> > curvDir;
50 bool curvatureComputed;
51 class Quadric
52 {
53 public:
54
Quadric()55 IGL_INLINE Quadric ()
56 {
57 a() = b() = c() = d() = e() = 1.0;
58 }
59
Quadric(double av,double bv,double cv,double dv,double ev)60 IGL_INLINE Quadric(double av, double bv, double cv, double dv, double ev)
61 {
62 a() = av;
63 b() = bv;
64 c() = cv;
65 d() = dv;
66 e() = ev;
67 }
68
a()69 IGL_INLINE double& a() { return data[0];}
b()70 IGL_INLINE double& b() { return data[1];}
c()71 IGL_INLINE double& c() { return data[2];}
d()72 IGL_INLINE double& d() { return data[3];}
e()73 IGL_INLINE double& e() { return data[4];}
74
75 double data[5];
76
evaluate(double u,double v)77 IGL_INLINE double evaluate(double u, double v)
78 {
79 return a()*u*u + b()*u*v + c()*v*v + d()*u + e()*v;
80 }
81
du(double u,double v)82 IGL_INLINE double du(double u, double v)
83 {
84 return 2.0*a()*u + b()*v + d();
85 }
86
dv(double u,double v)87 IGL_INLINE double dv(double u, double v)
88 {
89 return 2.0*c()*v + b()*u + e();
90 }
91
duv(double u,double v)92 IGL_INLINE double duv(double u, double v)
93 {
94 return b();
95 }
96
duu(double u,double v)97 IGL_INLINE double duu(double u, double v)
98 {
99 return 2.0*a();
100 }
101
dvv(double u,double v)102 IGL_INLINE double dvv(double u, double v)
103 {
104 return 2.0*c();
105 }
106
107
fit(const std::vector<Eigen::Vector3d> & VV)108 IGL_INLINE static Quadric fit(const std::vector<Eigen::Vector3d> &VV)
109 {
110 assert(VV.size() >= 5);
111 if (VV.size() < 5)
112 {
113 std::cerr << "ASSERT FAILED! fit function requires at least 5 points: Only " << VV.size() << " were given." << std::endl;
114 exit(0);
115 }
116
117 Eigen::MatrixXd A(VV.size(),5);
118 Eigen::MatrixXd b(VV.size(),1);
119 Eigen::MatrixXd sol(5,1);
120
121 for(unsigned int c=0; c < VV.size(); ++c)
122 {
123 double u = VV[c][0];
124 double v = VV[c][1];
125 double n = VV[c][2];
126
127 A(c,0) = u*u;
128 A(c,1) = u*v;
129 A(c,2) = v*v;
130 A(c,3) = u;
131 A(c,4) = v;
132
133 b(c) = n;
134 }
135
136 sol=A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
137
138 return Quadric(sol(0),sol(1),sol(2),sol(3),sol(4));
139 }
140 };
141
142 public:
143
144 Eigen::MatrixXd vertices;
145 // Face list of current mesh (#F x 3) or (#F x 4)
146 // The i-th row contains the indices of the vertices that forms the i-th face in ccw order
147 Eigen::MatrixXi faces;
148
149 std::vector<std::vector<int> > vertex_to_vertices;
150 std::vector<std::vector<int> > vertex_to_faces;
151 std::vector<std::vector<int> > vertex_to_faces_index;
152 Eigen::MatrixXd face_normals;
153 Eigen::MatrixXd vertex_normals;
154
155 /* Size of the neighborhood */
156 double sphereRadius;
157 int kRing;
158
159 bool localMode; /* Use local mode */
160 bool projectionPlaneCheck; /* Check collected vertices on tangent plane */
161 bool montecarlo;
162 unsigned int montecarloN;
163
164 searchType st; /* Use either a sphere search or a k-ring search */
165 normalType nt;
166
167 double lastRadius;
168 double scaledRadius;
169 std::string lastMeshName;
170
171 /* Benchmark related variables */
172 bool expStep; /* True if we want the radius to increase exponentially */
173 int step; /* If expStep==false, by how much rhe radius increases on every step */
174 int maxSize; /* The maximum limit of the radius in the benchmark */
175
176 IGL_INLINE CurvatureCalculator();
177 IGL_INLINE void init(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F);
178
179 IGL_INLINE void finalEigenStuff(int, const std::vector<Eigen::Vector3d>&, Quadric&);
180 IGL_INLINE void fitQuadric(const Eigen::Vector3d&, const std::vector<Eigen::Vector3d>& ref, const std::vector<int>& , Quadric *);
181 IGL_INLINE void applyProjOnPlane(const Eigen::Vector3d&, const std::vector<int>&, std::vector<int>&);
182 IGL_INLINE void getSphere(const int, const double, std::vector<int>&, int min);
183 IGL_INLINE void getKRing(const int, const double,std::vector<int>&);
184 IGL_INLINE Eigen::Vector3d project(const Eigen::Vector3d&, const Eigen::Vector3d&, const Eigen::Vector3d&);
185 IGL_INLINE void computeReferenceFrame(int, const Eigen::Vector3d&, std::vector<Eigen::Vector3d>&);
186 IGL_INLINE void getAverageNormal(int, const std::vector<int>&, Eigen::Vector3d&);
187 IGL_INLINE void getProjPlane(int, const std::vector<int>&, Eigen::Vector3d&);
188 IGL_INLINE void applyMontecarlo(const std::vector<int>&,std::vector<int>*);
189 IGL_INLINE void computeCurvature();
190 IGL_INLINE void printCurvature(const std::string& outpath);
191 IGL_INLINE double getAverageEdge();
192
rotateForward(double * v0,double * v1,double * v2)193 IGL_INLINE static int rotateForward (double *v0, double *v1, double *v2)
194 {
195 double t;
196
197 if (std::abs(*v2) >= std::abs(*v1) && std::abs(*v2) >= std::abs(*v0))
198 return 0;
199
200 t = *v0;
201 *v0 = *v2;
202 *v2 = *v1;
203 *v1 = t;
204
205 return 1 + rotateForward (v0, v1, v2);
206 }
207
rotateBackward(int nr,double * v0,double * v1,double * v2)208 IGL_INLINE static void rotateBackward (int nr, double *v0, double *v1, double *v2)
209 {
210 double t;
211
212 if (nr == 0)
213 return;
214
215 t = *v2;
216 *v2 = *v0;
217 *v0 = *v1;
218 *v1 = t;
219
220 rotateBackward (nr - 1, v0, v1, v2);
221 }
222
chooseMax(Eigen::Vector3d n,Eigen::Vector3d abc,double ab)223 IGL_INLINE static Eigen::Vector3d chooseMax (Eigen::Vector3d n, Eigen::Vector3d abc, double ab)
224 {
225 int max_i;
226 double max_sp;
227 Eigen::Vector3d nt[8];
228
229 n.normalize ();
230 abc.normalize ();
231
232 max_sp = - std::numeric_limits<double>::max();
233
234 for (int i = 0; i < 4; ++i)
235 {
236 nt[i] = n;
237 if (ab > 0)
238 {
239 switch (i)
240 {
241 case 0:
242 break;
243
244 case 1:
245 nt[i][2] = -n[2];
246 break;
247
248 case 2:
249 nt[i][0] = -n[0];
250 nt[i][1] = -n[1];
251 break;
252
253 case 3:
254 nt[i][0] = -n[0];
255 nt[i][1] = -n[1];
256 nt[i][2] = -n[2];
257 break;
258 }
259 }
260 else
261 {
262 switch (i)
263 {
264 case 0:
265 nt[i][0] = -n[0];
266 break;
267
268 case 1:
269 nt[i][1] = -n[1];
270 break;
271
272 case 2:
273 nt[i][0] = -n[0];
274 nt[i][2] = -n[2];
275 break;
276
277 case 3:
278 nt[i][1] = -n[1];
279 nt[i][2] = -n[2];
280 break;
281 }
282 }
283
284 if (nt[i].dot(abc) > max_sp)
285 {
286 max_sp = nt[i].dot(abc);
287 max_i = i;
288 }
289 }
290 return nt[max_i];
291 }
292
293 };
294
295 class comparer
296 {
297 public:
operator ()(const std::pair<int,double> & lhs,const std::pair<int,double> & rhs) const298 IGL_INLINE bool operator() (const std::pair<int, double>& lhs, const std::pair<int, double>&rhs) const
299 {
300 return lhs.second>rhs.second;
301 }
302 };
303
CurvatureCalculator()304 IGL_INLINE CurvatureCalculator::CurvatureCalculator()
305 {
306 this->localMode=true;
307 this->projectionPlaneCheck=true;
308 this->sphereRadius=5;
309 this->st=SPHERE_SEARCH;
310 this->nt=AVERAGE;
311 this->montecarlo=false;
312 this->montecarloN=0;
313 this->kRing=3;
314 this->curvatureComputed=false;
315 this->expStep=true;
316 }
317
init(const Eigen::MatrixXd & V,const Eigen::MatrixXi & F)318 IGL_INLINE void CurvatureCalculator::init(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F)
319 {
320 // Normalize vertices
321 vertices = V;
322
323 // vertices = vertices.array() - vertices.minCoeff();
324 // vertices = vertices.array() / vertices.maxCoeff();
325 // vertices = vertices.array() * (1.0/igl::avg_edge_length(V,F));
326
327 faces = F;
328 igl::adjacency_list(F, vertex_to_vertices);
329 igl::vertex_triangle_adjacency(V, F, vertex_to_faces, vertex_to_faces_index);
330 igl::per_face_normals(V, F, face_normals);
331 igl::per_vertex_normals(V, F, face_normals, vertex_normals);
332 }
333
fitQuadric(const Eigen::Vector3d & v,const std::vector<Eigen::Vector3d> & ref,const std::vector<int> & vv,Quadric * q)334 IGL_INLINE void CurvatureCalculator::fitQuadric(const Eigen::Vector3d& v, const std::vector<Eigen::Vector3d>& ref, const std::vector<int>& vv, Quadric *q)
335 {
336 std::vector<Eigen::Vector3d> points;
337 points.reserve (vv.size());
338
339 for (unsigned int i = 0; i < vv.size(); ++i) {
340
341 Eigen::Vector3d cp = vertices.row(vv[i]);
342
343 // vtang non e` il v tangente!!!
344 Eigen::Vector3d vTang = cp - v;
345
346 double x = vTang.dot(ref[0]);
347 double y = vTang.dot(ref[1]);
348 double z = vTang.dot(ref[2]);
349 points.push_back(Eigen::Vector3d (x,y,z));
350 }
351 if (points.size() < 5)
352 {
353 std::cerr << "ASSERT FAILED! fit function requires at least 5 points: Only " << points.size() << " were given." << std::endl;
354 *q = Quadric(0,0,0,0,0);
355 }
356 else
357 {
358 *q = Quadric::fit (points);
359 }
360 }
361
finalEigenStuff(int i,const std::vector<Eigen::Vector3d> & ref,Quadric & q)362 IGL_INLINE void CurvatureCalculator::finalEigenStuff(int i, const std::vector<Eigen::Vector3d>& ref, Quadric& q)
363 {
364
365 const double a = q.a();
366 const double b = q.b();
367 const double c = q.c();
368 const double d = q.d();
369 const double e = q.e();
370
371 // if (fabs(a) < 10e-8 || fabs(b) < 10e-8)
372 // {
373 // std::cout << "Degenerate quadric: " << i << std::endl;
374 // }
375
376 double E = 1.0 + d*d;
377 double F = d*e;
378 double G = 1.0 + e*e;
379
380 Eigen::Vector3d n = Eigen::Vector3d(-d,-e,1.0).normalized();
381
382 double L = 2.0 * a * n[2];
383 double M = b * n[2];
384 double N = 2 * c * n[2];
385
386
387 // ----------------- Eigen stuff
388 Eigen::Matrix2d m;
389 m << L*G - M*F, M*E-L*F, M*E-L*F, N*E-M*F;
390 m = m / (E*G-F*F);
391 Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> eig(m);
392
393 Eigen::Vector2d c_val = eig.eigenvalues();
394 Eigen::Matrix2d c_vec = eig.eigenvectors();
395
396 // std::cerr << "c_val:" << c_val << std::endl;
397 // std::cerr << "c_vec:" << c_vec << std::endl;
398
399 // std::cerr << "c_vec:" << c_vec(0) << " " << c_vec(1) << std::endl;
400
401 c_val = -c_val;
402
403 Eigen::Vector3d v1, v2;
404 v1[0] = c_vec(0);
405 v1[1] = c_vec(1);
406 v1[2] = 0; //d * v1[0] + e * v1[1];
407
408 v2[0] = c_vec(2);
409 v2[1] = c_vec(3);
410 v2[2] = 0; //d * v2[0] + e * v2[1];
411
412
413 // v1 = v1.normalized();
414 // v2 = v2.normalized();
415
416 Eigen::Vector3d v1global = ref[0] * v1[0] + ref[1] * v1[1] + ref[2] * v1[2];
417 Eigen::Vector3d v2global = ref[0] * v2[0] + ref[1] * v2[1] + ref[2] * v2[2];
418
419 v1global.normalize();
420 v2global.normalize();
421
422 v1global *= c_val(0);
423 v2global *= c_val(1);
424
425 if (c_val[0] > c_val[1])
426 {
427 curv[i]=std::vector<double>(2);
428 curv[i][0]=c_val(1);
429 curv[i][1]=c_val(0);
430 curvDir[i]=std::vector<Eigen::Vector3d>(2);
431 curvDir[i][0]=v2global;
432 curvDir[i][1]=v1global;
433 }
434 else
435 {
436 curv[i]=std::vector<double>(2);
437 curv[i][0]=c_val(0);
438 curv[i][1]=c_val(1);
439 curvDir[i]=std::vector<Eigen::Vector3d>(2);
440 curvDir[i][0]=v1global;
441 curvDir[i][1]=v2global;
442 }
443 // ---- end Eigen stuff
444 }
445
getKRing(const int start,const double r,std::vector<int> & vv)446 IGL_INLINE void CurvatureCalculator::getKRing(const int start, const double r, std::vector<int>&vv)
447 {
448 int bufsize=vertices.rows();
449 vv.reserve(bufsize);
450 std::list<std::pair<int,int> > queue;
451 bool* visited = (bool*)calloc(bufsize,sizeof(bool));
452 queue.push_back(std::pair<int,int>(start,0));
453 visited[start]=true;
454 while (!queue.empty())
455 {
456 int toVisit=queue.front().first;
457 int distance=queue.front().second;
458 queue.pop_front();
459 vv.push_back(toVisit);
460 if (distance<(int)r)
461 {
462 for (unsigned int i=0; i<vertex_to_vertices[toVisit].size(); ++i)
463 {
464 int neighbor=vertex_to_vertices[toVisit][i];
465 if (!visited[neighbor])
466 {
467 queue.push_back(std::pair<int,int> (neighbor,distance+1));
468 visited[neighbor]=true;
469 }
470 }
471 }
472 }
473 free(visited);
474 return;
475 }
476
477
getSphere(const int start,const double r,std::vector<int> & vv,int min)478 IGL_INLINE void CurvatureCalculator::getSphere(const int start, const double r, std::vector<int> &vv, int min)
479 {
480 int bufsize=vertices.rows();
481 vv.reserve(bufsize);
482 std::list<int>* queue= new std::list<int>();
483 bool* visited = (bool*)calloc(bufsize,sizeof(bool));
484 queue->push_back(start);
485 visited[start]=true;
486 Eigen::Vector3d me=vertices.row(start);
487 std::priority_queue<std::pair<int, double>, std::vector<std::pair<int, double> >, comparer >* extra_candidates= new std::priority_queue<std::pair<int, double>, std::vector<std::pair<int, double> >, comparer >();
488 while (!queue->empty())
489 {
490 int toVisit=queue->front();
491 queue->pop_front();
492 vv.push_back(toVisit);
493 for (unsigned int i=0; i<vertex_to_vertices[toVisit].size(); ++i)
494 {
495 int neighbor=vertex_to_vertices[toVisit][i];
496 if (!visited[neighbor])
497 {
498 Eigen::Vector3d neigh=vertices.row(neighbor);
499 double distance=(me-neigh).norm();
500 if (distance<r)
501 queue->push_back(neighbor);
502 else if ((int)vv.size()<min)
503 extra_candidates->push(std::pair<int,double>(neighbor,distance));
504 visited[neighbor]=true;
505 }
506 }
507 }
508 while (!extra_candidates->empty() && (int)vv.size()<min)
509 {
510 std::pair<int, double> cand=extra_candidates->top();
511 extra_candidates->pop();
512 vv.push_back(cand.first);
513 for (unsigned int i=0; i<vertex_to_vertices[cand.first].size(); ++i)
514 {
515 int neighbor=vertex_to_vertices[cand.first][i];
516 if (!visited[neighbor])
517 {
518 Eigen::Vector3d neigh=vertices.row(neighbor);
519 double distance=(me-neigh).norm();
520 extra_candidates->push(std::pair<int,double>(neighbor,distance));
521 visited[neighbor]=true;
522 }
523 }
524 }
525 free(extra_candidates);
526 free(queue);
527 free(visited);
528 }
529
project(const Eigen::Vector3d & v,const Eigen::Vector3d & vp,const Eigen::Vector3d & ppn)530 IGL_INLINE Eigen::Vector3d CurvatureCalculator::project(const Eigen::Vector3d& v, const Eigen::Vector3d& vp, const Eigen::Vector3d& ppn)
531 {
532 return (vp - (ppn * ((vp - v).dot(ppn))));
533 }
534
computeReferenceFrame(int i,const Eigen::Vector3d & normal,std::vector<Eigen::Vector3d> & ref)535 IGL_INLINE void CurvatureCalculator::computeReferenceFrame(int i, const Eigen::Vector3d& normal, std::vector<Eigen::Vector3d>& ref )
536 {
537
538 Eigen::Vector3d longest_v=Eigen::Vector3d(vertices.row(vertex_to_vertices[i][0]));
539
540 longest_v=(project(vertices.row(i),longest_v,normal)-Eigen::Vector3d(vertices.row(i))).normalized();
541
542 /* L'ultimo asse si ottiene come prodotto vettoriale tra i due
543 * calcolati */
544 Eigen::Vector3d y_axis=(normal.cross(longest_v)).normalized();
545 ref[0]=longest_v;
546 ref[1]=y_axis;
547 ref[2]=normal;
548 }
549
getAverageNormal(int j,const std::vector<int> & vv,Eigen::Vector3d & normal)550 IGL_INLINE void CurvatureCalculator::getAverageNormal(int j, const std::vector<int>& vv, Eigen::Vector3d& normal)
551 {
552 normal=(vertex_normals.row(j)).normalized();
553 if (localMode)
554 return;
555
556 for (unsigned int i=0; i<vv.size(); ++i)
557 {
558 normal+=vertex_normals.row(vv[i]).normalized();
559 }
560 normal.normalize();
561 }
562
getProjPlane(int j,const std::vector<int> & vv,Eigen::Vector3d & ppn)563 IGL_INLINE void CurvatureCalculator::getProjPlane(int j, const std::vector<int>& vv, Eigen::Vector3d& ppn)
564 {
565 int nr;
566 double a, b, c;
567 double nx, ny, nz;
568 double abcq;
569
570 a = b = c = 0;
571
572 if (localMode)
573 {
574 for (unsigned int i=0; i<vertex_to_faces.at(j).size(); ++i)
575 {
576 Eigen::Vector3d faceNormal=face_normals.row(vertex_to_faces.at(j).at(i));
577 a += faceNormal[0];
578 b += faceNormal[1];
579 c += faceNormal[2];
580 }
581 }
582 else
583 {
584 for (unsigned int i=0; i<vv.size(); ++i)
585 {
586 a+= vertex_normals.row(vv[i])[0];
587 b+= vertex_normals.row(vv[i])[1];
588 c+= vertex_normals.row(vv[i])[2];
589 }
590 }
591 nr = rotateForward (&a, &b, &c);
592 abcq = a*a + b*b + c*c;
593 nx = sqrt (a*a / abcq);
594 ny = sqrt (b*b / abcq);
595 nz = sqrt (1 - nx*nx - ny*ny);
596 rotateBackward (nr, &a, &b, &c);
597 rotateBackward (nr, &nx, &ny, &nz);
598
599 ppn = chooseMax (Eigen::Vector3d(nx, ny, nz), Eigen::Vector3d (a, b, c), a * b);
600 ppn.normalize();
601 }
602
603
getAverageEdge()604 IGL_INLINE double CurvatureCalculator::getAverageEdge()
605 {
606 double sum = 0;
607 int count = 0;
608
609 for (int i = 0; i<faces.rows(); ++i)
610 {
611 for (short unsigned j=0; j<3; ++j)
612 {
613 Eigen::Vector3d p1=vertices.row(faces.row(i)[j]);
614 Eigen::Vector3d p2=vertices.row(faces.row(i)[(j+1)%3]);
615
616 double l = (p1-p2).norm();
617
618 sum+=l;
619 ++count;
620 }
621 }
622
623 return (sum/(double)count);
624 }
625
626
applyProjOnPlane(const Eigen::Vector3d & ppn,const std::vector<int> & vin,std::vector<int> & vout)627 IGL_INLINE void CurvatureCalculator::applyProjOnPlane(const Eigen::Vector3d& ppn, const std::vector<int>& vin, std::vector<int> &vout)
628 {
629 for (std::vector<int>::const_iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
630 if (vertex_normals.row(*vpi) * ppn > 0.0)
631 vout.push_back(*vpi);
632 }
633
applyMontecarlo(const std::vector<int> & vin,std::vector<int> * vout)634 IGL_INLINE void CurvatureCalculator::applyMontecarlo(const std::vector<int>& vin, std::vector<int> *vout)
635 {
636 if (montecarloN >= vin.size ())
637 {
638 *vout = vin;
639 return;
640 }
641
642 float p = ((float) montecarloN) / (float) vin.size();
643 for (std::vector<int>::const_iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
644 {
645 float r;
646 if ((r = ((float)rand () / RAND_MAX)) < p)
647 {
648 vout->push_back(*vpi);
649 }
650 }
651 }
652
computeCurvature()653 IGL_INLINE void CurvatureCalculator::computeCurvature()
654 {
655 //CHECK che esista la mesh
656 const size_t vertices_count=vertices.rows();
657
658 if (vertices_count ==0)
659 return;
660
661 curvDir=std::vector< std::vector<Eigen::Vector3d> >(vertices_count);
662 curv=std::vector<std::vector<double> >(vertices_count);
663
664
665
666 scaledRadius=getAverageEdge()*sphereRadius;
667
668 std::vector<int> vv;
669 std::vector<int> vvtmp;
670 Eigen::Vector3d normal;
671
672 //double time_spent;
673 //double searchtime=0, ref_time=0, fit_time=0, final_time=0;
674
675 for (size_t i=0; i<vertices_count; ++i)
676 {
677 vv.clear();
678 vvtmp.clear();
679 Eigen::Vector3d me=vertices.row(i);
680 switch (st)
681 {
682 case SPHERE_SEARCH:
683 getSphere(i,scaledRadius,vv,6);
684 break;
685 case K_RING_SEARCH:
686 getKRing(i,kRing,vv);
687 break;
688 default:
689 fprintf(stderr,"Error: search type not recognized");
690 return;
691 }
692
693 if (vv.size()<6)
694 {
695 std::cerr << "Could not compute curvature of radius " << scaledRadius << std::endl;
696 return;
697 }
698
699
700 if (projectionPlaneCheck)
701 {
702 vvtmp.reserve (vv.size ());
703 applyProjOnPlane (vertex_normals.row(i), vv, vvtmp);
704 if (vvtmp.size() >= 6 && vvtmp.size()<vv.size())
705 vv = vvtmp;
706 }
707
708
709 switch (nt)
710 {
711 case AVERAGE:
712 getAverageNormal(i,vv,normal);
713 break;
714 case PROJ_PLANE:
715 getProjPlane(i,vv,normal);
716 break;
717 default:
718 fprintf(stderr,"Error: normal type not recognized");
719 return;
720 }
721 if (vv.size()<6)
722 {
723 std::cerr << "Could not compute curvature of radius " << scaledRadius << std::endl;
724 return;
725 }
726 if (montecarlo)
727 {
728 if(montecarloN<6)
729 break;
730 vvtmp.reserve(vv.size());
731 applyMontecarlo(vv,&vvtmp);
732 vv=vvtmp;
733 }
734
735 if (vv.size()<6)
736 return;
737 std::vector<Eigen::Vector3d> ref(3);
738 computeReferenceFrame(i,normal,ref);
739
740 Quadric q;
741 fitQuadric (me, ref, vv, &q);
742 finalEigenStuff(i,ref,q);
743 }
744
745 lastRadius=sphereRadius;
746 curvatureComputed=true;
747 }
748
printCurvature(const std::string & outpath)749 IGL_INLINE void CurvatureCalculator::printCurvature(const std::string& outpath)
750 {
751 using namespace std;
752 if (!curvatureComputed)
753 return;
754
755 std::ofstream of;
756 of.open(outpath.c_str());
757
758 if (!of)
759 {
760 fprintf(stderr, "Error: could not open output file %s\n", outpath.c_str());
761 return;
762 }
763
764 int vertices_count=vertices.rows();
765 of << vertices_count << endl;
766 for (int i=0; i<vertices_count; ++i)
767 {
768 of << curv[i][0] << " " << curv[i][1] << " " << curvDir[i][0][0] << " " << curvDir[i][0][1] << " " << curvDir[i][0][2] << " " <<
769 curvDir[i][1][0] << " " << curvDir[i][1][1] << " " << curvDir[i][1][2] << endl;
770 }
771
772 of.close();
773
774 }
775
776 template <
777 typename DerivedV,
778 typename DerivedF,
779 typename DerivedPD1,
780 typename DerivedPD2,
781 typename DerivedPV1,
782 typename DerivedPV2>
principal_curvature(const Eigen::PlainObjectBase<DerivedV> & V,const Eigen::PlainObjectBase<DerivedF> & F,Eigen::PlainObjectBase<DerivedPD1> & PD1,Eigen::PlainObjectBase<DerivedPD2> & PD2,Eigen::PlainObjectBase<DerivedPV1> & PV1,Eigen::PlainObjectBase<DerivedPV2> & PV2,unsigned radius,bool useKring)783 IGL_INLINE void igl::principal_curvature(
784 const Eigen::PlainObjectBase<DerivedV>& V,
785 const Eigen::PlainObjectBase<DerivedF>& F,
786 Eigen::PlainObjectBase<DerivedPD1>& PD1,
787 Eigen::PlainObjectBase<DerivedPD2>& PD2,
788 Eigen::PlainObjectBase<DerivedPV1>& PV1,
789 Eigen::PlainObjectBase<DerivedPV2>& PV2,
790 unsigned radius,
791 bool useKring)
792 {
793 if (radius < 2)
794 {
795 radius = 2;
796 std::cout << "WARNING: igl::principal_curvature needs a radius >= 2, fixing it to 2." << std::endl;
797 }
798
799 // Preallocate memory
800 PD1.resize(V.rows(),3);
801 PD2.resize(V.rows(),3);
802
803 // Preallocate memory
804 PV1.resize(V.rows(),1);
805 PV2.resize(V.rows(),1);
806
807 // Precomputation
808 CurvatureCalculator cc;
809 cc.init(V.template cast<double>(),F.template cast<int>());
810 cc.sphereRadius = radius;
811
812 if (useKring)
813 {
814 cc.kRing = radius;
815 cc.st = K_RING_SEARCH;
816 }
817
818 // Compute
819 cc.computeCurvature();
820
821 // Copy it back
822 for (unsigned i=0; i<V.rows(); ++i)
823 {
824 PD1.row(i) << cc.curvDir[i][0][0], cc.curvDir[i][0][1], cc.curvDir[i][0][2];
825 PD2.row(i) << cc.curvDir[i][1][0], cc.curvDir[i][1][1], cc.curvDir[i][1][2];
826 PD1.row(i).normalize();
827 PD2.row(i).normalize();
828
829 if (std::isnan(PD1(i,0)) || std::isnan(PD1(i,1)) || std::isnan(PD1(i,2)) || std::isnan(PD2(i,0)) || std::isnan(PD2(i,1)) || std::isnan(PD2(i,2)))
830 {
831 PD1.row(i) << 0,0,0;
832 PD2.row(i) << 0,0,0;
833 }
834
835 PV1(i) = cc.curv[i][0];
836 PV2(i) = cc.curv[i][1];
837
838 if (PD1.row(i) * PD2.row(i).transpose() > 10e-6)
839 {
840 std::cerr << "PRINCIPAL_CURVATURE: Something is wrong with vertex: " << i << std::endl;
841 PD1.row(i) *= 0;
842 PD2.row(i) *= 0;
843 }
844 }
845
846 }
847
848 #ifdef IGL_STATIC_LIBRARY
849 // Explicit template instantiation
850 // generated by autoexplicit.sh
851 template void igl::principal_curvature<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, unsigned int, bool);
852 template void igl::principal_curvature<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(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> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, unsigned int, bool);
853 template void igl::principal_curvature<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, unsigned int, bool);
854 #endif
855