curvature.hpp 26.6 KB
Newer Older
Pierre Kraemer's avatar
Pierre Kraemer committed
1 2 3
/*******************************************************************************
 * CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
 * version 0.1                                                                  *
4
 * Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg           *
Pierre Kraemer's avatar
Pierre Kraemer committed
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
 *                                                                              *
 * This library is free software; you can redistribute it and/or modify it      *
 * under the terms of the GNU Lesser General Public License as published by the *
 * Free Software Foundation; either version 2.1 of the License, or (at your     *
 * option) any later version.                                                   *
 *                                                                              *
 * This library is distributed in the hope that it will be useful, but WITHOUT  *
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or        *
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License  *
 * for more details.                                                            *
 *                                                                              *
 * You should have received a copy of the GNU Lesser General Public License     *
 * along with this library; if not, write to the Free Software Foundation,      *
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.           *
 *                                                                              *
20
 * Web site: http://cgogn.unistra.fr/                                           *
Pierre Kraemer's avatar
Pierre Kraemer committed
21 22 23 24 25 26
 * Contact information: cgogn@unistra.fr                                        *
 *                                                                              *
 *******************************************************************************/

#include "Algo/Geometry/localFrame.h"
#include "Geometry/matrix.h"
27 28
#include "Topology/generic/traversor/traversorCell.h"
#include "Topology/generic/traversor/traversor2.h"
Sylvain Thery's avatar
Sylvain Thery committed
29

