qem.h 9.19 KB
Newer Older
Pierre Kraemer's avatar
Pierre Kraemer committed
1
/*******************************************************************************
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
 * CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
 * version 0.1                                                                  *
 * Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg           *
 *                                                                              *
 * 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.           *
 *                                                                              *
 * Web site: http://cgogn.unistra.fr/                                           *
 * Contact information: cgogn@unistra.fr                                        *
 *                                                                              *
 *******************************************************************************/
Pierre Kraemer's avatar
Pierre Kraemer committed
24

25
/*! \file qem.h
26 27 28
* Header file for Quadric Error Metric classes.
*/

Pierre Kraemer's avatar
Pierre Kraemer committed
29 30 31
#ifndef __QEM__
#define __QEM__

32
#include "Utils/os_spec.h" // allow compilation under windows
33
#include <cmath>
34

Pierre Kraemer's avatar
Pierre Kraemer committed
35 36
#include "Geometry/vector_gen.h"
#include "Geometry/matrix.h"
37
#include "Geometry/tensor.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
38 39
#include "Geometry/plane_3d.h"

40 41
// Eigen includes
#include <Eigen/Dense>
42

43
#define CONST_VAL -5212368.54127 // random value
44

45 46 47
/*! \namespace CGoGN
 * \brief namespace for all elements composing the CGoGN library
 */
