curvature.hpp 17.9 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/traversorCell.h"
#include "Topology/generic/traversor2.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
29
#include "Algo/Selection/collector.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

extern "C"
{
	#include "C_BLAS_LAPACK/INCLUDE/f2c.h"
	#include "C_BLAS_LAPACK/INCLUDE/clapack.h"
}
#undef max
#undef min

namespace CGoGN
{

namespace Algo
{

45 46 47
namespace Surface
{

Pierre Kraemer's avatar
Pierre Kraemer committed
48 49 50 51
namespace Geometry
{

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
52 53
void computeCurvatureVertices_QuadraticFitting(
	typename PFP::MAP& map,
54 55 56 57 58 59
	const VertexAttribute<typename PFP::VEC3>& position,
	const VertexAttribute<typename PFP::VEC3>& normal,
	VertexAttribute<typename PFP::REAL>& kmax,
	VertexAttribute<typename PFP::REAL>& kmin,
	VertexAttribute<typename PFP::VEC3>& Kmax,
	VertexAttribute<typename PFP::VEC3>& Kmin,
Pierre Kraemer's avatar
Pierre Kraemer committed
60
	const FunctorSelect& select)
Pierre Kraemer's avatar
Pierre Kraemer committed
61
{
Pierre Kraemer's avatar
Pierre Kraemer committed
62
	TraversorV<typename PFP::MAP> t(map, select) ;
63
	for(Dart d = t.begin(); d != t.end(); d = t.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
64
		computeCurvatureVertex_QuadraticFitting<PFP>(map, d, position, normal, kmax, kmin, Kmax, Kmin) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
65 66 67
}

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

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

Pierre Kraemer's avatar
Pierre Kraemer committed
84 85 86 87 88 89
	MATRIX33 localFrame = Algo::Geometry::vertexLocalFrame<PFP>(map, dart, position, n) ;
	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, dart, 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_(&e, &f, &g, &maxC, &minC, &dirX, &dirY) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
94
	/*int res = */slaev2_(&a, &b, &c, &kmax_v, &kmin_v, &Kmax_x, &Kmax_y) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
95

Pierre Kraemer's avatar
Pierre Kraemer committed
96 97 98
	VEC3 Kmax_v(Kmax_x, Kmax_y, 0.0f) ;
	Kmax_v = invLocalFrame * Kmax_v ;
	VEC3 Kmin_v = n ^ Kmax_v ;
Pierre Kraemer's avatar
Pierre Kraemer committed
99

Pierre Kraemer's avatar
Pierre Kraemer committed
100
	if (kmax_v < kmin_v)
101
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
102 103 104 105
		kmax[dart] = -kmax_v ;
		kmin[dart] = -kmin_v ;
		Kmax[dart] = Kmax_v ;
		Kmin[dart] = Kmin_v ;
106 107 108
	}
	else
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
109 110 111 112
		kmax[dart] = -kmin_v ;
		kmin[dart] = -kmax_v ;
		Kmax[dart] = Kmin_v ;
		Kmin[dart] = Kmax_v ;
113
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
114 115 116
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
117
void vertexQuadraticFitting(
Pierre Kraemer's avatar
Pierre Kraemer committed
118 119 120
	typename PFP::MAP& map,
	Dart dart,
	typename PFP::MATRIX33& localFrame,
121 122
	const VertexAttribute<typename PFP::VEC3>& position,
	const VertexAttribute<typename PFP::VEC3>& normal,
Pierre Kraemer's avatar
Pierre Kraemer committed
123
	float& a, float& b, float& c, float& d, float& e)
Pierre Kraemer's avatar
Pierre Kraemer committed
124
{
Pierre Kraemer's avatar
Pierre Kraemer committed
125
	typename PFP::VEC3 p = position[dart] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
126 127 128 129

	LinearSolver<CPUSolverTraits> solver(5) ;
	solver.set_least_squares(true) ;
	solver.begin_system() ;
130 131
	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
132
	{
133
		typename PFP::VEC3 v = position[it] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
134
		quadraticFittingAddVertexPos<PFP>(v, p, localFrame, solver) ;
135
		typename PFP::VEC3 n = normal[it] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
136
		quadraticFittingAddVertexNormal<PFP>(v, n, p, localFrame, solver) ;
137
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
138 139 140 141 142 143 144 145 146 147 148
	solver.end_system() ;
	solver.solve() ;

	a = solver.variable(0).value() ;
	b = solver.variable(1).value() ;
	c = solver.variable(2).value() ;
	d = solver.variable(3).value() ;
	e = solver.variable(4).value() ;
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
149
void quadraticFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
Pierre Kraemer's avatar
Pierre Kraemer committed
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
{
	typename PFP::VEC3 vec = v - p ;
	vec = localFrame * vec ;
	solver.begin_row() ;

	solver.add_coefficient(0, vec[0]*vec[0]) ;
	solver.add_coefficient(1, vec[0]*vec[1]) ;
	solver.add_coefficient(2, vec[1]*vec[1]) ;
	solver.add_coefficient(3, vec[0]) ;
	solver.add_coefficient(4, vec[1]) ;

	solver.set_right_hand_side(vec[2]) ;
	solver.end_row() ;
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
166
void quadraticFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
Pierre Kraemer's avatar
Pierre Kraemer committed
167 168 169 170 171 172
{
	typename PFP::VEC3 vec = v - p ;
	vec = localFrame * vec ;
	typename PFP::VEC3 norm = localFrame * n ;
	solver.begin_row() ;

Pierre Kraemer's avatar
Pierre Kraemer committed
173
	solver.add_coefficient(0, 2.0f * vec[0]) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
174 175 176 177
	solver.add_coefficient(1, vec[1]) ;
	solver.add_coefficient(2, 0) ;
	solver.add_coefficient(3, 1.0f) ;
	solver.add_coefficient(4, 0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
178
	solver.set_right_hand_side(-1.0f * norm[0] / norm[2]) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
179 180 181 182 183
	solver.end_row() ;

	solver.begin_row() ;
	solver.add_coefficient(0, 0) ;
	solver.add_coefficient(1, vec[0]) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
184
	solver.add_coefficient(2, 2.0f * vec[1]) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
185 186
	solver.add_coefficient(3, 0) ;
	solver.add_coefficient(4, 1.0f) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
187
	solver.set_right_hand_side(-1.0f * norm[1] / norm[2]) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
188 189 190 191
	solver.end_row() ;
}
/*
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
192
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
193 194 195 196 197 198 199 200
{
	gmtl::Matrix33f localFrame, invLocalFrame ;
	Algo::Geometry::vertexLocalFrame<PFP>(m_map,dart,normal,localFrame) ;
	gmtl::invertFull(invLocalFrame, localFrame) ;
	gmtl::Vec3f p = m_map.getVertexEmb(dart)->getPosition() ;
	solverC->reset(false) ;
	solverC->set_least_squares(true) ;
	solverC->begin_system() ;
201 202
	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
203 204
	{
		// 1-ring vertices
205
		gmtl::Vec3f v = m_map.getVertexEmb(it)->getPosition() ;
Pierre Kraemer's avatar
Pierre Kraemer committed
206
		cubicFittingAddVertexPos(v,p,localFrame) ;
207
		gmtl::Vec3f n = m_normalsV[m_map.getVertexEmb(it)->getLabel()] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
208
		cubicFittingAddVertexNormal(v,n,p,localFrame) ;
209
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227
	solverC->end_system() ;
	solverC->solve() ;
	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
228
void cubicFittingAddVertexPos(gmtl::Vec3f& v, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
Pierre Kraemer's avatar
Pierre Kraemer committed
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
{
	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
249
void cubicFittingAddVertexNormal(gmtl::Vec3f& v, gmtl::Vec3f& n, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
Pierre Kraemer's avatar
Pierre Kraemer committed
250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281
{
	gmtl::Vec3f vec = v - p ;
	vec = localFrame * vec ;
	gmtl::Vec3f norm = localFrame * n ;
	solverC->begin_row() ;

	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
282 283 284 285 286

template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
287 288 289 290 291 292 293 294
	const VertexAttribute<typename PFP::VEC3>& position,
	const VertexAttribute<typename PFP::VEC3>& normal,
	const EdgeAttribute<typename PFP::REAL>& edgeangle,
	VertexAttribute<typename PFP::REAL>& kmax,
	VertexAttribute<typename PFP::REAL>& kmin,
	VertexAttribute<typename PFP::VEC3>& Kmax,
	VertexAttribute<typename PFP::VEC3>& Kmin,
	VertexAttribute<typename PFP::VEC3>& Knormal,
295
	const FunctorSelect& select, unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
296
{
Pierre Kraemer's avatar
Pierre Kraemer committed
297
	TraversorV<typename PFP::MAP> t(map, select) ;
298
	for(Dart d = t.begin(); d != t.end(); d = t.next())
299
		computeCurvatureVertex_NormalCycles<PFP>(map, d, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal,thread) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
300 301 302 303 304 305 306
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles(
	typename PFP::MAP& map,
	Dart dart,
	typename PFP::REAL radius,
307 308 309 310 311 312 313
	const VertexAttribute<typename PFP::VEC3>& position,
	const VertexAttribute<typename PFP::VEC3>& normal,
	const EdgeAttribute<typename PFP::REAL>& edgeangle,
	VertexAttribute<typename PFP::REAL>& kmax,
	VertexAttribute<typename PFP::REAL>& kmin,
	VertexAttribute<typename PFP::VEC3>& Kmax,
	VertexAttribute<typename PFP::VEC3>& Kmin,
314
	VertexAttribute<typename PFP::VEC3>& Knormal, unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
315
{
316 317 318
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL ;

319
	// collect the normal cycle tensor
320
	Algo::Selection::Collector_WithinSphere<PFP> neigh(map, position, radius, thread) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
321
	neigh.collectAll(dart) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
322
	neigh.computeArea() ;
323

Pierre Kraemer's avatar
Pierre Kraemer committed
324
	VEC3 center = position[dart] ;
325 326 327

	typename PFP::MATRIX33 tensor(0) ;

328
	// collect edges inside the neighborhood
329 330 331
	const std::vector<Dart>& vd1 = neigh.getInsideEdges() ;
	for (std::vector<Dart>::const_iterator it = vd1.begin(); it != vd1.end(); ++it)
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
332
		const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(map, *it, position) ;
333
		tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) ;
334
	}
335
	// collect edges crossing the neighborhood's border
336 337 338
	const std::vector<Dart>& vd2 = neigh.getBorder() ;
	for (std::vector<Dart>::const_iterator it = vd2.begin(); it != vd2.end(); ++it)
	{
339 340 341
		const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(map, *it, position) ;
		REAL alpha ;
		Algo::Geometry::intersectionSphereEdge<PFP>(map, center, radius, *it, position, alpha) ;
342
		tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) * alpha ;
343 344 345 346
	}

	tensor /= neigh.getArea() ;

347 348
	// solve eigen problem
	Eigen::Matrix3f m3 ;
349
	m3 << tensor(0,0) , tensor(0,1) , tensor(0,2) , tensor(1,0) , tensor(1,1) , tensor(1,2) , tensor(2,0) , tensor(2,1) , tensor(2,2) ;
350 351 352 353 354 355 356 357 358 359 360
	Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> solver (m3);
	Eigen::Vector3f ev = solver.eigenvalues();
	Eigen::Matrix3f evec = solver.eigenvectors();

	// sort eigen components : ev[s[0]] has minimal absolute value ; kmin = ev[s[1]] <= ev[s[2]] = kmax
	int s[3] = {0, 1, 2} ;
	int tmp ;
	if (abs(ev[s[2]]) < abs(ev[s[1]])) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
	if (abs(ev[s[1]]) < abs(ev[s[0]])) { tmp = s[0] ; s[0] = s[1] ; s[1] = tmp ; }
	if (ev[s[2]] < ev[s[1]]) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }

361
	// set curvatures from sorted eigen components
362 363 364 365 366 367 368 369
	// warning : Kmin and Kmax are switched w.r.t. kmin and kmax
	// normal direction : minimal absolute eigen value
	VEC3& dirNormal = Knormal[dart] ;
	dirNormal[0] = evec(0,s[0]);
	dirNormal[1] = evec(1,s[0]);
	dirNormal[2] = evec(2,s[0]);
	if (dirNormal * normal[dart] < 0) dirNormal *= -1; // change orientation
	// min curvature
370 371 372 373
	kmin[dart] = ev[s[1]] ;
	VEC3& dirMin = Kmin[dart] ;
	dirMin[0] = evec(0,s[2]);
	dirMin[1] = evec(1,s[2]);
374 375 376
	dirMin[2] = evec(2,s[2]);
	// max curvature
	kmax[dart] = ev[s[2]] ;
377 378 379
	VEC3& dirMax = Kmax[dart] ;
	dirMax[0] = evec(0,s[1]);
	dirMax[1] = evec(1,s[1]);
380
	dirMax[2] = evec(2,s[1]);
Pierre Kraemer's avatar
Pierre Kraemer committed
381 382
}

383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418

namespace Parallel
{

template <typename PFP>
class FunctorComputeCurvatureVertices_NormalCycles: public FunctorMapThreaded<typename PFP::MAP >
{
	typename PFP::REAL m_radius;
	const VertexAttribute<typename PFP::VEC3>& m_position;
	const VertexAttribute<typename PFP::VEC3>& m_normal;
	const EdgeAttribute<typename PFP::REAL>& m_edgeangle;
	VertexAttribute<typename PFP::REAL>& m_kmax;
	VertexAttribute<typename PFP::REAL>& m_kmin;
	VertexAttribute<typename PFP::VEC3>& m_Kmax;
	VertexAttribute<typename PFP::VEC3>& m_Kmin;
	VertexAttribute<typename PFP::VEC3>& m_Knormal;
public:
	 FunctorComputeCurvatureVertices_NormalCycles( typename PFP::MAP& map,
		typename PFP::REAL radius,
		const VertexAttribute<typename PFP::VEC3>& position,
		const VertexAttribute<typename PFP::VEC3>& normal,
		const EdgeAttribute<typename PFP::REAL>& edgeangle,
		VertexAttribute<typename PFP::REAL>& kmax,
		VertexAttribute<typename PFP::REAL>& kmin,
		VertexAttribute<typename PFP::VEC3>& Kmax,
		VertexAttribute<typename PFP::VEC3>& Kmin,
		VertexAttribute<typename PFP::VEC3>& Knormal):
	  FunctorMapThreaded<typename PFP::MAP>(map),
	  m_radius(radius),
	  m_position(position),
	  m_normal(normal),
	  m_edgeangle(edgeangle),
	  m_kmax(kmax),
	  m_kmin(kmin),
	  m_Kmax(Kmax),
	  m_Kmin(Kmin),
Sylvain Thery's avatar
Sylvain Thery committed
419
	  m_Knormal(Knormal)
420 421
	 { }

422
	void run(Dart d, unsigned int threadID)
423
	{
Sylvain Thery's avatar
Sylvain Thery committed
424
		computeCurvatureVertex_NormalCycles<PFP>(this->m_map, d, m_radius, m_position, m_normal, m_edgeangle, m_kmax, m_kmin, m_Kmax, m_Kmin, m_Knormal, threadID) ;
425 426 427 428 429 430 431 432 433 434 435 436 437 438 439
	}
};

template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
	const VertexAttribute<typename PFP::VEC3>& position,
	const VertexAttribute<typename PFP::VEC3>& normal,
	const EdgeAttribute<typename PFP::REAL>& edgeangle,
	VertexAttribute<typename PFP::REAL>& kmax,
	VertexAttribute<typename PFP::REAL>& kmin,
	VertexAttribute<typename PFP::VEC3>& Kmax,
	VertexAttribute<typename PFP::VEC3>& Kmin,
	VertexAttribute<typename PFP::VEC3>& Knormal,
440
	const FunctorSelect& select, unsigned int nbth)
441
{
Sylvain Thery's avatar
Sylvain Thery committed
442 443 444 445
	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
	if (!map. template isOrbitEmbedded<VERTEX>())
	{
		CellMarkerNoUnmark<VERTEX> cm(map);
446
		map. template initAllOrbitsEmbedding<VERTEX>();
Sylvain Thery's avatar
Sylvain Thery committed
447 448 449 450
	}
	if (!map. template isOrbitEmbedded<EDGE>())
	{
		CellMarkerNoUnmark<EDGE> cm(map);
451
		map. template initAllOrbitsEmbedding<EDGE>();
Sylvain Thery's avatar
Sylvain Thery committed
452 453 454 455
	}
	if (!map. template isOrbitEmbedded<FACE>())
	{
		CellMarkerNoUnmark<FACE> cm(map);
456
		map. template initAllOrbitsEmbedding<FACE>();
Sylvain Thery's avatar
Sylvain Thery committed
457 458
	}

459
	FunctorComputeCurvatureVertices_NormalCycles<PFP> funct(map, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal);
460
	Algo::Parallel::foreach_cell<typename PFP::MAP,VERTEX>(map, funct, nbth, true, select);
461 462
}

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


template <typename PFP>
class FunctorComputeCurvatureVertices_QuadraticFitting: public FunctorMapThreaded<typename PFP::MAP >
{
	const VertexAttribute<typename PFP::VEC3>& m_position;
	const VertexAttribute<typename PFP::VEC3>& m_normal;
	VertexAttribute<typename PFP::REAL>& m_kmax;
	VertexAttribute<typename PFP::REAL>& m_kmin;
	VertexAttribute<typename PFP::VEC3>& m_Kmax;
	VertexAttribute<typename PFP::VEC3>& m_Kmin;
public:
	 FunctorComputeCurvatureVertices_QuadraticFitting( typename PFP::MAP& map,
		const VertexAttribute<typename PFP::VEC3>& position,
		const VertexAttribute<typename PFP::VEC3>& normal,
		VertexAttribute<typename PFP::REAL>& kmax,
		VertexAttribute<typename PFP::REAL>& kmin,
		VertexAttribute<typename PFP::VEC3>& Kmax,
		VertexAttribute<typename PFP::VEC3>& Kmin):
	  FunctorMapThreaded<typename PFP::MAP>(map),
	  m_position(position),
	  m_normal(normal),
	  m_kmax(kmax),
	  m_kmin(kmin),
	  m_Kmax(Kmax),
	  m_Kmin(Kmin)
	 { }

	void run(Dart d, unsigned int threadID)
	{
		computeCurvatureVertex_QuadraticFitting<PFP>(this->m_map, d, m_position, m_normal, m_kmax, m_kmin, m_Kmax, m_Kmin) ;
	}
};


template <typename PFP>
void computeCurvatureVertices_QuadraticFitting(
	typename PFP::MAP& map,
	const VertexAttribute<typename PFP::VEC3>& position,
	const VertexAttribute<typename PFP::VEC3>& normal,
	VertexAttribute<typename PFP::REAL>& kmax,
	VertexAttribute<typename PFP::REAL>& kmin,
	VertexAttribute<typename PFP::VEC3>& Kmax,
	VertexAttribute<typename PFP::VEC3>& Kmin,
507
	const FunctorSelect& select, unsigned int nbth)
Sylvain Thery's avatar
Sylvain Thery committed
508 509
{
	FunctorComputeCurvatureVertices_QuadraticFitting<PFP> funct(map, position, normal, kmax, kmin, Kmax, Kmin);
510
	Algo::Parallel::foreach_cell<typename PFP::MAP,VERTEX>(map, funct, nbth, true, select);
511 512
}

Sylvain Thery's avatar
Sylvain Thery committed
513
} // namespace Parallel
514 515


Pierre Kraemer's avatar
Pierre Kraemer committed
516 517
} // namespace Geometry

518 519
}

Pierre Kraemer's avatar
Pierre Kraemer committed
520 521 522
} // namespace Algo

} // namespace CGoGN