curvature.hpp 27.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,
Pierre Kraemer's avatar
Pierre Kraemer committed
45
46
47
48
49
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmax,
Sylvain Thery's avatar
Sylvain Thery committed
50
51
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmin,
	unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
52
{
Sylvain Thery's avatar
Sylvain Thery committed
53
54
55
56
57
58
59
60
61

	if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread==0))
	{
		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
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, thread);
Sylvain Thery's avatar
Sylvain Thery committed
65
66
67
68

//	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
69
70
71
}

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

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

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

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

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

	// 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
106
107
	REAL kmax_v = ev[0];
	REAL kmin_v = ev[1];
Pierre Kraemer's avatar
Pierre Kraemer committed
108
109
110

	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
111
	Kmax_v = invLocalFrame * Kmax_v ;
Pierre Kraemer's avatar
Pierre Kraemer committed
112
113
114
	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
115

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

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

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

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

174
175
176
177
178
179
180
181
	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
182
183
184
}

template <typename PFP>
185
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
186
187
188
189
190
{
	typename PFP::VEC3 vec = v - p ;
	vec = localFrame * vec ;
	typename PFP::VEC3 norm = localFrame * n ;

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

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

275
	solverC->begin_row() ;
Pierre Kraemer's avatar
Pierre Kraemer committed
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
301
	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
302
303
304
305
306

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

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

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

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

355
	// collect the normal cycle tensor
356
	Selection::Collector_WithinSphere<PFP> neigh(map, position, radius, thread) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
357
	neigh.collectAll(v) ;
358

359
	MATRIX tensor(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
360
	neigh.computeNormalCyclesTensor(position, edgeangle, tensor);
361

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

Pierre Kraemer's avatar
Pierre Kraemer committed
367
	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v],thread);
Sauvage's avatar
Sauvage committed
368
369
370
371
372
373
374
375
376
377

//	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;
//	}
378
379
}

Sauvage's avatar
Sauvage committed
380
381
382
383
template <typename PFP>
void computeCurvatureVertices_NormalCycles_Projected(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
Pierre Kraemer's avatar
Pierre Kraemer committed
384
385
386
387
388
389
390
391
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Knormal,
392
	unsigned int thread)
Sauvage's avatar
Sauvage committed
393
{
Sylvain Thery's avatar
Sylvain Thery committed
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
	if ((CGoGN::Parallel::NumberOfThreads > 1) && (thread==0))
	{
		Parallel::computeCurvatureVertices_NormalCycles_Projected<PFP>(map, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal);
		return;
	}

	foreach_cell<VERTEX>(map, [&] (Vertex v)
	{
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal,thread) ;
	}
	,FORCE_CELL_MARKING,thread);

//	TraversorV<typename PFP::MAP> t(map) ;
//	for(Vertex v = t.begin(); v.dart != t.end(); v = t.next())
//		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal, thread) ;
Sauvage's avatar
Sauvage committed
409
410
411
412
413
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles_Projected(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
414
	Vertex v,
Sauvage's avatar
Sauvage committed
415
	typename PFP::REAL radius,
Pierre Kraemer's avatar
Pierre Kraemer committed
416
417
418
419
420
421
422
423
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Knormal,
424
	unsigned int thread)
Sauvage's avatar
Sauvage committed
425
{
426
	typedef typename PFP::REAL REAL ;
Sauvage's avatar
Sauvage committed
427
428
429
430
	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;
431

432
	// collect the normal cycle tensor
433
	Selection::Collector_WithinSphere<PFP> neigh(map, position, radius, thread) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
434
	neigh.collectAll(v) ;
435

Sauvage's avatar
Sauvage committed
436
	MATRIX tensor(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
437
	neigh.computeNormalCyclesTensor(position, edgeangle, tensor);
438

Sauvage's avatar
Sauvage committed
439
	// project the tensor
Pierre Kraemer's avatar
Pierre Kraemer committed
440
	normalCycles_ProjectTensor<PFP>(tensor, normal[v], thread);
441

Sauvage's avatar
Sauvage committed
442
	// solve eigen problem
Pierre Kraemer's avatar
Pierre Kraemer committed
443
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
Sauvage's avatar
Sauvage committed
444
445
446
	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());

Pierre Kraemer's avatar
Pierre Kraemer committed
447
	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v],thread);
Sauvage's avatar
Sauvage committed
448
449
}

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

template <typename PFP>
void computeCurvatureVertex_NormalCycles(
Pierre Kraemer's avatar
Pierre Kraemer committed
471
472
	Vertex v,
	Algo::Surface::Selection::Collector<PFP>& neigh,
Pierre Kraemer's avatar
Pierre Kraemer committed
473
474
475
476
477
478
479
480
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Knormal,
481
	unsigned int thread)