48 49
namespace CGoGN {

50 51 52
/*! \namespace Utils
 * \brief namespace for tool classes used by CGoGN and its applications
 */
53
namespace Utils {
54

55 56 57 58 59
/*! \class Quadric
 *
 * \brief Quadric for computing the quadric error metric (QEM)
 * introduced by Garland and Heckbert in 1997.
 */
Pierre Kraemer's avatar
Pierre Kraemer committed
60 61 62 63
template <typename REAL>
class Quadric
{
public:
64
	/**
65 66
	 * \brief get CGoGN name of current type
	 *
67 68 69 70 71 72
	 * \return name of the CGoGN type
	 */
	static std::string CGoGNnameOfType()
	{
		return "Quadric" ;
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
73

74 75 76
	typedef Geom::Vector<3,REAL> VEC3 ;
	typedef Geom::Vector<4,REAL> VEC4 ;
	typedef Geom::Matrix<4,4,double> MATRIX44 ; // double is crucial here !
Pierre Kraemer's avatar
Pierre Kraemer committed
77

78 79 80 81 82
	/*!
	 * \brief Default constructor
	 *
	 * Initializes empty members
	 */
83
	Quadric() ;
Pierre Kraemer's avatar
Pierre Kraemer committed
84

85
	Quadric(int i) ;
86 87 88 89 90 91 92 93

	/*!
	 * \brief Constructor building a quadric given three points (defining a plane);
	 *
	 * \param p1 first point
	 * \param p2 second point
	 * \param p3 third point
	 */
94
	Quadric(VEC3& p1, VEC3& p2, VEC3& p3) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
95

96 97 98
	/*!
	 * \brief set members to zero
	 */
99
	void zero() ;
100

101 102 103 104 105
	/*!
	 * \brief affectation operator (by copy)
	 *
	 * \param q the Quadric to copy
	 */
106
	void operator= (const Quadric<REAL>& q) ;
107 108 109 110 111 112

	/*!
	 * \brief sum of Quadric operator
	 *
	 * \param q the Quadric to sum
	 */
113
	Quadric& operator+= (const Quadric<REAL>& q) ;
114 115 116 117 118 119

	/*!
	 * \brief substract of Quadric operator
	 *
	 * \param q the Quadric to substract
	 */
120
	Quadric& operator -= (const Quadric<REAL>& q) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
121

122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
	/*!
	 * \brief scalar product operator
	 *
	 * \param v the scalar to multiply the Quadric with
	 */
	Quadric& operator *= (const REAL& v) ;

	/*!
	 * \brief scalar division operator
	 *
	 * \param v the scalar to divide the Quadric with
	 */
	Quadric& operator /= (const REAL& v) ;

	/*!
	 * `brief error evaluation operator
	 *
	 * \param v a point expressed in homogeneous coordinates in space
	 *
	 * \return the error
	 */
143
	REAL operator() (const VEC4& v) const ;
144 145 146 147 148 149 150 151

	/*!
	 * `brief error evaluation operator
	 *
	 * \param v a point in space
	 *
	 * \param the error
	 */
152
	REAL operator() (const VEC3& v) const ;
Pierre Kraemer's avatar
Pierre Kraemer committed
153

154 155 156 157 158 159 160 161
	/*!
	 * \brief Write to stream operator
	 *
	 * \param out the stream to write to
	 * \param q the Quadric to write in the stream
	 *
	 * \return the stream reference
	 */
Pierre Kraemer's avatar
Pierre Kraemer committed
162 163 164 165
	friend std::ostream& operator<<(std::ostream& out, const Quadric<REAL>& q)
	{
		out << q.A ;
		return out ;
166
	} ;
Pierre Kraemer's avatar
Pierre Kraemer committed
167

168 169 170 171 172 173 174 175
	/*!
	 * \brief Read from stream operator
	 *
	 * \param int the stream to read from
	 * \param q the Quadric to write the data that has been read in
	 *
	 * \return the stream reference
	 */
Pierre Kraemer's avatar
Pierre Kraemer committed
176 177 178 179
	friend std::istream& operator>>(std::istream& in, Quadric<REAL>& q)
	{
		in >> q.A ;
		return in ;
180
	} ;
Pierre Kraemer's avatar
Pierre Kraemer committed
181

182 183 184 185 186 187 188
	/*!
	 * \brief Method to deduce a position in space that minimizes the error
	 *
	 * \param v the ideal position (if it can be computed)
	 *
	 * \return true if the ideal position has been computed correctly
	 */
189
	bool findOptimizedPos(VEC3& v) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
190 191

private:
192 193 194 195 196 197 198 199 200
	MATRIX44 A ; /*!< The Quadric matrix */

	/*!
	 * \brief method to evaluate the error at a given point in space (homogeneous coordinates)
	 *
	 * \param v the given point
	 *
	 * \return the error
	 */
201
	REAL evaluate(const VEC4& v) const ;
Pierre Kraemer's avatar
Pierre Kraemer committed
202

203 204 205 206 207 208 209 210
	/*!
	 * \brief method to deduce an optimal position in space (homogeneous coordinates)
	 * w.r.t. the current Quadric.
	 *
	 * \param v will contain the optimal position (if it can be computed)
	 *
	 * \return true if an optimal position was correctly computed
	 */
211
	bool optimize(VEC4& v) const ;
Pierre Kraemer's avatar
Pierre Kraemer committed
212 213
} ;

214 215 216 217 218 219
/*! \class QuadricNd
 * \brief extension of Quadric (which is 3D) to nD points.
 * This was published by Garland and Heckbert in 1998 and is meant to define a quadric
 * for a nD-space which contains geometry (3D) + other attributes like color (3D),
 * normals (3D) or texture coordinates (2D) for instance.
 */
220 221 222 223
template <typename REAL, unsigned int N>
class QuadricNd
{
public:
224 225 226 227 228
	/**
	 * \brief get CGoGN name of current type
	 *
	 * \return name of the CGoGN type
	 */
229 230 231 232
	static std::string CGoGNnameOfType()
	{
		return "QuadricNd" ;
	}
233 234 235 236

	typedef Geom::Vector<N,REAL> VECN ;
	typedef Geom::Vector<N+1,REAL> VECNp ;

237
	QuadricNd() ;
238

239
	QuadricNd(int i) ;
240

241
	QuadricNd(const VECN& p1_r, const VECN& p2_r, const VECN& p3_r) ;
242

243
	void zero() ;
244

245
	void operator= (const QuadricNd<REAL,N>& q) ;
246

247
	QuadricNd& operator+= (const QuadricNd<REAL,N>& q) ;
248

249
	QuadricNd& operator -= (const QuadricNd<REAL,N>& q) ;
250

251
	QuadricNd& operator *= (REAL v) ;
252

253
	QuadricNd& operator /= (REAL v) ;
254

255
	REAL operator() (const VECNp& v) const ;
256

257
	REAL operator() (const VECN& v) const ;
258 259 260

