diff --git a/include/Algo/Geometry/curvature.hpp b/include/Algo/Geometry/curvature.hpp index 386e2a7f024a3c8c44708e136cffc177c2f34f62..4fd71c8b3076fd925e9c3980ab15462046ceee06 100644 --- a/include/Algo/Geometry/curvature.hpp +++ b/include/Algo/Geometry/curvature.hpp @@ -317,7 +317,6 @@ void computeCurvatureVertex_NormalCycles( // collect the normal cycle tensor Algo::Selection::Collector_WithinSphere neigh(map, position, radius, thread) ; -// Algo::Selection::Collector_OneRing neigh(map, thread) ; neigh.collectAll(dart) ; MATRIX tensor(0) ; @@ -329,6 +328,16 @@ void computeCurvatureVertex_NormalCycles( const MATRIX& evec = Utils::convertRef(solver.eigenvectors()); normalCycles_SortAndSetEigenComponents(ev,evec,kmax[dart],kmin[dart],Kmax[dart],Kmin[dart],Knormal[dart],normal[dart],thread); + +// if (dart.index % 15000 == 0) +// { +// CGoGNout << solver.eigenvalues() << CGoGNendl; +// CGoGNout << solver.eigenvectors() << CGoGNendl; +// normalCycles_SortTensor(tensor); +// solver.compute(Utils::convertRef(tensor)); +// CGoGNout << solver.eigenvalues() << CGoGNendl; +// CGoGNout << solver.eigenvectors() << CGoGNendl; +// } } template @@ -398,12 +407,12 @@ void normalCycles_SortAndSetEigenComponents( unsigned int thread=0) { // sort eigen components : ev[inormal] has minimal absolute value ; kmin = ev[imin] <= ev[imax] = kmax - int inormal=0, imin, imax, tmp ; + int inormal=0, imin, imax ; if (abs(e_val[1]) < abs(e_val[inormal])) inormal = 1; if (abs(e_val[2]) < abs(e_val[inormal])) inormal = 2; imin = (inormal + 1) % 3; imax = (inormal + 2) % 3; - if (e_val[imax] < e_val[imin]) { tmp = imin ; imin = imax ; imax = tmp ; } + if (e_val[imax] < e_val[imin]) { int tmp = imin ; imin = imax ; imax = tmp ; } // set curvatures from sorted eigen components // warning : Kmin and Kmax are switched w.r.t. kmin and kmax @@ -424,7 +433,6 @@ void normalCycles_SortAndSetEigenComponents( Kmax[2] = e_vec(2,imin); } -/* template void normalCycles_SortTensor( Geom::Matrix<3,3,typename PFP::REAL> & tensor, unsigned int thread=0) { @@ -436,13 +444,27 @@ void normalCycles_SortTensor( Geom::Matrix<3,3,typename PFP::REAL> & tensor, uns // compute eigen components Eigen::SelfAdjointEigenSolver solver (Utils::convertRef(tensor)); - const VEC3& ev = Utils::convertRef(solver.eigenvalues()); - const MATRIX& evec = Utils::convertRef(solver.eigenvectors()); - - // switch kmin and kmax + const VEC3& e_val = Utils::convertRef(solver.eigenvalues()); + const MATRIX& e_vec = Utils::convertRef(solver.eigenvectors()); + // switch kmin and kmax w.r.t. Kmin and Kmax + int inormal=0, imin, imax ; + if (abs(e_val[1]) < abs(e_val[inormal])) inormal = 1; + if (abs(e_val[2]) < abs(e_val[inormal])) inormal = 2; + imin = (inormal + 1) % 3; + imax = (inormal + 2) % 3; + if (e_val[imax] < e_val[imin]) { int tmp = imin ; imin = imax ; imax = tmp ; } + + tensor = e_vec; + int i; REAL v; + i = inormal; v = e_val[inormal]; + tensor(0,i) *= v; tensor(1,i) *= v; tensor(2,i) *= v; + i = imin; v = e_val[imax]; + tensor(0,i) *= v; tensor(1,i) *= v; tensor(2,i) *= v; + i = imax; v = e_val[imin]; + tensor(0,i) *= v; tensor(1,i) *= v; tensor(2,i) *= v; + tensor = tensor*e_vec.transposed(); } -*/ namespace Parallel {