curvature.hpp 26 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

Sylvain Thery's avatar
Sylvain Thery committed
30 31 32 33 34 35 36 37 38
extern "C"
{
#include "C_BLAS_LAPACK/INCLUDE/f2c.h"
#include "C_BLAS_LAPACK/INCLUDE/clapack.h"
}
#undef max
#undef min


Pierre Kraemer's avatar
Pierre Kraemer committed
39 40 41 42 43 44
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

84
	MATRIX33 localFrame = Algo::Geometry::vertexLocalFrame<PFP>(map, dart, 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, 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_(&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 96 97 98 99 100 101 102 103 104 105 106
	Eigen::Matrix<REAL,2,2> m;
	m << a, b, b, c;

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

	const Eigen::Matrix<REAL,2,1>& ev = solver.eigenvalues();
	kmax_v = ev[0];
	kmin_v = ev[1];

	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[dart] = -kmax_v ;
		kmin[dart] = -kmin_v ;
		Kmax[dart] = Kmax_v ;
		Kmin[dart] = Kmin_v ;
118 119 120
	}
	else
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
121 122 123 124
		kmax[dart] = -kmin_v ;
		kmin[dart] = -kmax_v ;
		Kmax[dart] = Kmin_v ;
		Kmin[dart] = 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 131 132
	typename PFP::MAP& map,
	Dart dart,
	typename PFP::MATRIX33& localFrame,
133 134
	const VertexAttribute<typename PFP::VEC3>& position,
	const VertexAttribute<typename PFP::VEC3>& normal,
Pierre Kraemer's avatar
Pierre Kraemer committed
135
	float& a, float& b, float& c, float& d, float& e)