	friend std::ostream& operator<<(std::ostream& out, const QuadricNd<REAL,N>& q)
	{
261
		out << "(" << q.A << ", " << q.b << ", " << q.c << ")" ;
262
		return out ;
263
	} ;
264 265 266

	friend std::istream& operator>>(std::istream& in, QuadricNd<REAL,N>& q)
	{
267 268 269
		in >> q.A ;
		in >> q.b ;
		in >> q.c ;
270
		return in ;
271
	} ;
272

273
	bool findOptimizedVec(VECN& v) ;
274 275

private:
276 277 278 279
	// Double computation is crucial for stability
	Geom::Matrix<N,N,double> A ;
	Geom::Vector<N,double> b ;
	double c ;
280

281
	REAL evaluate(const VECN& v) const ;
282

283
	bool optimize(VECN& v) const ;
284 285
} ;

286 287 288 289
template <typename REAL>
class QuadricHF
{
public:
290 291 292 293 294
	/**
	 * \brief get CGoGN name of current type
	 *
	 * \return name of the CGoGN type
	 */
295 296
	static std::string CGoGNnameOfType() { return "QuadricHF" ; }

297
	typedef Geom::Vector<3,REAL> VEC3 ;
298

299
	QuadricHF() ;
300

301
	QuadricHF(int i) ;
302

303
	QuadricHF(const std::vector<VEC3>& v, const REAL& gamma, const REAL& alpha) ;
304

305
	QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REAL& alpha) ;
306

307
	~QuadricHF() ;
308

309 310 311 312 313 314 315 316 317 318 319 320 321 322
	void zero() ;

	QuadricHF& operator= (const QuadricHF<REAL>& q) ;

	QuadricHF& operator+= (const QuadricHF<REAL>& q) ;
	QuadricHF& operator -= (const QuadricHF<REAL>& q) ;

	QuadricHF& operator *= (const REAL& v) ;

	QuadricHF& operator /= (const REAL& v) ;

	//	REAL operator() (const VECNp& v) const ;

	REAL operator() (const std::vector<VEC3>& coefs) const ;
323 324 325

	friend std::ostream& operator<<(std::ostream& out, const QuadricHF<REAL>& q)
	{
326
		// TODO out << "(" << q.m_A << ", " << q.m_b << ", " << q.m_c << ")" ;
327
		return out ;
328
	} ;
329 330 331

	friend std::istream& operator>>(std::istream& in, QuadricHF<REAL>& q)
	{
332
		// TODO
333 334 335
		//		in >> q.A ;
		//		in >> q.b ;
		//		in >> q.c ;
336
		return in ;
337
	} ;
338

339
	bool findOptimizedCoefs(std::vector<VEC3>& coefs) ;
340 341 342

private:
	// Double computation is crucial for stability
343
	Eigen::MatrixXd m_A ;
344 345
	Eigen::VectorXd m_b[3] ;
	double m_c[3] ;
346

347
	REAL evaluate(const std::vector<VEC3>& coefs) const ;
348

349
	bool optimize(std::vector<VEC3>& coefs) const ;
350

351
	Geom::Matrix33d rotateMatrix(const REAL& gamma) ;
352

353
	Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ;
354

355 356
	Eigen::MatrixXd buildIntegralMatrix_A(const REAL& alpha, unsigned int size) ;
	Eigen::MatrixXd buildIntegralMatrix_B(const REAL& alpha, unsigned int size) ;
357

358
	Eigen::MatrixXd buildIntegralMatrix_C(const REAL& alpha, unsigned int size) ;
359

360
	void fillTensor(Geom::Tensor3d& T) ;
361

362 363 364
	Geom::Tensor3d* tensorsFromCoefs(const std::vector<VEC3>& coefs) ;

	std::vector<VEC3> coefsFromTensors(Geom::Tensor3d* A) ;
365 366
} ;

367
} // Utils
368 369 370

} // CGOGN

371 372
#include "qem.hpp"

Pierre Kraemer's avatar
Pierre Kraemer committed
373
#endif