482
483
484
485
486
{
	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
487
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;
488
489

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

492
	MATRIX tensor(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
493
	neigh.computeNormalCyclesTensor(position, edgeangle, tensor);
494

495
	// solve eigen problem
Pierre Kraemer's avatar
Pierre Kraemer committed
496
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
497
498
499
	const VEC3& ev = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& evec = Utils::convertRef<MATRIX>(solver.eigenvectors());

Pierre Kraemer's avatar
Pierre Kraemer committed
500
	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v],thread);
501
502
}

Sauvage's avatar
Sauvage committed
503
504
505
template <typename PFP>
void computeCurvatureVertices_NormalCycles_Projected(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
506
	Algo::Surface::Selection::Collector<PFP>& neigh,
Pierre Kraemer's avatar
Pierre Kraemer committed
507
508
509
510
511
512
513
514
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Knormal,
515
	unsigned int thread)
Sauvage's avatar
Sauvage committed
516
{
517
	TraversorV<typename PFP::MAP> t(map) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
518
519
	for(Vertex v = t.begin(); v != t.end(); v = t.next())
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, neigh, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal,thread) ;
Sauvage's avatar
Sauvage committed
520
521
522
523
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles_Projected(
Pierre Kraemer's avatar
Pierre Kraemer committed
524
525
	Vertex v,
	Algo::Surface::Selection::Collector<PFP>& neigh,
Pierre Kraemer's avatar
Pierre Kraemer committed
526
527
528
529
530
531
532
533
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Knormal,
534
	unsigned int thread)
Sauvage's avatar
Sauvage committed
535
536
537
538
539
540
541
542
{
	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
543
	neigh.collectAll(v) ;
Sauvage's avatar
Sauvage committed
544
545

	MATRIX tensor(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
546
	neigh.computeNormalCyclesTensor(position, edgeangle, tensor);
Sauvage's avatar
Sauvage committed
547
548

	// project the tensor
Pierre Kraemer's avatar
Pierre Kraemer committed
549
	normalCycles_ProjectTensor<PFP>(tensor, normal[v], thread);
Sauvage's avatar
Sauvage committed
550
551

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

Pierre Kraemer's avatar
Pierre Kraemer committed
556
	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[v],kmin[v],Kmax[v],Kmin[v],Knormal[v],normal[v],thread);
Sauvage's avatar
Sauvage committed
557
}
558

559
560
561
562
563
564
565
566
567
568
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,
	const typename PFP::VEC3& normal,
Sylvain Thery's avatar
Sylvain Thery committed
569
	unsigned int /*thread*/)
570
{
Sauvage's avatar
Sauvage committed
571
	// sort eigen components : ev[inormal] has minimal absolute value ; kmin = ev[imin] <= ev[imax] = kmax
Pierre Kraemer's avatar
Pierre Kraemer committed
572
	int inormal = 0, imin, imax ;
Pierre Kraemer's avatar
Pierre Kraemer committed
573
574
	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
575
576
	imin = (inormal + 1) % 3;
	imax = (inormal + 2) % 3;
Sauvage's avatar
Sauvage committed
577
	if (e_val[imax] < e_val[imin]) { int tmp = imin ; imin = imax ; imax = tmp ; }
578

579
	// set curvatures from sorted eigen components
580
581
	// warning : Kmin and Kmax are switched w.r.t. kmin and kmax
	// normal direction : minimal absolute eigen value
Sauvage's avatar
Sauvage committed
582
583
584
	Knormal[0] = e_vec(0,inormal);
	Knormal[1] = e_vec(1,inormal);
	Knormal[2] = e_vec(2,inormal);
585
	if (Knormal * normal < 0) Knormal *= -1; // change orientation
586
	// min curvature
Sauvage's avatar
Sauvage committed
587
588
589
590
	kmin = e_val[imin] ;
	Kmin[0] = e_vec(0,imax);
	Kmin[1] = e_vec(1,imax);
	Kmin[2] = e_vec(2,imax);
591
	// max curvature
Sauvage's avatar
Sauvage committed
592
593
594
595
	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
596
597
}

Sauvage's avatar
Sauvage committed
598
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
599
void normalCycles_SortTensor(Geom::Matrix<3,3,typename PFP::REAL>& tensor, unsigned int /*thread*/)
Sauvage's avatar
Sauvage committed
600
601
602
603
604
605
606
607
{
	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
608
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver(Utils::convertRef<E_MATRIX>(tensor));
Sauvage's avatar
Sauvage committed
609
610
	const VEC3& e_val = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& e_vec = Utils::convertRef<MATRIX>(solver.eigenvectors());
Sauvage's avatar
Sauvage committed
611

Sauvage's avatar
Sauvage committed
612
	// switch kmin and kmax w.r.t. Kmin and Kmax
Pierre Kraemer's avatar
Pierre Kraemer committed
613
	int inormal = 0, imin, imax ;
Pierre Kraemer's avatar
Pierre Kraemer committed
614
615
	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
616
617
618
619
620
621
622
623
624
625
626
627
628
	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
629
630
}