Pierre Kraemer's avatar
Pierre Kraemer committed
30 31 32 33 34 35
namespace CGoGN
{

namespace Algo
{

36 37 38
namespace Surface
{

Pierre Kraemer's avatar
Pierre Kraemer committed
39 40 41 42
namespace Geometry
{

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
43 44
void computeCurvatureVertices_QuadraticFitting(
	typename PFP::MAP& map,
45 46 47 48 49
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
Sylvain Thery's avatar
Sylvain Thery committed
50
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin)
Pierre Kraemer's avatar
Pierre Kraemer committed
51
{
52

Sylvain Thery's avatar
Sylvain Thery committed
53
	if (CGoGN::Parallel::NumberOfThreads > 1)
54 55 56 57 58 59 60
	{
		Parallel::computeCurvatureVertices_QuadraticFitting<PFP>(map, position, normal, kmax, kmin, Kmax, Kmin);
		return;
	}

	foreach_cell<VERTEX>(map, [&] (Vertex v)
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
61
		computeCurvatureVertex_QuadraticFitting<PFP>(map, v, position, normal, kmax, kmin, Kmax, Kmin) ;
62
	}
Sylvain Thery's avatar
Sylvain Thery committed
63
	, FORCE_CELL_MARKING);
64 65 66 67

//	TraversorV<typename PFP::MAP> t(map) ;
//	for(Vertex v = t.begin(); v != t.end(); v = t.next())
//		computeCurvatureVertex_QuadraticFitting<PFP>(map, v, position, normal, kmax, kmin, Kmax, Kmin) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
68 69 70
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
71 72
void computeCurvatureVertex_QuadraticFitting(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
73
	Vertex v,
74 75 76 77 78 79
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin)
Pierre Kraemer's avatar
Pierre Kraemer committed
80 81 82 83 84
{
	typedef typename PFP::REAL REAL ;
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::MATRIX33 MATRIX33 ;

Pierre Kraemer's avatar
Pierre Kraemer committed
85
	VEC3 n = normal[v] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
86

Pierre Kraemer's avatar
Pierre Kraemer committed
87
	MATRIX33 localFrame = Algo::Geometry::vertexLocalFrame<PFP>(map, v, position, n) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
88 89 90 91 92
	MATRIX33 invLocalFrame ;
	localFrame.invert(invLocalFrame) ;

	REAL a, b, c, d, e;
	//vertexCubicFitting(map,dart,localFrame,a,b,c,d,e,f,g,h,i) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
93
	vertexQuadraticFitting<PFP>(map, v, localFrame, position, normal, a, b, c, d, e) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
94

Pierre Kraemer's avatar
Pierre Kraemer committed
95
//	REAL kmax_v, kmin_v, Kmax_x, Kmax_y ;
Pierre Kraemer's avatar
Pierre Kraemer committed
96
//	/*int res = */slaev2_(&a, &b, &c, &kmax_v, &kmin_v, &Kmax_x, &Kmax_y) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
97

Pierre Kraemer's avatar
Pierre Kraemer committed
98
	Eigen::Matrix<REAL,2,2> m;
99
	m << 2*a, b, b, 2*c;
Pierre Kraemer's avatar
Pierre Kraemer committed
100 101 102 103 104

	// solve eigen problem
	Eigen::SelfAdjointEigenSolver<Eigen::Matrix<REAL,2,2> > solver(m);

	const Eigen::Matrix<REAL,2,1>& ev = solver.eigenvalues();
Pierre Kraemer's avatar
Pierre Kraemer committed
105 106
	REAL kmax_v = ev[0];
	REAL kmin_v = ev[1];
Pierre Kraemer's avatar
Pierre Kraemer committed
107 108 109

	const Eigen::Matrix<REAL,2,2>& evec = solver.eigenvectors();
	VEC3 Kmax_v(evec(0,0), evec(1,0), 0.0f) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
110
	Kmax_v = invLocalFrame * Kmax_v ;
Pierre Kraemer's avatar
Pierre Kraemer committed
111 112 113
	VEC3 Kmin_v(evec(0,1), evec(1,1), 0.0f) ;
	Kmin_v = invLocalFrame * Kmin_v ;
//	VEC3 Kmin_v = n ^ Kmax_v ;
Pierre Kraemer's avatar
Pierre Kraemer committed
114

Pierre Kraemer's avatar
Pierre Kraemer committed
115
	if (kmax_v < kmin_v)
116
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
117 118 119 120
		kmax[v] = -kmax_v ;
		kmin[v] = -kmin_v ;
		Kmax[v] = Kmax_v ;
		Kmin[v] = Kmin_v ;
121 122 123
	}
	else
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
124 125 126 127
		kmax[v] = -kmin_v ;
		kmin[v] = -kmax_v ;
		Kmax[v] = Kmin_v ;
		Kmin[v] = Kmax_v ;
128
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
129 130 131
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
132
void vertexQuadraticFitting(
Pierre Kraemer's avatar
Pierre Kraemer committed
133
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
134
	Vertex v,
Pierre Kraemer's avatar
Pierre Kraemer committed
135
	typename PFP::MATRIX33& localFrame,
136 137
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
138
    typename PFP::REAL& a, typename PFP::REAL& b, typename PFP::REAL& c, typename PFP::REAL& d, typename PFP::REAL& e)
Pierre Kraemer's avatar
Pierre Kraemer committed
139
{
Pierre Kraemer's avatar
Pierre Kraemer committed
140
	typename PFP::VEC3 p = position[v] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
141

142 143 144 145 146 147
	NLContext nlContext = nlNewContext() ;
	nlMakeCurrent(nlContext) ;
	nlSolverParameteri(NL_NB_VARIABLES, 5) ;
	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE) ;
	nlBegin(NL_SYSTEM) ;
	nlBegin(NL_MATRIX) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
148 149
	foreach_adjacent2<EDGE>(map, v, [&] (Vertex it) {
		typename PFP::VEC3 itp = position[it] ;
150
        quadraticFittingAddVertexPos<PFP>(p, itp, localFrame) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
151 152 153
		typename PFP::VEC3 itn = normal[it] ;
		quadraticFittingAddVertexNormal<PFP>(itp, itn, p, localFrame) ;
	});
154 155 156 157 158 159 160 161 162 163 164
	nlEnd(NL_MATRIX) ;
	nlEnd(NL_SYSTEM) ;
	nlSolve() ;

	a = nlGetVariable(0) ;
	b = nlGetVariable(1) ;
	c = nlGetVariable(2) ;
	d = nlGetVariable(3) ;
	e = nlGetVariable(4) ;