Pierre Kraemer's avatar
Pierre Kraemer committed
136
{
Pierre Kraemer's avatar
Pierre Kraemer committed
137
	typename PFP::VEC3 p = position[dart] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
138

139
	LinearSolver<typename PFP::REAL> solver(5) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
140 141
	solver.set_least_squares(true) ;
	solver.begin_system() ;
142 143
	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
144
	{
145
		typename PFP::VEC3 v = position[it] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
146
		quadraticFittingAddVertexPos<PFP>(v, p, localFrame, solver) ;
147
		typename PFP::VEC3 n = normal[it] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
148
		quadraticFittingAddVertexNormal<PFP>(v, n, p, localFrame, solver) ;
149
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
150 151 152 153 154 155 156 157 158 159 160
	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>
161
void quadraticFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<typename PFP::REAL>& solver)
Pierre Kraemer's avatar
Pierre Kraemer committed
162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
{
	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>
178
void quadraticFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<typename PFP::REAL>& solver)
Pierre Kraemer's avatar
Pierre Kraemer committed
179 180 181 182 183 184
{
	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
185
	solver.add_coefficient(0, 2.0f * vec[0]) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
186 187 188 189
	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
190
	solver.set_right_hand_side(-1.0f * norm[0] / norm[2]) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
191 192 193 194 195
	solver.end_row() ;

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

template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
	typename PFP::REAL radius,
299 300 301 302 303 304 305 306
	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,
307 308
	const FunctorSelect& select,
	unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
309
{
Pierre Kraemer's avatar
Pierre Kraemer committed
310
	TraversorV<typename PFP::MAP> t(map, select) ;
311
	for(Dart d = t.begin(); d != t.end(); d = t.next())
312
		computeCurvatureVertex_NormalCycles<PFP>(map, d, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal,thread) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
313 314 315 316 317 318 319
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles(
	typename PFP::MAP& map,
	Dart dart,
	typename PFP::REAL radius,
320 321 322 323 324 325 326
	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,
327 328
	VertexAttribute<typename PFP::VEC3>& Knormal,
	unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
329
{
330 331
	typedef typename PFP::REAL REAL ;
	typedef typename PFP::VEC3 VEC3 ;
332 333
	typedef Geom::Matrix<3,3,REAL> MATRIX;
	typedef Eigen::Matrix<REAL,3,1> E_VEC3;
Sauvage's avatar
Sauvage committed
334
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;
335

336
	// collect the normal cycle tensor
337
	Selection::Collector_WithinSphere<PFP> neigh(map, position, radius, thread) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
338
	neigh.collectAll(dart) ;
339

340
	MATRIX tensor(0) ;
341
	neigh.computeNormalCyclesTensor(position, edgeangle,tensor);
342

343
	// solve eigen problem
344 345 346
	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());
347

348
	normalCycles_SortAndSetEigenComponents<PFP>(ev,evec,kmax[dart],kmin[dart],Kmax[dart],Kmin[dart],Knormal[dart],normal[dart],thread);
Sauvage's avatar
Sauvage committed
349 350 351 352 353 354 355 356 357 358

//	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;
//	}
359 360
}

Sauvage's avatar
Sauvage committed
361 362 363 364 365 366 367 368 369 370 371 372
template <typename PFP>
void computeCurvatureVertices_NormalCycles_Projected(
	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,
373 374
	const FunctorSelect& select,
	unsigned int thread)
Sauvage's avatar
Sauvage committed
375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
{
	TraversorV<typename PFP::MAP> t(map, select) ;
	for(Dart d = t.begin(); d != t.end(); d = t.next())
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, d, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal,thread) ;
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles_Projected(
	typename PFP::MAP& map,
	Dart dart,
	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,
393 394
	VertexAttribute<typename PFP::VEC3>& Knormal,
	unsigned int thread)
Sauvage's avatar
Sauvage committed
395
{
396
	typedef typename PFP::REAL REAL ;
Sauvage's avatar
Sauvage committed
397 398 399 400
	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;
401

402
	// collect the normal cycle tensor
403
	Selection::Collector_WithinSphere<PFP> neigh(map, position, radius, thread) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
404
	neigh.collectAll(dart) ;
405

Sauvage's avatar
Sauvage committed
406 407
	MATRIX tensor(0) ;
	neigh.computeNormalCyclesTensor(position, edgeangle,tensor);
408

Sauvage's avatar
Sauvage committed
409 410
	// project the tensor
	normalCycles_ProjectTensor<PFP>(tensor,normal[dart],thread);
411

Sauvage's avatar
Sauvage committed
412 413 414 415 416 417 418 419
	// 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[dart],kmin[dart],Kmax[dart],Kmin[dart],Knormal[dart],normal[dart],thread);
}

420 421 422
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
	typename PFP::MAP& map,
423
	Selection::Collector<PFP> & neigh,
424 425 426 427 428 429 430 431
	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,
432 433
	const FunctorSelect& select,
	unsigned int thread)
434 435 436 437 438 439 440 441 442 443
{
	TraversorV<typename PFP::MAP> t(map, select) ;
	for(Dart d = t.begin(); d != t.end(); d = t.next())
		computeCurvatureVertex_NormalCycles<PFP>(map, d, neigh, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal,thread) ;
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles(
	typename PFP::MAP& map,
	Dart dart,
444
	Selection::Collector<PFP> & neigh,
445 446 447 448 449 450 451
	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,
452 453
	VertexAttribute<typename PFP::VEC3>& Knormal,
	unsigned int thread)
454 455 456 457 458
{
	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
459
	typedef Eigen::Matrix<REAL,3,3,Eigen::RowMajor> E_MATRIX;
460 461 462

	// collect the normal cycle tensor
	neigh.collectAll(dart) ;
463

464 465
	MATRIX tensor(0) ;
	neigh.computeNormalCyclesTensor(position, edgeangle,tensor);
466

467
	// solve eigen problem
468 469 470 471 472 473 474
	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[dart],kmin[dart],Kmax[dart],Kmin[dart],Knormal[dart],normal[dart],thread);
}

Sauvage's avatar
Sauvage committed
475 476 477
template <typename PFP>
void computeCurvatureVertices_NormalCycles_Projected(
	typename PFP::MAP& map,
478
	Selection::Collector<PFP> & neigh,
Sauvage's avatar
Sauvage committed
479 480 481 482 483 484 485 486
	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,
487 488
	const FunctorSelect& select,
	unsigned int thread)
Sauvage's avatar
Sauvage committed
489 490 491 492 493 494 495 496 497 498
{
	TraversorV<typename PFP::MAP> t(map, select) ;
	for(Dart d = t.begin(); d != t.end(); d = t.next())
		computeCurvatureVertex_NormalCycles_Projected<PFP>(map, d, neigh, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal,thread) ;
}

template <typename PFP>
void computeCurvatureVertex_NormalCycles_Projected(
	typename PFP::MAP& map,
	Dart dart,
499
	Selection::Collector<PFP> & neigh,
Sauvage's avatar
Sauvage committed
500 501 502 503 504 505 506
	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,
507 508
	VertexAttribute<typename PFP::VEC3>& Knormal,
	unsigned int thread)
Sauvage's avatar
Sauvage committed
509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531
{
	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(dart) ;

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

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

	// 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[dart],kmin[dart],Kmax[dart],Kmin[dart],Knormal[dart],normal[dart],thread);
}
532

533 534 535 536 537 538 539 540 541 542
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,
543
	unsigned int thread)
544
{
Sauvage's avatar
Sauvage committed
545
	// sort eigen components : ev[inormal] has minimal absolute value ; kmin = ev[imin] <= ev[imax] = kmax
Sauvage's avatar
Sauvage committed
546
	int inormal=0, imin, imax ;
Sauvage's avatar
Sauvage committed
547 548 549 550
	if (abs(e_val[1]) < abs(e_val[inormal])) inormal = 1;
	if (abs(e_val[2]) < abs(e_val[inormal])) inormal = 2;
	imin = (inormal + 1) % 3;
	imax = (inormal + 2) % 3;
Sauvage's avatar
Sauvage committed
551
	if (e_val[imax] < e_val[imin]) { int tmp = imin ; imin = imax ; imax = tmp ; }
552

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

Sauvage's avatar
Sauvage committed
572
template <typename PFP>
573
void normalCycles_SortTensor(Geom::Matrix<3,3,typename PFP::REAL> & tensor, unsigned int thread)
Sauvage's avatar
Sauvage committed
574 575 576 577 578 579 580 581 582
{
	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
	Eigen::SelfAdjointEigenSolver<E_MATRIX> solver (Utils::convertRef<E_MATRIX>(tensor));
Sauvage's avatar
Sauvage committed
583 584
	const VEC3& e_val = Utils::convertRef<VEC3>(solver.eigenvalues());
	const MATRIX& e_vec = Utils::convertRef<MATRIX>(solver.eigenvectors());
Sauvage's avatar
Sauvage committed
585

Sauvage's avatar
Sauvage committed
586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602
	// switch kmin and kmax w.r.t. Kmin and Kmax
	int inormal=0, imin, imax ;
	if (abs(e_val[1]) < abs(e_val[inormal])) inormal = 1;
	if (abs(e_val[2]) < abs(e_val[inormal])) inormal = 2;
	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
603 604
}

Sauvage's avatar
Sauvage committed
605
template <typename PFP>
606
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
607 608 609 610 611
{
	Geom::Matrix<3,3,typename PFP::REAL> proj;
	proj.identity();
	proj -= Geom::transposed_vectors_mult(normal_vector,normal_vector);
	tensor = proj * tensor * proj;
Pierre Kraemer's avatar
Pierre Kraemer committed
612 613
}

614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649

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
650
	  m_Knormal(Knormal)
651 652
	 { }

653
	void run(Dart d, unsigned int threadID)
654
	{
Sylvain Thery's avatar
Sylvain Thery committed
655
		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) ;
656 657 658 659 660 661 662 663 664 665 666 667 668 669 670
	}
};

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,
671 672
	const FunctorSelect& select,
	unsigned int nbth)
