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
{
Sylvain Thery's avatar
Sylvain Thery committed
52

Sylvain Thery's avatar
Sylvain Thery committed
53
	if (CGoGN::Parallel::NumberOfThreads > 1)
Sylvain Thery's avatar
Sylvain Thery committed
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) ;
Sylvain Thery's avatar
Sylvain Thery committed
62
	}
Sylvain Thery's avatar
Sylvain Thery committed
63
	, FORCE_CELL_MARKING);
Sylvain Thery's avatar
Sylvain Thery committed
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
	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,
309
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
310 311 312 313
	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
314
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Pierre Kraemer's avatar
Pierre Kraemer committed
315
{
Sylvain Thery's avatar
Sylvain Thery committed
316
	if (CGoGN::Parallel::NumberOfThreads > 1)
Sylvain Thery's avatar
Sylvain Thery committed
317
	{
318
		Parallel::computeCurvatureVertices_NormalCycles<PFP>(map, radius, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal);
Sylvain Thery's avatar
Sylvain Thery committed
319 320 321 322 323
		return;
	}

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

template <typename PFP>
void computeCurvatureVertex_NormalCycles(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
332
	Vertex v,
Pierre Kraemer's avatar
Pierre Kraemer committed
333
	typename PFP::REAL radius,
334 335 336
	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,
337
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
338 339 340 341
	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
342
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Pierre Kraemer's avatar
Pierre Kraemer committed
343
{
344 345
	typedef typename PFP::REAL REAL ;
	typedef typename PFP::VEC3 VEC3 ;
346 347
	typedef Geom::Matrix<3,3,REAL> MATRIX;
	typedef Eigen::Matrix<REAL,3,1> E_VEC3;
Sauvage's avatar
Sauvage committed
348
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;
349

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

354
	MATRIX tensor(0) ;
355
	neigh.computeNormalCyclesTensor(position, edgeangle, edgearea, tensor);
356

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

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

Sauvage's avatar
Sauvage committed
365 366 367 368
template <typename PFP>
void computeCurvatureVertices_NormalCycles_Projected(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
369 370 371
	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,
372
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
373 374 375 376
	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
377
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Sauvage's avatar
Sauvage committed
378
{
379 380 381 382 383
//	if (CGoGN::Parallel::NumberOfThreads > 1)
//	{
//		Parallel::computeCurvatureVertices_NormalCycles_Projected<PFP>(map, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal);
//		return;
//	}
Sylvain Thery's avatar
Sylvain Thery committed
384 385 386

	foreach_cell<VERTEX>(map, [&] (Vertex v)
	{
387
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, radius, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal) ;
Sylvain Thery's avatar
Sylvain Thery committed
388
	}
Sylvain Thery's avatar
Sylvain Thery committed
389
	,FORCE_CELL_MARKING);
Sauvage's avatar
Sauvage committed
390 391 392 393 394
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles_Projected(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
395
	Vertex v,
Sauvage's avatar
Sauvage committed
396
	typename PFP::REAL radius,
397 398 399
	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,
400
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
401 402 403 404
	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
405
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Sauvage's avatar
Sauvage committed
406
{
407
	typedef typename PFP::REAL REAL ;
Sauvage's avatar
Sauvage committed
408 409 410 411
	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;
412

413
	// collect the normal cycle tensor
Sylvain Thery's avatar
Sylvain Thery committed
414
	Selection::Collector_WithinSphere<PFP> neigh(map, position, radius) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
415
	neigh.collectAll(v) ;
416

Sauvage's avatar
Sauvage committed
417
	MATRIX tensor(0) ;
418
	neigh.computeNormalCyclesTensor(position, edgeangle, edgearea, tensor);
419

Sauvage's avatar
Sauvage committed
420
	// project the tensor
Sylvain Thery's avatar
Sylvain Thery committed
421
	normalCycles_ProjectTensor<PFP>(tensor, normal[v]);
422

Sauvage's avatar
Sauvage committed
423
	// solve eigen problem
Pierre Kraemer's avatar
Pierre Kraemer committed
424
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
Sauvage's avatar
Sauvage committed
425 426 427
	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());

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

431 432 433
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
434
	Algo::Surface::Selection::Collector<PFP>& neigh,
435 436 437
	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,
438
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
439 440 441 442
	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
443
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
444
{
445
	TraversorV<typename PFP::MAP> t(map) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
446
	for(Vertex v = t.begin(); v != t.end(); v = t.next())
447
		computeCurvatureVertex_NormalCycles<PFP>(map, v, neigh, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal) ;
448 449 450 451
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles(
Pierre Kraemer's avatar
Pierre Kraemer committed
452 453
	Vertex v,
	Algo::Surface::Selection::Collector<PFP>& neigh,
454 455 456
	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,
457
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
458 459 460 461
	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
462
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
463 464 465 466 467
{
	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
468
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;
469 470

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

473
	MATRIX tensor(0) ;
474
	neigh.computeNormalCyclesTensor(position, edgeangle, edgearea, tensor);
475

476
	// solve eigen problem
Pierre Kraemer's avatar
Pierre Kraemer committed
477
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
478 479 480
	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());

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

Sauvage's avatar
Sauvage committed
484 485 486
template <typename PFP>
void computeCurvatureVertices_NormalCycles_Projected(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
487
	Algo::Surface::Selection::Collector<PFP>& neigh,
488 489 490
	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,
491
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
492 493 494 495
	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
496
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Sauvage's avatar
Sauvage committed
497
{
498
	TraversorV<typename PFP::MAP> t(map) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
499
	for(Vertex v = t.begin(); v != t.end(); v = t.next())
500
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, neigh, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal) ;
Sauvage's avatar
Sauvage committed
501 502 503 504
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles_Projected(
Pierre Kraemer's avatar
Pierre Kraemer committed
505 506
	Vertex v,
	Algo::Surface::Selection::Collector<PFP>& neigh,
507 508 509
	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,
510
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
511 512 513 514
	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
515
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Knormal)
Sauvage's avatar
Sauvage committed
516 517 518 519 520 521 522 523
{
	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
524
	neigh.collectAll(v) ;
Sauvage's avatar
Sauvage committed
525 526

	MATRIX tensor(0) ;
527
	neigh.computeNormalCyclesTensor(position, edgeangle, edgearea, tensor);
Sauvage's avatar
Sauvage committed
528 529

	// project the tensor
Sylvain Thery's avatar
Sylvain Thery committed
530
	normalCycles_ProjectTensor<PFP>(tensor, normal[v]);
Sauvage's avatar
Sauvage committed
531 532

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

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

540 541 542 543 544 545 546 547 548
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
549
	const typename PFP::VEC3& normal)
550
{
Sauvage's avatar
Sauvage committed
551
	// sort eigen components : ev[inormal] has minimal absolute value ; kmin = ev[imin] <= ev[imax] = kmax
Pierre Kraemer's avatar
Pierre Kraemer committed
552
	int inormal = 0, imin, imax ;
Pierre Kraemer's avatar
Pierre Kraemer committed
553 554
	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
555 556
	imin = (inormal + 1) % 3;
	imax = (inormal + 2) % 3;
Sauvage's avatar
Sauvage committed
557
	if (e_val[imax] < e_val[imin]) { int tmp = imin ; imin = imax ; imax = tmp ; }
558

559
	// set curvatures from sorted eigen components
560 561
	// warning : Kmin and Kmax are switched w.r.t. kmin and kmax
	// normal direction : minimal absolute eigen value
Sauvage's avatar
Sauvage committed
562 563 564
	Knormal[0] = e_vec(0,inormal);
	Knormal[1] = e_vec(1,inormal);
	Knormal[2] = e_vec(2,inormal);
565
	if (Knormal * normal < 0) Knormal *= -1; // change orientation
566
	// min curvature
Sauvage's avatar
Sauvage committed
567 568 569 570
	kmin = e_val[imin] ;
	Kmin[0] = e_vec(0,imax);
	Kmin[1] = e_vec(1,imax);
	Kmin[2] = e_vec(2,imax);
571
	// max curvature
Sauvage's avatar
Sauvage committed
572 573 574 575
	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
576 577
}

Sauvage's avatar
Sauvage committed
578
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
579
void normalCycles_SortTensor(Geom::Matrix<3,3,typename PFP::REAL>& tensor)
Sauvage's avatar
Sauvage committed
580 581 582 583 584 585 586 587
{
	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
588
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
Sauvage's avatar
Sauvage committed
589 590
	const VEC3& e_val = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& e_vec = Utils::convertRef<MATRIX>(solver.eigenvectors());
Sauvage's avatar
Sauvage committed
591

Sauvage's avatar
Sauvage committed
592
	// switch kmin and kmax w.r.t. Kmin and Kmax
Pierre Kraemer's avatar
Pierre Kraemer committed
593
	int inormal = 0, imin, imax ;
Pierre Kraemer's avatar
Pierre Kraemer committed
594 595
	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
596 597 598 599 600 601 602 603 604 605 606 607 608
	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
609 610
}

Sauvage's avatar
Sauvage committed
611
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
612
void normalCycles_ProjectTensor(Geom::Matrix<3,3,typename PFP::REAL>& tensor, const typename PFP::VEC3& normal_vector)
Sauvage's avatar
Sauvage committed
613 614 615
{
	Geom::Matrix<3,3,typename PFP::REAL> proj;
	proj.identity();
Pierre Kraemer's avatar
Pierre Kraemer committed
616
	proj -= Geom::transposed_vectors_mult(normal_vector, normal_vector);
Sauvage's avatar
Sauvage committed
617
	tensor = proj * tensor * proj;
Pierre Kraemer's avatar
Pierre Kraemer committed
618 619
}

620 621 622

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

624 625 626 627
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
628 629 630
	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,
631
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
632 633 634 635 636
	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)
637
{
Sylvain Thery's avatar
Sylvain Thery committed
638
	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
Pierre Kraemer's avatar
Pierre Kraemer committed
639
	if (!map.template isOrbitEmbedded<VERTEX>())
Pierre Kraemer's avatar
Pierre Kraemer committed
640 641
		Algo::Topo::initAllOrbitsEmbedding<VERTEX>(map);

Pierre Kraemer's avatar
Pierre Kraemer committed
642
	if (!map.template isOrbitEmbedded<EDGE>())
Pierre Kraemer's avatar
Pierre Kraemer committed
643 644
		Algo::Topo::initAllOrbitsEmbedding<EDGE>(map);

Pierre Kraemer's avatar
Pierre Kraemer committed
645
	if (!map.template isOrbitEmbedded<FACE>())
Pierre Kraemer's avatar
Pierre Kraemer committed
646
		Algo::Topo::initAllOrbitsEmbedding<FACE>(map);
Sylvain Thery's avatar
Sylvain Thery committed
647

Sylvain Thery's avatar
Sylvain Thery committed
648
	CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
Sylvain Thery's avatar
Sylvain Thery committed
649
	{
650
		computeCurvatureVertex_NormalCycles<PFP>(map, v, radius, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
651
	}, FORCE_CELL_MARKING);
652 653
}

Sylvain Thery's avatar
Sylvain Thery committed
654
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
655
void computeCurvatureVertices_NormalCycles_Projected(
Sylvain Thery's avatar
Sylvain Thery committed
656
	typename PFP::MAP& map,
Sylvain Thery's avatar
Sylvain Thery committed
657
	typename PFP::REAL radius,
658 659 660
	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,
661
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
662 663 664 665 666
	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)
Sylvain Thery's avatar
Sylvain Thery committed
667
{
Pierre Kraemer's avatar
Pierre Kraemer committed
668
//	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
Sylvain Thery's avatar
Sylvain Thery committed
669
	if (!map.template isOrbitEmbedded<VERTEX>())
Pierre Kraemer's avatar
Pierre Kraemer committed
670 671
		Algo::Topo::initAllOrbitsEmbedding<VERTEX>(map);

Sylvain Thery's avatar
Sylvain Thery committed
672
	if (!map.template isOrbitEmbedded<EDGE>())
Pierre Kraemer's avatar
Pierre Kraemer committed
673 674
		Algo::Topo::initAllOrbitsEmbedding<EDGE>(map);

Sylvain Thery's avatar
Sylvain Thery committed
675
	if (!map.template isOrbitEmbedded<FACE>())
Pierre Kraemer's avatar
Pierre Kraemer committed
676
		Algo::Topo::initAllOrbitsEmbedding<FACE>(map);
Sylvain Thery's avatar
Sylvain Thery committed
677

Sylvain Thery's avatar
Sylvain Thery committed
678
	CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
Sylvain Thery's avatar
Sylvain Thery committed
679
	{
680
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, radius, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
681
	}, FORCE_CELL_MARKING);
Sylvain Thery's avatar
Sylvain Thery committed
682 683 684 685 686
}

template <typename PFP>
void computeCurvatureVertices_QuadraticFitting(
	typename PFP::MAP& map,
687 688 689 690 691 692
	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)
Sylvain Thery's avatar
Sylvain Thery committed
693
{
Sylvain Thery's avatar
Sylvain Thery committed
694
	CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
Sylvain Thery's avatar
Sylvain Thery committed
695
	{
Sylvain Thery's avatar
Sylvain Thery committed
696
		computeCurvatureVertex_QuadraticFitting<PFP>(map, v, position, normal, kmax, kmin, Kmax, Kmin) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
697
	}, FORCE_CELL_MARKING);
698 699
}

Sylvain Thery's avatar
Sylvain Thery committed
700
} // namespace Parallel
701 702


Pierre Kraemer's avatar
Pierre Kraemer committed
703 704
} // namespace Geometry

Pierre Kraemer's avatar
Pierre Kraemer committed
705
} // namespace Surface
706

Pierre Kraemer's avatar
Pierre Kraemer committed
707 708 709
} // namespace Algo

} // namespace CGoGN