curvature.hpp 29.7 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
	// TODO: nl not thread safe
Sylvain Thery's avatar
Sylvain Thery committed
53

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

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

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
68 69
void computeCurvatureVertex_QuadraticFitting(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
70
	Vertex v,
71 72 73 74 75 76
	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
77 78 79 80 81
{
	typedef typename PFP::REAL REAL ;
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::MATRIX33 MATRIX33 ;

Pierre Kraemer's avatar
Pierre Kraemer committed
82
	VEC3 n = normal[v] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
83

Pierre Kraemer's avatar
Pierre Kraemer committed
84
	MATRIX33 localFrame = Algo::Geometry::vertexLocalFrame<PFP>(map, v, position, n) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
85 86 87 88 89
	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
90
	vertexQuadraticFitting<PFP>(map, v, localFrame, position, normal, a, b, c, d, e) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
91

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

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

	// 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
102 103
	REAL kmax_v = ev[0];
	REAL kmin_v = ev[1];
Pierre Kraemer's avatar
Pierre Kraemer committed
104 105 106

	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
107
	Kmax_v = invLocalFrame * Kmax_v ;
Pierre Kraemer's avatar
Pierre Kraemer committed
108 109 110
	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
111

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

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
129
void vertexQuadraticFitting(
Pierre Kraemer's avatar
Pierre Kraemer committed
130
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
131
	Vertex v,
Pierre Kraemer's avatar
Pierre Kraemer committed
132
	typename PFP::MATRIX33& localFrame,
133 134
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
135
    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
136
{
Pierre Kraemer's avatar
Pierre Kraemer committed
137
	typename PFP::VEC3 p = position[v] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
138

139 140 141 142 143 144
	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
145 146
	foreach_adjacent2<EDGE>(map, v, [&] (Vertex it) {
		typename PFP::VEC3 itp = position[it] ;
147
        quadraticFittingAddVertexPos<PFP>(p, itp, localFrame) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
148 149 150
		typename PFP::VEC3 itn = normal[it] ;
		quadraticFittingAddVertexNormal<PFP>(itp, itn, p, localFrame) ;
	});
151 152 153 154 155 156 157 158 159 160 161
	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
162 163 164
}

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

170 171 172 173 174 175 176 177
	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
178 179 180
}

template <typename PFP>
181
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
182 183 184 185 186
{
	typename PFP::VEC3 vec = v - p ;
	vec = localFrame * vec ;
	typename PFP::VEC3 norm = localFrame * n ;

187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203
	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
204 205 206
}
/*
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
207
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
208 209
{
	gmtl::Matrix33f localFrame, invLocalFrame ;
210
	Geometry::vertexLocalFrame<PFP>(m_map,dart,normal,localFrame) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
211 212 213 214 215
	gmtl::invertFull(invLocalFrame, localFrame) ;
	gmtl::Vec3f p = m_map.getVertexEmb(dart)->getPosition() ;
	solverC->reset(false) ;
	solverC->set_least_squares(true) ;
	solverC->begin_system() ;
216 217
	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
218 219
	{
		// 1-ring vertices
220
		gmtl::Vec3f v = m_map.getVertexEmb(it)->getPosition() ;
Pierre Kraemer's avatar
Pierre Kraemer committed
221
		cubicFittingAddVertexPos(v,p,localFrame) ;
222
		gmtl::Vec3f n = m_normalsV[m_map.getVertexEmb(it)->getLabel()] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
223
		cubicFittingAddVertexNormal(v,n,p,localFrame) ;
224
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
225 226
	solverC->end_system() ;
	solverC->solve() ;
227

Pierre Kraemer's avatar
Pierre Kraemer committed
228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
	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
244
void cubicFittingAddVertexPos(gmtl::Vec3f& v, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
Pierre Kraemer's avatar
Pierre Kraemer committed
245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264
{
	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
265
void cubicFittingAddVertexNormal(gmtl::Vec3f& v, gmtl::Vec3f& n, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
Pierre Kraemer's avatar
Pierre Kraemer committed
266 267 268 269 270
{
	gmtl::Vec3f vec = v - p ;
	vec = localFrame * vec ;
	gmtl::Vec3f norm = localFrame * n ;

271
	solverC->begin_row() ;
Pierre Kraemer's avatar
Pierre Kraemer committed
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297
	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
298 299 300 301 302

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

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

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

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

351
	MATRIX tensor(0) ;
352
	normalCycles_computeTensor(neigh, position, edgeangle, edgearea, tensor);
353

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

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

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

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

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

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

Sauvage's avatar
Sauvage committed
414
	MATRIX tensor(0) ;
415
	normalCycles_computeTensor(neigh, position, edgeangle, edgearea, tensor);
416

Sauvage's avatar
Sauvage committed
417
	// project the tensor
Sylvain Thery's avatar
Sylvain Thery committed
418
	normalCycles_ProjectTensor<PFP>(tensor, normal[v]);
419

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

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

428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535
//template <typename PFP>
//void computeCurvatureVertices_NormalCycles(
//	typename PFP::MAP& map,
//	Algo::Surface::Selection::Collector<PFP>& neigh,
//	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,
//	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
//	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)
//{
//	TraversorV<typename PFP::MAP> t(map) ;
//	for(Vertex v = t.begin(); v != t.end(); v = t.next())
//		computeCurvatureVertex_NormalCycles<PFP>(map, v, neigh, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal) ;
//}

//template <typename PFP>
//void computeCurvatureVertex_NormalCycles(
//	Vertex v,
//	Algo::Surface::Selection::Collector<PFP>& neigh,
//	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,
//	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
//	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)
//{
//	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
//	neigh.collectAll(v) ;

//	MATRIX tensor(0) ;
//	neigh.computeNormalCyclesTensor(position, edgeangle, edgearea, tensor);

//	// solve eigen problem
//	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
//	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
//	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());

//	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v]);
//}

//template <typename PFP>
//void computeCurvatureVertices_NormalCycles_Projected(
//	typename PFP::MAP& map,
//	Algo::Surface::Selection::Collector<PFP>& neigh,
//	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,
//	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
//	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)
//{
//	TraversorV<typename PFP::MAP> t(map) ;
//	for(Vertex v = t.begin(); v != t.end(); v = t.next())
//		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, neigh, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal) ;
//}

//template <typename PFP>
//void computeCurvatureVertex_NormalCycles_Projected(
//	Vertex v,
//	Algo::Surface::Selection::Collector<PFP>& neigh,
//	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,
//	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
//	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)
//{
//	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
//	neigh.collectAll(v) ;

//	MATRIX tensor(0) ;
//	neigh.computeNormalCyclesTensor(position, edgeangle, edgearea, tensor);

//	// project the tensor
//	normalCycles_ProjectTensor<PFP>(tensor, normal[v]);

//	// solve eigen problem
//	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
//	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
//	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());

//	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v]);
//}
536 537

template <typename PFP>
538 539
void normalCycles_computeTensor(
	Algo::Surface::Selection::Collector<PFP>& col,
540
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
541
	Geom::Matrix<3,3,typename PFP::REAL>& tensor)
542
{
543
	tensor.zero();
544

545 546 547 548 549 550 551
	// collect edges inside the neighborhood
	for (Edge e : col.getInsideEdges())
	{
		typename PFP::REAL edgeangle = Algo::Surface::Geometry::computeAngleBetweenNormalsOnEdge<PFP>(col.getMap(), e, position);
		typename PFP::VEC3 ev = Algo::Surface::Geometry::vectorOutOfDart<PFP>(col.getMap(), e, position);
		tensor += Geom::transposed_vectors_mult(ev,ev) * edgeangle * (1 / ev.norm());
	}
552

553 554 555 556 557 558 559 560
	// collect edges on the border
	for (Dart d : col.getBorder())
	{
		typename PFP::REAL edgeangle = Algo::Surface::Geometry::computeAngleBetweenNormalsOnEdge<PFP>(col.getMap(), d, position);
		typename PFP::REAL alpha = col.borderEdgeRatio(d, position);
		typename PFP::VEC3 ev = Algo::Surface::Geometry::vectorOutOfDart<PFP>(col.getMap(), d, position);
		tensor += Geom::transposed_vectors_mult(ev,ev) * edgeangle * (1 / ev.norm()) * alpha;
	}
561

562
	tensor /= col.computeArea(position);
563 564
}

Sauvage's avatar
Sauvage committed
565
template <typename PFP>
566 567
void normalCycles_computeTensor(
	Algo::Surface::Selection::Collector<PFP>& col,
568 569
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
570
	Geom::Matrix<3,3,typename PFP::REAL>& tensor)
Sauvage's avatar
Sauvage committed
571
{
572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589
	tensor.zero();

	// collect edges inside the neighborhood
	for (Edge e : col.getInsideEdges())
	{
		typename PFP::VEC3 ev = Algo::Surface::Geometry::vectorOutOfDart<PFP>(col.getMap(), e, position);
		tensor += Geom::transposed_vectors_mult(ev,ev) * edgeangle[e] * (1 / ev.norm());
	}

	// collect edges on the border
	for (Dart d : col.getBorder())
	{
		typename PFP::REAL alpha = col.borderEdgeRatio(d, position);
		typename PFP::VEC3 ev = Algo::Surface::Geometry::vectorOutOfDart<PFP>(col.getMap(), d, position);
		tensor += Geom::transposed_vectors_mult(ev,ev) * edgeangle[d] * (1 / ev.norm()) * alpha;
	}

	tensor /= col.computeArea(position);
Sauvage's avatar
Sauvage committed
590 591 592
}

template <typename PFP>
593 594
void normalCycles_computeTensor(
	Algo::Surface::Selection::Collector<PFP>& col,
595 596
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
597
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
598
	Geom::Matrix<3,3,typename PFP::REAL>& tensor)
Sauvage's avatar
Sauvage committed
599
{
600
	tensor.zero();
Sauvage's avatar
Sauvage committed
601

602 603 604 605 606 607
	// collect edges inside the neighborhood
	for (Edge e : col.getInsideEdges())
	{
		typename PFP::VEC3 ev = Algo::Surface::Geometry::vectorOutOfDart<PFP>(col.getMap(), e, position);
		tensor += Geom::transposed_vectors_mult(ev,ev) * edgeangle[e] * (1 / ev.norm());
	}
Sauvage's avatar
Sauvage committed
608

609 610 611 612 613 614 615
	// collect edges on the border
	for (Dart d : col.getBorder())
	{
		typename PFP::REAL alpha = col.borderEdgeRatio(d, position);
		typename PFP::VEC3 ev = Algo::Surface::Geometry::vectorOutOfDart<PFP>(col.getMap(), d, position);
		tensor += Geom::transposed_vectors_mult(ev,ev) * edgeangle[d] * (1 / ev.norm()) * alpha;
	}
Sauvage's avatar
Sauvage committed
616

617
	tensor /= col.computeArea(position, edgearea);
Sauvage's avatar
Sauvage committed
618
}
619

620 621 622 623 624 625 626 627 628
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
629
	const typename PFP::VEC3& normal)
630
{
Sauvage's avatar
Sauvage committed
631
	// sort eigen components : ev[inormal] has minimal absolute value ; kmin = ev[imin] <= ev[imax] = kmax
632
	int inormal = 0, imin, imax;
Pierre Kraemer's avatar
Pierre Kraemer committed
633 634
	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
635 636
	imin = (inormal + 1) % 3;
	imax = (inormal + 2) % 3;
637
	if (e_val[imax] < e_val[imin]) { int tmp = imin; imin = imax; imax = tmp; }
638

639
	// set curvatures from sorted eigen components
640 641
	// warning : Kmin and Kmax are switched w.r.t. kmin and kmax
	// normal direction : minimal absolute eigen value
Sauvage's avatar
Sauvage committed
642 643 644
	Knormal[0] = e_vec(0,inormal);
	Knormal[1] = e_vec(1,inormal);
	Knormal[2] = e_vec(2,inormal);
645
	if (Knormal * normal < 0) Knormal *= -1; // change orientation
646
	// min curvature
Sauvage's avatar
Sauvage committed
647 648 649 650
	kmin = e_val[imin] ;
	Kmin[0] = e_vec(0,imax);
	Kmin[1] = e_vec(1,imax);
	Kmin[2] = e_vec(2,imax);
651
	// max curvature
Sauvage's avatar
Sauvage committed
652 653 654 655
	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
656 657
}

Sauvage's avatar
Sauvage committed
658
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
659
void normalCycles_SortTensor(Geom::Matrix<3,3,typename PFP::REAL>& tensor)
Sauvage's avatar
Sauvage committed
660 661 662 663 664 665 666 667
{
	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
668
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
Sauvage's avatar
Sauvage committed
669 670
	const VEC3& e_val = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& e_vec = Utils::convertRef<MATRIX>(solver.eigenvectors());
Sauvage's avatar
Sauvage committed
671

Sauvage's avatar
Sauvage committed
672
	// switch kmin and kmax w.r.t. Kmin and Kmax
Pierre Kraemer's avatar
Pierre Kraemer committed
673
	int inormal = 0, imin, imax ;
Pierre Kraemer's avatar
Pierre Kraemer committed
674 675
	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
676 677 678 679 680 681 682 683 684 685 686 687 688
	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
689 690
}

Sauvage's avatar
Sauvage committed
691
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
692
void normalCycles_ProjectTensor(Geom::Matrix<3,3,typename PFP::REAL>& tensor, const typename PFP::VEC3& normal_vector)
Sauvage's avatar
Sauvage committed
693 694 695
{
	Geom::Matrix<3,3,typename PFP::REAL> proj;
	proj.identity();
Pierre Kraemer's avatar
Pierre Kraemer committed
696
	proj -= Geom::transposed_vectors_mult(normal_vector, normal_vector);
Sauvage's avatar
Sauvage committed
697
	tensor = proj * tensor * proj;
Pierre Kraemer's avatar
Pierre Kraemer committed
698 699
}

700 701 702

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

704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719
//template <typename PFP>
//void computeCurvatureVertices_QuadraticFitting(
//	typename PFP::MAP& map,
//	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)
//{
//	CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
//	{
//		computeCurvatureVertex_QuadraticFitting<PFP>(map, v, position, normal, kmax, kmin, Kmax, Kmin) ;
//	}, FORCE_CELL_MARKING);
//}

720 721 722 723
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
724 725 726
	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,
727
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
728 729 730 731 732
	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)
733
{
Sylvain Thery's avatar
Sylvain Thery committed
734
	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
Pierre Kraemer's avatar
Pierre Kraemer committed
735
	if (!map.template isOrbitEmbedded<VERTEX>())
Pierre Kraemer's avatar
Pierre Kraemer committed
736 737
		Algo::Topo::initAllOrbitsEmbedding<VERTEX>(map);

Pierre Kraemer's avatar
Pierre Kraemer committed
738
	if (!map.template isOrbitEmbedded<EDGE>())
Pierre Kraemer's avatar
Pierre Kraemer committed
739 740
		Algo::Topo::initAllOrbitsEmbedding<EDGE>(map);

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

Sylvain Thery's avatar
Sylvain Thery committed
744
	CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
Sylvain Thery's avatar
Sylvain Thery committed
745
	{
746
		computeCurvatureVertex_NormalCycles<PFP>(map, v, radius, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
747
	}, FORCE_CELL_MARKING);
748 749
}

Sylvain Thery's avatar
Sylvain Thery committed
750
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
751
void computeCurvatureVertices_NormalCycles_Projected(
Sylvain Thery's avatar
Sylvain Thery committed
752
	typename PFP::MAP& map,
Sylvain Thery's avatar
Sylvain Thery committed
753
	typename PFP::REAL radius,
754 755 756
	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,
757
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
758 759 760 761 762
	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
763
{
764
	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
Sylvain Thery's avatar
Sylvain Thery committed
765
	if (!map.template isOrbitEmbedded<VERTEX>())
Pierre Kraemer's avatar
Pierre Kraemer committed
766 767
		Algo::Topo::initAllOrbitsEmbedding<VERTEX>(map);

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

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

Sylvain Thery's avatar
Sylvain Thery committed
774
	CGoGN::Parallel::foreach_cell<VERTEX>(map, [&] (Vertex v, unsigned int /*thr*/)
Sylvain Thery's avatar
Sylvain Thery committed
775
	{
776
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, radius, position, normal, edgeangle, edgearea, kmax, kmin, Kmax, Kmin, Knormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
777
	}, FORCE_CELL_MARKING);
Sylvain Thery's avatar
Sylvain Thery committed
778 779
}

Sylvain Thery's avatar
Sylvain Thery committed
780
} // namespace Parallel
781 782


Pierre Kraemer's avatar
Pierre Kraemer committed
783 784
} // namespace Geometry

Pierre Kraemer's avatar
Pierre Kraemer committed
785
} // namespace Surface
786

Pierre Kraemer's avatar
Pierre Kraemer committed
787 788 789
} // namespace Algo

} // namespace CGoGN