	nlDeleteContext(nlContext) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
165 166 167
}

template <typename PFP>
168
void quadraticFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame)
Pierre Kraemer's avatar
Pierre Kraemer committed
169 170 171 172
{
	typename PFP::VEC3 vec = v - p ;
	vec = localFrame * vec ;

173 174 175 176 177 178 179 180
	nlRowParameterd(NL_RIGHT_HAND_SIDE, vec[2]) ;
	nlBegin(NL_ROW) ;
	nlCoefficient(0, vec[0]*vec[0]) ;
	nlCoefficient(1, vec[0]*vec[1]) ;
	nlCoefficient(2, vec[1]*vec[1]) ;
	nlCoefficient(3, vec[0]) ;
	nlCoefficient(4, vec[1]) ;
	nlEnd(NL_ROW) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
181 182 183
}

template <typename PFP>
184
void quadraticFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame)
Pierre Kraemer's avatar
Pierre Kraemer committed
185 186 187 188 189
{
	typename PFP::VEC3 vec = v - p ;
	vec = localFrame * vec ;
	typename PFP::VEC3 norm = localFrame * n ;

190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206
	nlRowParameterd(NL_RIGHT_HAND_SIDE, -1.0f * norm[0]) ;
	nlBegin(NL_ROW);
	nlCoefficient(0, 2.0f * vec[0] * norm[2]) ;
	nlCoefficient(1, vec[1] * norm[2]) ;
	nlCoefficient(2, 0) ;
	nlCoefficient(3, 1.0f * norm[2]) ;
	nlCoefficient(4, 0) ;
	nlEnd(NL_ROW) ;

	nlRowParameterd(NL_RIGHT_HAND_SIDE, -1.0f * norm[1]) ;
	nlBegin(NL_ROW);
	nlCoefficient(0, 0) ;
	nlCoefficient(1, vec[0] * norm[2]) ;
	nlCoefficient(2, 2.0f * vec[1] * norm[2]) ;
	nlCoefficient(3, 0) ;
	nlCoefficient(4, 1.0f * norm[2]) ;
	nlEnd(NL_ROW) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
207 208 209
}
/*
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
210
void vertexCubicFitting(Dart dart, gmtl::Vec3f& normal, float& a, float& b, float& c, float& d, float& e, float& f, float& g, float& h, float& i)
Pierre Kraemer's avatar
Pierre Kraemer committed
211 212
{
	gmtl::Matrix33f localFrame, invLocalFrame ;
213
	Geometry::vertexLocalFrame<PFP>(m_map,dart,normal,localFrame) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
214 215 216 217 218
	gmtl::invertFull(invLocalFrame, localFrame) ;
	gmtl::Vec3f p = m_map.getVertexEmb(dart)->getPosition() ;
	solverC->reset(false) ;
	solverC->set_least_squares(true) ;
	solverC->begin_system() ;
219 220
	Traversor2VVaE<typename PFP::MAP> tav(map, dart) ;
	for(Dart it = tav.begin(); it != tav.end(); it = tav.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
221 222
	{
		// 1-ring vertices
223
		gmtl::Vec3f v = m_map.getVertexEmb(it)->getPosition() ;
Pierre Kraemer's avatar
Pierre Kraemer committed
224
		cubicFittingAddVertexPos(v,p,localFrame) ;
225
		gmtl::Vec3f n = m_normalsV[m_map.getVertexEmb(it)->getLabel()] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
226
		cubicFittingAddVertexNormal(v,n,p,localFrame) ;
227
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
228 229
	solverC->end_system() ;
	solverC->solve() ;
230

Pierre Kraemer's avatar
Pierre Kraemer committed
231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
	a = solverC->variable(0).value() ;
	b = solverC->variable(1).value() ;
	c = solverC->variable(2).value() ;
	d = solverC->variable(3).value() ;
	e = solverC->variable(4).value() ;
	f = solverC->variable(5).value() ;
	g = solverC->variable(6).value() ;
	h = solverC->variable(7).value() ;
	i = solverC->variable(8).value() ;

//	normal = gmtl::Vec3f(-h, -i, 1.0f) ;
//	gmtl::normalize(normal) ;
//	normal = invLocalFrame * normal ;
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
247
void cubicFittingAddVertexPos(gmtl::Vec3f& v, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
Pierre Kraemer's avatar
Pierre Kraemer committed
248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
{
	gmtl::Vec3f vec = v - p ;
	vec = localFrame * vec ;
	solverC->begin_row() ;

	solverC->add_coefficient(0, vec[0]*vec[0]*vec[0]) ;
	solverC->add_coefficient(1, vec[0]*vec[0]*vec[1]) ;
	solverC->add_coefficient(2, vec[0]*vec[1]*vec[1]) ;
	solverC->add_coefficient(3, vec[1]*vec[1]*vec[1]) ;
	solverC->add_coefficient(4, vec[0]*vec[0]) ;
	solverC->add_coefficient(5, vec[0]*vec[1]) ;
	solverC->add_coefficient(6, vec[1]*vec[1]) ;
	solverC->add_coefficient(7, vec[0]) ;
	solverC->add_coefficient(8, vec[1]) ;

	solverC->set_right_hand_side(vec[2]) ;
	solverC->end_row() ;
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
268
void cubicFittingAddVertexNormal(gmtl::Vec3f& v, gmtl::Vec3f& n, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
Pierre Kraemer's avatar
Pierre Kraemer committed
269 270 271 272 273
{
	gmtl::Vec3f vec = v - p ;
	vec = localFrame * vec ;
	gmtl::Vec3f norm = localFrame * n ;

274
	solverC->begin_row() ;
Pierre Kraemer's avatar
Pierre Kraemer committed
275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
	solverC->add_coefficient(0, 3.0f*vec[0]*vec[0]) ;
	solverC->add_coefficient(1, 2.0f*vec[0]*vec[1]) ;
	solverC->add_coefficient(2, vec[1]*vec[1]) ;
	solverC->add_coefficient(3, 0) ;
	solverC->add_coefficient(4, 2.0f*vec[0]) ;
	solverC->add_coefficient(5, vec[1]) ;
	solverC->add_coefficient(6, 0) ;
	solverC->add_coefficient(7, 1.0f) ;
	solverC->add_coefficient(8, 0) ;
	solverC->set_right_hand_side(-1.0f*norm[0]/norm[2]) ;
	solverC->end_row() ;

	solverC->begin_row() ;
	solverC->add_coefficient(0, 0) ;
	solverC->add_coefficient(1, vec[0]*vec[0]) ;
	solverC->add_coefficient(2, 2.0f*vec[0]*vec[1]) ;
	solverC->add_coefficient(3, 3.0f*vec[1]*vec[1]) ;
	solverC->add_coefficient(4, 0) ;
	solverC->add_coefficient(5, vec[0]) ;
	solverC->add_coefficient(6, 2.0f*vec[1]) ;
	solverC->add_coefficient(7, 0) ;
	solverC->add_coefficient(8, 1.0f) ;
	solverC->set_right_hand_side(-1.0f*norm[1]/norm[2]) ;
	solverC->end_row() ;
}
*/
Pierre Kraemer's avatar
Pierre Kraemer committed
301 302 303 304 305

template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
306 307 308 309 310 311 312
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
313
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Pierre Kraemer's avatar
Pierre Kraemer committed
314
{
Sylvain Thery's avatar
Sylvain Thery committed
315
	if (CGoGN::Parallel::NumberOfThreads > 1)
316 317 318 319 320 321 322
	{
		Parallel::computeCurvatureVertices_NormalCycles<PFP>(map, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal);
		return;
	}

	foreach_cell<VERTEX>(map, [&] (Vertex v)
	{
Sylvain Thery's avatar
Sylvain Thery committed
323
		computeCurvatureVertex_NormalCycles<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ;
324
	}
Sylvain Thery's avatar
Sylvain Thery committed
325
	,FORCE_CELL_MARKING);
326 327 328

//	TraversorV<typename PFP::MAP> t(map) ;
//	for(Vertex v = t.begin(); v != t.end(); v = t.next())
Sylvain Thery's avatar
Sylvain Thery committed
329
//		computeCurvatureVertex_NormalCycles<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
330 331 332 333 334
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
335
	Vertex v,
Pierre Kraemer's avatar
Pierre Kraemer committed
336
	typename PFP::REAL radius,
337 338 339 340 341 342 343
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
344
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Pierre Kraemer's avatar
Pierre Kraemer committed
345
{
346 347
	typedef typename PFP::REAL REAL ;
	typedef typename PFP::VEC3 VEC3 ;
348 349
	typedef Geom::Matrix<3,3,REAL> MATRIX;
	typedef Eigen::Matrix<REAL,3,1> E_VEC3;
Sauvage's avatar
Sauvage committed
350
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;
351

352
	// collect the normal cycle tensor
Sylvain Thery's avatar
Sylvain Thery committed
353
	Selection::Collector_WithinSphere<PFP> neigh(map, position, radius) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
354
	neigh.collectAll(v) ;
355

356
	MATRIX tensor(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
357
	neigh.computeNormalCyclesTensor(position, edgeangle, tensor);
358

359
	// solve eigen problem
Pierre Kraemer's avatar
Pierre Kraemer committed
360
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
361 362
	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());
363

Sylvain Thery's avatar
Sylvain Thery committed
364
	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v]);
Sauvage's avatar
Sauvage committed
365 366 367 368 369 370 371 372 373 374

//	if (dart.index % 15000 == 0)
//	{
//		CGoGNout << solver.eigenvalues() << CGoGNendl;
//		CGoGNout << solver.eigenvectors() << CGoGNendl;
//		normalCycles_SortTensor<PFP>(tensor);
//		solver.compute(Utils::convertRef<E_MATRIX>(tensor));
//		CGoGNout << solver.eigenvalues() << CGoGNendl;
//		CGoGNout << solver.eigenvectors() << CGoGNendl;
//	}
375 376
}