673
{
Sylvain Thery's avatar
Sylvain Thery committed
674 675 676 677
	// WAHOO BIG PROBLEM WITH LAZZY EMBEDDING !!!
	if (!map. template isOrbitEmbedded<VERTEX>())
	{
		CellMarkerNoUnmark<VERTEX> cm(map);
678
		map. template initAllOrbitsEmbedding<VERTEX>();
Sylvain Thery's avatar
Sylvain Thery committed
679 680 681 682
	}
	if (!map. template isOrbitEmbedded<EDGE>())
	{
		CellMarkerNoUnmark<EDGE> cm(map);
683
		map. template initAllOrbitsEmbedding<EDGE>();
Sylvain Thery's avatar
Sylvain Thery committed
684 685 686 687
	}
	if (!map. template isOrbitEmbedded<FACE>())
	{
		CellMarkerNoUnmark<FACE> cm(map);
688
		map. template initAllOrbitsEmbedding<FACE>();
Sylvain Thery's avatar
Sylvain Thery committed
689 690
	}

691
	FunctorComputeCurvatureVertices_NormalCycles<PFP> funct(map, radius, position, normal, edgeangle, kmax, kmin, Kmax, Kmin, Knormal);
692
	Algo::Parallel::foreach_cell<typename PFP::MAP,VERTEX>(map, funct, nbth, true, select);
693 694
}

Sylvain Thery's avatar
Sylvain Thery committed
695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738


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,
739 740
	const FunctorSelect& select,
	unsigned int nbth)
Sylvain Thery's avatar
Sylvain Thery committed
741 742
{
	FunctorComputeCurvatureVertices_QuadraticFitting<PFP> funct(map, position, normal, kmax, kmin, Kmax, Kmin);
743
	Algo::Parallel::foreach_cell<typename PFP::MAP,VERTEX>(map, funct, nbth, true, select);
744 745
}

Sylvain Thery's avatar
Sylvain Thery committed
746
} // namespace Parallel
747 748


Pierre Kraemer's avatar
Pierre Kraemer committed
749 750
} // namespace Geometry

751 752
}

Pierre Kraemer's avatar
Pierre Kraemer committed
753 754 755
} // namespace Algo

} // namespace CGoGN