qem.h 6.29 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
26
27
28
29
30
31
/**
* @file qem.h
* Header file for Quadric Error Metric classes.
*/



Pierre Kraemer's avatar
Pierre Kraemer committed
32
33
34
#ifndef __QEM__
#define __QEM__

35
#include "Utils/os_spec.h" // allow compilation under windows
36
#include <cmath>
37

Pierre Kraemer's avatar
Pierre Kraemer committed
38
39
#include "Geometry/vector_gen.h"
#include "Geometry/matrix.h"
40
#include "Geometry/tensor.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
41
42
#include "Geometry/plane_3d.h"

43
44
// Eigen includes
#include <Eigen/Dense>
45

46
#define CONST_VAL -5212368.54127 // random value
47

48
49
namespace CGoGN {

50
namespace Utils {
51

Pierre Kraemer's avatar
Pierre Kraemer committed
52
53
54
55
56
57
template <typename REAL>
class Quadric
{
public:
	static std::string CGoGNnameOfType() { return "Quadric" ; }

58
59
60
	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
61

62
63
	Quadric() ;
	Quadric(int i) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
64

65
	Quadric(VEC3& p1, VEC3& p2, VEC3& p3) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
66

67
	void zero() ;
68

69
70
71
72
73
	void operator= (const Quadric<REAL>& q) ;
	Quadric& operator+= (const Quadric<REAL>& q) ;
	Quadric& operator -= (const Quadric<REAL>& q) ;
	Quadric& operator *= (REAL v) ;
	Quadric& operator /= (REAL v) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
74

75
76
	REAL operator() (const VEC4& v) const ;
	REAL operator() (const VEC3& v) const ;
Pierre Kraemer's avatar
Pierre Kraemer committed
77
78
79
80
81

	friend std::ostream& operator<<(std::ostream& out, const Quadric<REAL>& q)
	{
		out << q.A ;
		return out ;
82
	} ;
Pierre Kraemer's avatar
Pierre Kraemer committed
83
84
85
86
87

	friend std::istream& operator>>(std::istream& in, Quadric<REAL>& q)
	{
		in >> q.A ;
		return in ;
88
	} ;
Pierre Kraemer's avatar
Pierre Kraemer committed
89

90
	bool findOptimizedPos(VEC3& v) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
91
92
93
94

private:
	MATRIX44 A ;

95
	REAL evaluate(const VEC4& v) const ;
Pierre Kraemer's avatar
Pierre Kraemer committed
96

97
	bool optimize(VEC4& v) const ;
Pierre Kraemer's avatar
Pierre Kraemer committed
98
99
} ;

100
101
102
103
104
105
106
107
108
template <typename REAL, unsigned int N>
class QuadricNd
{
public:
	static std::string CGoGNnameOfType() { return "QuadricNd" ; }

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

109
	QuadricNd() ;
110

111
	QuadricNd(int i) ;
112

113
	QuadricNd(const VECN& p1_r, const VECN& p2_r, const VECN& p3_r) ;
114

115
	void zero() ;
116

117
	void operator= (const QuadricNd<REAL,N>& q) ;
118

119
	QuadricNd& operator+= (const QuadricNd<REAL,N>& q) ;
120

121
	QuadricNd& operator -= (const QuadricNd<REAL,N>& q) ;
122

123
	QuadricNd& operator *= (REAL v) ;
124

125
	QuadricNd& operator /= (REAL v) ;
126

127
	REAL operator() (const VECNp& v) const ;
128

129
	REAL operator() (const VECN& v) const ;
130
131
132

	friend std::ostream& operator<<(std::ostream& out, const QuadricNd<REAL,N>& q)
	{
133
		out << "(" << q.A << ", " << q.b << ", " << q.c << ")" ;
134
		return out ;
135
	} ;
136
137
138

	friend std::istream& operator>>(std::istream& in, QuadricNd<REAL,N>& q)
	{
139
140
141
		in >> q.A ;
		in >> q.b ;
		in >> q.c ;
142
		return in ;
143
	} ;
144

145
	bool findOptimizedVec(VECN& v) ;
146
147

private:
148
149
150
151
	// Double computation is crucial for stability
	Geom::Matrix<N,N,double> A ;
	Geom::Vector<N,double> b ;
	double c ;
152

153
	REAL evaluate(const VECN& v) const ;
154

155
	bool optimize(VECN& v) const ;
156
157
} ;

158
159
160
161
162
163
template <typename REAL>
class QuadricHF
{
public:
	static std::string CGoGNnameOfType() { return "QuadricHF" ; }

164
	typedef Geom::Vector<3,REAL> VEC3 ;
165

166
	QuadricHF() ;
167

168
	QuadricHF(int i) ;
169

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

172
	QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REAL& alpha) ;
173

174
	~QuadricHF() ;
175

176
177
178
179
180
181
182
183
184
185
186
187
188
189
	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 ;
190
191
192

	friend std::ostream& operator<<(std::ostream& out, const QuadricHF<REAL>& q)
	{
193
		// TODO out << "(" << q.m_A << ", " << q.m_b << ", " << q.m_c << ")" ;
194
		return out ;
195
	} ;
196
197
198

	friend std::istream& operator>>(std::istream& in, QuadricHF<REAL>& q)
	{
199
		// TODO
200
201
202
		//		in >> q.A ;
		//		in >> q.b ;
		//		in >> q.c ;
203
		return in ;
204
	} ;
205

206
	bool findOptimizedCoefs(std::vector<VEC3>& coefs) ;
207
208
209

private:
	// Double computation is crucial for stability
210
	Eigen::MatrixXd m_A ;
211
212
	Eigen::VectorXd m_b[3] ;
	double m_c[3] ;
213

214
	REAL evaluate(const std::vector<VEC3>& coefs) const ;
215

216
	bool optimize(std::vector<VEC3>& coefs) const ;
217

218
	Geom::Matrix33d rotateMatrix(const REAL& gamma) ;
219

220
	Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ;
221

222
223
	Eigen::MatrixXd buildIntegralMatrix_A(const REAL& alpha, unsigned int size) ;
	Eigen::MatrixXd buildIntegralMatrix_B(const REAL& alpha, unsigned int size) ;
224

225
	Eigen::MatrixXd buildIntegralMatrix_C(const REAL& alpha, unsigned int size) ;
226

227
	void fillTensor(Geom::Tensor3d& T) ;
228

229
230
231
	Geom::Tensor3d* tensorsFromCoefs(const std::vector<VEC3>& coefs) ;

	std::vector<VEC3> coefsFromTensors(Geom::Tensor3d* A) ;
232
233
} ;

234
} // Utils
235
236
237

} // CGOGN

238
239
#include "qem.hpp"

Pierre Kraemer's avatar
Pierre Kraemer committed
240
#endif