Sauvage's avatar
Sauvage committed
377 378 379 380
template <typename PFP>
void computeCurvatureVertices_NormalCycles_Projected(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
381 382 383 384 385 386 387
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
388
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Sauvage's avatar
Sauvage committed
389
{
Sylvain Thery's avatar
Sylvain Thery committed
390
	if (CGoGN::Parallel::NumberOfThreads > 1)
391 392 393 394 395 396 397
	{
		Parallel::computeCurvatureVertices_NormalCycles_Projected<PFP>(map, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal);
		return;
	}

	foreach_cell<VERTEX>(map, [&] (Vertex v)
	{
Sylvain Thery's avatar
Sylvain Thery committed
398
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ;
399
	}
Sylvain Thery's avatar
Sylvain Thery committed
400
	,FORCE_CELL_MARKING);
401 402 403

//	TraversorV<typename PFP::MAP> t(map) ;
//	for(Vertex v = t.begin(); v.dart != t.end(); v = t.next())
Sylvain Thery's avatar
Sylvain Thery committed
404
//		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ;
Sauvage's avatar
Sauvage committed
405 406 407 408 409
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles_Projected(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
410
	Vertex v,
Sauvage's avatar
Sauvage committed
411
	typename PFP::REAL radius,
412 413 414 415 416 417 418
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
419
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Sauvage's avatar
Sauvage committed
420
{
421
	typedef typename PFP::REAL REAL ;
Sauvage's avatar
Sauvage committed
422 423 424 425
	typedef typename PFP::VEC3 VEC3 ;
	typedef Geom::Matrix<3,3,REAL> MATRIX;
	typedef Eigen::Matrix<REAL,3,1> E_VEC3;
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;
426

427
	// collect the normal cycle tensor
Sylvain Thery's avatar
Sylvain Thery committed
428
	Selection::Collector_WithinSphere<PFP> neigh(map, position, radius) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
429
	neigh.collectAll(v) ;
430

Sauvage's avatar
Sauvage committed
431
	MATRIX tensor(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
432
	neigh.computeNormalCyclesTensor(position, edgeangle, tensor);
433

Sauvage's avatar
Sauvage committed
434
	// project the tensor
Sylvain Thery's avatar
Sylvain Thery committed
435
	normalCycles_ProjectTensor<PFP>(tensor, normal[v]);
436

Sauvage's avatar
Sauvage committed
437
	// solve eigen problem
Pierre Kraemer's avatar
Pierre Kraemer committed
438
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
Sauvage's avatar
Sauvage committed
439 440 441
	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());

Sylvain Thery's avatar
Sylvain Thery committed
442
	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v]);
Sauvage's avatar
Sauvage committed
443 444
}

445 446 447
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
448
	Algo::Surface::Selection::Collector<PFP>& neigh,
449 450 451 452 453 454 455
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
456
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
457
{
458
	TraversorV<typename PFP::MAP> t(map) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
459
	for(Vertex v = t.begin(); v != t.end(); v = t.next())
Sylvain Thery's avatar
Sylvain Thery committed
460
		computeCurvatureVertex_NormalCycles<PFP>(map, v, neigh, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ;
461 462 463 464
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles(
Pierre Kraemer's avatar
Pierre Kraemer committed
465 466
	Vertex v,
	Algo::Surface::Selection::Collector<PFP>& neigh,
467 468 469 470 471 472 473
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
474
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
475 476 477 478 479
{
	typedef typename PFP::REAL REAL ;
	typedef typename PFP::VEC3 VEC3 ;
	typedef Geom::Matrix<3,3,REAL> MATRIX;
	typedef Eigen::Matrix<REAL,3,1> E_VEC3;
Sauvage's avatar
Sauvage committed
480
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;
481 482

	// collect the normal cycle tensor
Pierre Kraemer's avatar
Pierre Kraemer committed
483
	neigh.collectAll(v) ;
484

485
	MATRIX tensor(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
486
	neigh.computeNormalCyclesTensor(position, edgeangle, tensor);
487

488
	// solve eigen problem
Pierre Kraemer's avatar
Pierre Kraemer committed
489
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
490 491 492
	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());

Sylvain Thery's avatar
Sylvain Thery committed
493
	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v]);
494 495
}

Sauvage's avatar
Sauvage committed
496 497 498
template <typename PFP>
void computeCurvatureVertices_NormalCycles_Projected(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
499
	Algo::Surface::Selection::Collector<PFP>& neigh,
500 501 502 503 504 505 506
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
507
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Sauvage's avatar
Sauvage committed
508
{
509
	TraversorV<typename PFP::MAP> t(map) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
510
	for(Vertex v = t.begin(); v != t.end(); v = t.next())
Sylvain Thery's avatar
Sylvain Thery committed
511
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, neigh, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ;
Sauvage's avatar
Sauvage committed
512 513 514 515
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles_Projected(
Pierre Kraemer's avatar
Pierre Kraemer committed
516 517
	Vertex v,
	Algo::Surface::Selection::Collector<PFP>& neigh,
518 519 520 521 522 523 524
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
525
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Sauvage's avatar
Sauvage committed
526 527 528 529 530 531 532 533
{
	typedef typename PFP::REAL REAL ;
	typedef typename PFP::VEC3 VEC3 ;
	typedef Geom::Matrix<3,3,REAL> MATRIX;
	typedef Eigen::Matrix<REAL,3,1> E_VEC3;
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;

	// collect the normal cycle tensor
Pierre Kraemer's avatar
Pierre Kraemer committed
534
	neigh.collectAll(v) ;
Sauvage's avatar
Sauvage committed
535 536

	MATRIX tensor(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
537
	neigh.computeNormalCyclesTensor(position, edgeangle, tensor);
Sauvage's avatar
Sauvage committed
538 539

	// project the tensor
Sylvain Thery's avatar
Sylvain Thery committed
540
	normalCycles_ProjectTensor<PFP>(tensor, normal[v]);
Sauvage's avatar
Sauvage committed
541 542

	// solve eigen problem
Pierre Kraemer's avatar
Pierre Kraemer committed
543
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
Sauvage's avatar
Sauvage committed
544 545 546
	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());

Sylvain Thery's avatar
Sylvain Thery committed
547
	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v]);
Sauvage's avatar
Sauvage committed
548
}
549

550 551 552 553 554 555 556 557 558
template <typename PFP>
void normalCycles_SortAndSetEigenComponents(
	const typename PFP::VEC3& e_val,
	const Geom::Matrix<3,3,typename PFP::REAL> & e_vec,
	typename PFP::REAL& kmax,
	typename PFP::REAL& kmin,
	typename PFP::VEC3& Kmax,
	typename PFP::VEC3& Kmin,
	typename PFP::VEC3& Knormal,
Sylvain Thery's avatar
Sylvain Thery committed
559
	const typename PFP::VEC3& normal)
560
{
Sauvage's avatar
Sauvage committed
561
	// sort eigen components : ev[inormal] has minimal absolute value ; kmin = ev[imin] <= ev[imax] = kmax
Pierre Kraemer's avatar
Pierre Kraemer committed
562
	int inormal = 0, imin, imax ;
Pierre Kraemer's avatar
Pierre Kraemer committed
563 564
	if (fabs(e_val[1]) < fabs(e_val[inormal])) inormal = 1;
	if (fabs(e_val[2]) < fabs(e_val[inormal])) inormal = 2;
Sauvage's avatar
Sauvage committed
565 566
	imin = (inormal + 1) % 3;
	imax = (inormal + 2) % 3;
Sauvage's avatar
Sauvage committed
567
	if (e_val[imax] < e_val[imin]) { int tmp = imin ; imin = imax ; imax = tmp ; }
568

569
	// set curvatures from sorted eigen components
570 571
	// warning : Kmin and Kmax are switched w.r.t. kmin and kmax
	// normal direction : minimal absolute eigen value
Sauvage's avatar
Sauvage committed
572 573 574
	Knormal[0] = e_vec(0,inormal);
	Knormal[1] = e_vec(1,inormal);
	Knormal[2] = e_vec(2,inormal);
575
	if (Knormal * normal < 0) Knormal *= -1; // change orientation
576
	// min curvature
Sauvage's avatar
Sauvage committed
577 578 579 580
	kmin = e_val[imin] ;
	Kmin[0] = e_vec(0,imax);
	Kmin[1] = e_vec(1,imax);
	Kmin[2] = e_vec(2,imax);
581
	// max curvature
Sauvage's avatar
Sauvage committed
582 583 584 585
	kmax = e_val[imax] ;
	Kmax[0] = e_vec(0,imin);
	Kmax[1] = e_vec(1,imin);
	Kmax[2] = e_vec(2,imin);
Pierre Kraemer's avatar
Pierre Kraemer committed
586 587
}

Sauvage's avatar
Sauvage committed
588
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
589
void normalCycles_SortTensor(Geom::Matrix<3,3,typename PFP::REAL>& tensor)
Sauvage's avatar
Sauvage committed
590 591 592 593 594 595 596 597
{
	typedef typename PFP::REAL REAL ;
	typedef typename PFP::VEC3 VEC3 ;
	typedef Geom::Matrix<3,3,REAL> MATRIX;
	typedef Eigen::Matrix<REAL,3,1> E_VEC3;
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;

	// compute eigen components
Pierre Kraemer's avatar
Pierre Kraemer committed
598
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
Sauvage's avatar
Sauvage committed
599 600
	const VEC3& e_val = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& e_vec = Utils::convertRef<MATRIX>(solver.eigenvectors());
Sauvage's avatar
Sauvage committed
601

Sauvage's avatar
Sauvage committed
602
	// switch kmin and kmax w.r.t. Kmin and Kmax
Pierre Kraemer's avatar
Pierre Kraemer committed
603
	int inormal = 0, imin, imax ;
Pierre Kraemer's avatar
Pierre Kraemer committed
604 605
	if (fabs(e_val[1]) < fabs(e_val[inormal])) inormal = 1;
	if (fabs(e_val[2]) < fabs(e_val[inormal])) inormal = 2;
Sauvage's avatar
Sauvage committed
606 607 608 609 610 611 612 613 614 615 616 617 618
	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();
Sauvage's avatar
Sauvage committed
619 620
}

Sauvage's avatar
Sauvage committed
621
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
622
void normalCycles_ProjectTensor(Geom::Matrix<3,3,typename PFP::REAL>& tensor, const typename PFP::VEC3& normal_vector)
Sauvage's avatar
Sauvage committed
623 624 625
{
	Geom::Matrix<3,3,typename PFP::REAL> proj;
	proj.identity();
Pierre Kraemer's avatar
Pierre Kraemer committed
626
	proj -= Geom::transposed_vectors_mult(normal_vector, normal_vector);
Sauvage's avatar
Sauvage committed
627
	tensor = proj * tensor * proj;
Pierre Kraemer's avatar
Pierre Kraemer committed
628 629
}

630 631 632

namespace Parallel
{
Pierre Kraemer's avatar
Pierre Kraemer committed
633

634 635 636 637
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
638 639 640 641 642 643 644 645
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
646
{
647
	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
648
	if (!map.template isOrbitEmbedded<VERTEX>())
Pierre Kraemer's avatar
Pierre Kraemer committed
649 650
		Algo::Topo::initAllOrbitsEmbedding<VERTEX>(map);

651
	if (!map.template isOrbitEmbedded<EDGE>())
Pierre Kraemer's avatar
Pierre Kraemer committed
652 653
		Algo::Topo::initAllOrbitsEmbedding<EDGE>(map);

654
	if (!map.template isOrbitEmbedded<FACE>())
Pierre Kraemer's avatar
Pierre Kraemer committed
655
		Algo::Topo::initAllOrbitsEmbedding<FACE>(map);
656

Sylvain Thery's avatar
Sylvain Thery committed
657
	CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
Sylvain Thery's avatar
Sylvain Thery committed
658
	{
Sylvain Thery's avatar
Sylvain Thery committed
659
		computeCurvatureVertex_NormalCycles<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
660
	}, FORCE_CELL_MARKING);
661 662
}

663
template <typename PFP>
664
void computeCurvatureVertices_NormalCycles_Projected(
665
	typename PFP::MAP& map,
666
	typename PFP::REAL radius,
667 668 669 670 671 672 673 674
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
675
{
Pierre Kraemer's avatar
Pierre Kraemer committed
676
//	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
677
	if (!map.template isOrbitEmbedded<VERTEX>())
Pierre Kraemer's avatar
Pierre Kraemer committed
678 679
		Algo::Topo::initAllOrbitsEmbedding<VERTEX>(map);

680
	if (!map.template isOrbitEmbedded<EDGE>())
Pierre Kraemer's avatar
Pierre Kraemer committed
681 682
		Algo::Topo::initAllOrbitsEmbedding<EDGE>(map);

683
	if (!map.template isOrbitEmbedded<FACE>())
Pierre Kraemer's avatar
Pierre Kraemer committed
684
		Algo::Topo::initAllOrbitsEmbedding<FACE>(map);
685

Sylvain Thery's avatar
Sylvain Thery committed
686
	CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
687
	{
Sylvain Thery's avatar
Sylvain Thery committed
688
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
689
	}, FORCE_CELL_MARKING);
690 691 692 693 694
}

template <typename PFP>
void computeCurvatureVertices_QuadraticFitting(
	typename PFP::MAP& map,
695 696 697 698 699 700
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmin)
701
{
Sylvain Thery's avatar
Sylvain Thery committed
702
	CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
Sylvain Thery's avatar
Sylvain Thery committed
703
	{
Sylvain Thery's avatar
Sylvain Thery committed
704
		computeCurvatureVertex_QuadraticFitting<PFP>(map, v, position, normal, kmax, kmin, Kmax, Kmin) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
705
	}, FORCE_CELL_MARKING);
Sylvain Thery's avatar
Sylvain Thery committed
706

707 708
}

709
} // namespace Parallel
710 711


Pierre Kraemer's avatar
Pierre Kraemer committed
712 713
} // namespace Geometry

714
} // namespace Surface
715

Pierre Kraemer's avatar
Pierre Kraemer committed
716 717 718
} // namespace Algo

} // namespace CGoGN