Sauvage's avatar
Sauvage committed
631
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
632
void normalCycles_ProjectTensor(Geom::Matrix<3,3,typename PFP::REAL>& tensor, const typename PFP::VEC3& normal_vector, unsigned int thread)
Sauvage's avatar
Sauvage committed
633
634
635
{
	Geom::Matrix<3,3,typename PFP::REAL> proj;
	proj.identity();
Pierre Kraemer's avatar
Pierre Kraemer committed
636
	proj -= Geom::transposed_vectors_mult(normal_vector, normal_vector);
Sauvage's avatar
Sauvage committed
637
	tensor = proj * tensor * proj;
Pierre Kraemer's avatar
Pierre Kraemer committed
638
639
}

640
641
642

namespace Parallel
{
Sylvain Thery's avatar
Sylvain Thery committed
643

644
645
646
647
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
Pierre Kraemer's avatar
Pierre Kraemer committed
648
649
650
651
652
653
654
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal,
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeangle,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
655
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Knormal)
656
{
Sylvain Thery's avatar
Sylvain Thery committed
657
	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
Pierre Kraemer's avatar
Pierre Kraemer committed
658
	if (!map.template isOrbitEmbedded<VERTEX>())
Sylvain Thery's avatar
Sylvain Thery committed
659
660
		Algo::Topo::initAllOrbitsEmbedding<VERTEX>(map);

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

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

Sylvain Thery's avatar
Sylvain Thery committed
667
	CGoGN::Parallel::foreach_cell<VERTEX>(map,[&](Vertex v, unsigned int threadID)
Sylvain Thery's avatar
Sylvain Thery committed
668
669
	{
		computeCurvatureVertex_NormalCycles<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal, threadID) ;
Sylvain Thery's avatar
Sylvain Thery committed
670
	},true,FORCE_CELL_MARKING);
671
672
}

Sylvain Thery's avatar
Sylvain Thery committed
673
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
674
void computeCurvatureVertices_NormalCycles_Projected(
Sylvain Thery's avatar
Sylvain Thery committed
675
	typename PFP::MAP& map,
Sylvain Thery's avatar
Sylvain Thery committed
676
	typename PFP::REAL radius,
Pierre Kraemer's avatar
Pierre Kraemer committed
677
678
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal,
Sylvain Thery's avatar
Sylvain Thery committed
679
	const EdgeAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& edgeangle,
Pierre Kraemer's avatar
Pierre Kraemer committed
680
681
682
683
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmin,
Sylvain Thery's avatar
Sylvain Thery committed
684
685
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Knormal)
{
Sylvain Thery's avatar
Sylvain Thery committed
686
//	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
Sylvain Thery's avatar
Sylvain Thery committed
687
	if (!map.template isOrbitEmbedded<VERTEX>())
Sylvain Thery's avatar
Sylvain Thery committed
688
689
		Algo::Topo::initAllOrbitsEmbedding<VERTEX>(map);

Sylvain Thery's avatar
Sylvain Thery committed
690
	if (!map.template isOrbitEmbedded<EDGE>())
Sylvain Thery's avatar
Sylvain Thery committed
691
692
		Algo::Topo::initAllOrbitsEmbedding<EDGE>(map);

Sylvain Thery's avatar
Sylvain Thery committed
693
	if (!map.template isOrbitEmbedded<FACE>())
Sylvain Thery's avatar
Sylvain Thery committed
694
		Algo::Topo::initAllOrbitsEmbedding<FACE>(map);
Sylvain Thery's avatar
Sylvain Thery committed
695

Sylvain Thery's avatar
Sylvain Thery committed
696
	CGoGN::Parallel::foreach_cell<VERTEX>(map,[&](Vertex v, unsigned int threadID)
Sylvain Thery's avatar
Sylvain Thery committed
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
	{
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, v, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal, threadID) ;
	},true,FORCE_CELL_MARKING);
}



template <typename PFP>
void computeCurvatureVertices_QuadraticFitting(
	typename PFP::MAP& map,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmax,
	VertexAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& kmin,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmax,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& Kmin)
Sylvain Thery's avatar
Sylvain Thery committed
713
{
Sylvain Thery's avatar
Sylvain Thery committed
714
	CGoGN::Parallel::foreach_cell<VERTEX>(map,[&](Vertex v, unsigned int threadID)
Sylvain Thery's avatar
Sylvain Thery committed
715
716
	{
		computeCurvatureVertex_QuadraticFitting<PFP>(map, v, position, normal, kmax, kmin, Kmax, Kmin, threadID) ;
Sylvain Thery's avatar
Sylvain Thery committed
717
	},true,FORCE_CELL_MARKING);
Sylvain Thery's avatar
Sylvain Thery committed
718

719
720
}

Sylvain Thery's avatar
Sylvain Thery committed
721
} // namespace Parallel
722
723


Pierre Kraemer's avatar
Pierre Kraemer committed
724
725
} // namespace Geometry

Pierre Kraemer's avatar
Pierre Kraemer committed
726
} // namespace Surface
727

Pierre Kraemer's avatar
Pierre Kraemer committed
728
729
730
} // namespace Algo

} // namespace CGoGN