lightfieldApproximator.hpp 15.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
/*******************************************************************************
* 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                                        *
*                                                                              *
*******************************************************************************/

namespace CGoGN
{

namespace Algo
{

namespace Decimation
{

/************************************************************************************
35
36
 *                         FRAME INTERPOLATION APPROXIMATOR                         *
 ************************************************************************************/
37
38

template <typename PFP>
39
void Approximator_FrameInterpolation<PFP>::approximate(Dart d)
40
{
41
42
43
	const Dart dd = this->m_map.phi2(d) ;

	// if QEM placed new vertex at one of the previous : place same frame
44
45
46
47
48
49
50
51
52
53
54
55
56
//	if (this->m_approxposition[d] == this->m_position[d]) // new Position is position of vertex d
//	{
//		for (unsigned int i = 0 ; i < 3 ; ++i)
//			this->m_approx[i][d] = this->m_attrV[i]->operator[](d) ;
//		//std::cout << "fallback to p1" << std::endl ;
//	}
//	else if (this->m_approxposition[d] == this->m_position[dd]) // new Position is position of vertex dd
//	{
//		for (unsigned int i = 0 ; i < 3 ; ++i)
//			this->m_approx[i][d] = this->m_attrV[i]->operator[](dd) ;
//		// std::cout << "fallback to p2" << std::endl ;
//	}
//	else
57
58
59
60
61
62
63
64
65
	{
		// Create two segments : v0-v1 and v0-v
		VEC3 v0v1 = this->m_position[dd] ;
		v0v1 -= this->m_position[d] ;

		VEC3 v0v = this->m_approxposition[d] ;
		v0v -= this->m_position[d] ;

		// Orthogonal projection of v0-v onto v0-v1 : get coefficient t
66
		REAL t = this->m_position[d] == this->m_position[dd] ? 1. : (v0v1 * v0v) / v0v1.norm() ;
67
68
		t = std::max (std::min (t , REAL(1)) , REAL(0) ) ; // clamp it to [0,1]

69
70
71
72
73
		VEC3& normal1 = this->m_attrV[2]->operator[](d) ;
		VEC3& normal2 = this->m_attrV[2]->operator[](dd) ;

		//Geom::Vector<3,double> n1 = normal1 ;
		//Geom::Vector<3,double> n2 = normal2 ;
74
75
76

		VEC3 newN = slerp(normal1,normal2,t) ; // spherical interpolation
		newN.normalize() ;
77
		// assert (!isnan(newN[0])) ;
78

79
80
81
82
83
		VEC3 newT ;
		if (normal2 != normal1)
			newT = normal2 ^ normal1 ; // i is perpendicular to newNormal
		else
			newT = normal1 ^ VEC3(1,2,3) ; // second random vector
84
85
86
87
88
89
90
91
92
93
		newT.normalize() ;

		VEC3 newB = newN ^ newT ; // (i,j,n) is a direct frame
		newB.normalize() ;

		// save into m_approx
		this->m_approx[0][d] = newT ;
		this->m_approx[1][d] = newB ;
		this->m_approx[2][d] = newN ;
	}
94
95
}

96
97
98
/************************************************************************************
 *                      HEMIFUNCTION COEFS APPROXIMATOR                             *
 ************************************************************************************/
99
100

template <typename PFP>
101
bool Approximator_HemiFuncCoefs<PFP>::init()
102
{
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
	// get frames
	m_frameT = this->m_map.template getAttribute<VEC3, VERTEX>("frameT") ;
	m_frameB = this->m_map.template getAttribute<VEC3, VERTEX>("frameB") ;
	m_frameN = this->m_map.template getAttribute<VEC3, VERTEX>("frameN") ;

	assert((m_frameT.isValid() && m_frameB.isValid() && m_frameN.isValid()) || !"Frame embeddings are not set") ;

	m_newFrameT = this->m_map.template getAttribute<VEC3, EDGE>("approx_frameT") ;
	m_newFrameB = this->m_map.template getAttribute<VEC3, EDGE>("approx_frameB") ;
	m_newFrameN = this->m_map.template getAttribute<VEC3, EDGE>("approx_frameN") ;

	assert((m_newFrameT.isValid() && m_newFrameB.isValid() && m_frameN.isValid())
			|| !"New frame embeddings are not set") ;

	if(this->m_predictor)
	{
		return false ;
	}

	return m_frameT.isValid() && m_frameB.isValid() && m_frameN.isValid() && m_newFrameT.isValid() && m_newFrameB.isValid() && m_newFrameN.isValid() && (m_nbCoefs > 0) ;
123
}
124
125
126
127
128
129

template <typename PFP>
void Approximator_HemiFuncCoefs<PFP>::approximate(Dart d)
{
	const Dart& dd = this->m_map.phi1(d) ;

130
131
132
133
	// get coefs of v1 and v2
	std::vector<VEC3> coefs1, coefs2 ;
	coefs1.resize(m_nbCoefs) ; coefs2.resize(m_nbCoefs) ;
	for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
134
	{
135
136
		coefs1[i] = m_coefs[i]->operator[](d) ;
		coefs2[i] = m_coefs[i]->operator[](dd) ;
137
	}
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

	const VEC3& T = m_newFrameT[d] ;
	const VEC3& T1 = m_frameT[d] ;
	const VEC3& T2 = m_frameT[dd] ;
	const VEC3& B1 = m_frameB[d] ;
	const VEC3& B2 = m_frameB[dd] ;
	const VEC3& N = m_newFrameN[d] ;
	const VEC3& N1 = m_frameN[d] ;
	const VEC3& N2 = m_frameN[dd] ;

	//assert(abs(T * N1) < 1e-6 || !"Approximator_FrameInterpolation<PFP>::approximate: T is not located in the first tangent plane") ;
	//assert(abs(T * N2) < 1e-6 || !"Approximator_FrameInterpolation<PFP>::approximate: T is not located in the second tangent plane") ;

	// Compute D1' and D2'
	VEC3 B1prime = N1 ^ T ;
	B1prime.normalize() ;
	VEC3 B2prime = N2 ^ T ;
	B2prime.normalize() ;

	// Rotation dans sens trigo dans le plan tangent autour de N (T1 --> T)
	const REAL gamma1 = ((B1 * T) > 0 ? 1 : -1) * acos( std::max(std::min(1.0f, T1 * T ), -1.0f)) ; // angle positif ssi
	const REAL gamma2 = ((B2 * T) > 0 ? 1 : -1) * acos( std::max(std::min(1.0f, T2 * T ), -1.0f)) ; // -PI/2 < angle(i,j1) < PI/2  ssi i*j1 > 0
	// Rotation dans le sens trigo autour de l'axe T (N1 --> N)
	const REAL alpha1 = ((N * B1prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N1), -1.0f) ) ; // angle positif ssi
	const REAL alpha2 = ((N * B2prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N2), -1.0f) ) ; // PI/2 < angle(j1',n) < -PI/2 ssi j1'*n < 0

	assert (m_quadricHF.isValid() || !"LightfieldApproximator:approximate quadricHF is not valid") ;

166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
//	if (T == T1)
//	{
//		// std::cout << "Fallback to p1" << std::endl ;
//
//		m_quadricHF[d] = Utils::QuadricHF<REAL>(coefs1, gamma1, 0) ;
//		m_quadricHF[d] += Utils::QuadricHF<REAL>(coefs2, gamma2, alpha2 - alpha1) ;
//
//		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
//			this->m_approx[i][d] = m_coefs[i]->operator[](d) ;
//	}
//	else if (T == T2)
//	{
//		// std::cout << "Fallback to p2" << std::endl ;
//
//		m_quadricHF[d] = Utils::QuadricHF<REAL>(coefs1, gamma1, alpha1 - alpha2) ;
//		m_quadricHF[d] += Utils::QuadricHF<REAL>(coefs2, gamma2, 0) ;
//
//		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
//			this->m_approx[i][d] = m_coefs[i]->operator[](dd) ;
//	}
//	else
187
	{
188
189
190
191
192
193
194
195
196
197
198
199
200
//		std::cout << "N1: " << N1 << std::endl ;
//		std::cout << "N2: " << N2 << std::endl ;
//		std::cout << "N: " << N << std::endl ;
//		std::cout << "T: " << T << std::endl ;
//		std::cout << "T1: " << T1 << std::endl ;
//		std::cout << "T2: " << T2 << std::endl ;
//		std::cout << "gamma1 " << gamma1 << std::endl ;
//		std::cout << "gamma2 " << gamma2 << std::endl ;
//		std::cout << "alpha1 " << alpha1 << std::endl ;
//		std::cout << "alpha2 " << alpha2 << std::endl ;
		// Create and sum quadrics
		m_quadricHF[d] = Utils::QuadricHF<REAL>(coefs1, gamma1, alpha1) ;
		m_quadricHF[d] += Utils::QuadricHF<REAL>(coefs2, gamma2, alpha2) ;
201

202
203
204
205
		std::vector<VEC3> coefs ;
		// Compute new function
		bool opt = m_quadricHF[d].findOptimizedCoefs(coefs) ;
		// copy result
206
		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
207
			this->m_approx[i][d] = coefs[i] ;
208
209

		if (!opt)
Kenneth Vanhoey's avatar
Kenneth Vanhoey committed
210
			std::cerr << "LightfieldApproximator::Inversion failed: coefs are averaged" << std::endl ;
211
	}
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
}


/************************************************************************************
 *           FRAME INTERPOLATION APPROXIMATOR FOR HALF EDGE CONTRACTION             *
 ************************************************************************************/

template <typename PFP>
void Approximator_FrameInterpolationHalfEdge<PFP>::approximate(Dart d)
{
	const Dart dd = this->m_map.phi2(d) ;

	// if QEM placed new vertex at one of the previous : place same frame
//	if (this->m_approxposition[d] == this->m_position[d]) // new Position is position of vertex d
//	{
//		for (unsigned int i = 0 ; i < 3 ; ++i)
//			this->m_approx[i][d] = this->m_attrV[i]->operator[](d) ;
//		//std::cout << "fallback to p1" << std::endl ;
//	}
//	else if (this->m_approxposition[d] == this->m_position[dd]) // new Position is position of vertex dd
//	{
//		for (unsigned int i = 0 ; i < 3 ; ++i)
//			this->m_approx[i][d] = this->m_attrV[i]->operator[](dd) ;
//		// std::cout << "fallback to p2" << std::endl ;
//	}
//	else
238
	{
239
240
241
		// Create two segments : v0-v1 and v0-v
		VEC3 v0v1 = this->m_position[dd] ;
		v0v1 -= this->m_position[d] ;
242

243
244
		VEC3 v0v = this->m_approxposition[d] ;
		v0v -= this->m_position[d] ;
245

246
247
248
249
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
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
		// Orthogonal projection of v0-v onto v0-v1 : get coefficient t
		REAL t = this->m_position[d] == this->m_position[dd] ? 1. : (v0v1 * v0v) / v0v1.norm() ;
		t = std::max (std::min (t , REAL(1)) , REAL(0) ) ; // clamp it to [0,1]

		VEC3& normal1 = this->m_attrV[2]->operator[](d) ;
		VEC3& normal2 = this->m_attrV[2]->operator[](dd) ;

		// Keep normal
		VEC3 newN = normal1 ;
		// assert (!isnan(newN[0])) ;

		// but do the rotation in the tangent plane
		VEC3 newT ;
		if (normal2 != normal1)
			newT = normal2 ^ normal1 ; // i is perpendicular to newNormal
		else
			newT = normal1 ^ VEC3(1,2,3) ; // second random vector
		newT.normalize() ;

		VEC3 newB = newN ^ newT ; // (i,j,n) is a direct frame
		newB.normalize() ;

		// save into m_approx
		this->m_approx[0][d] = newT ;
		this->m_approx[1][d] = newB ;
		this->m_approx[2][d] = newN ;
	}
}

/************************************************************************************
 *            HEMIFUNCTION COEFS APPROXIMATOR FOR HALF EDGE CONTRACTION             *
 ************************************************************************************/

template <typename PFP>
bool Approximator_HemiFuncCoefsHalfEdge<PFP>::init()
{
	// get frames
	m_frameT = this->m_map.template getAttribute<VEC3, VERTEX>("frameT") ;
	m_frameB = this->m_map.template getAttribute<VEC3, VERTEX>("frameB") ;
	m_frameN = this->m_map.template getAttribute<VEC3, VERTEX>("frameN") ;

	assert((m_frameT.isValid() && m_frameB.isValid() && m_frameN.isValid()) || !"Frame embeddings are not set") ;

	m_newFrameT = this->m_map.template getAttribute<VEC3, DART>("approx_frameT") ;
	m_newFrameB = this->m_map.template getAttribute<VEC3, DART>("approx_frameB") ;
	m_newFrameN = this->m_map.template getAttribute<VEC3, DART>("approx_frameN") ;

	assert((m_newFrameT.isValid() && m_newFrameB.isValid() && m_frameN.isValid())
			|| !"New frame embeddings are not set") ;

	if(this->m_predictor)
	{
		return false ;
299
	}
300
301
302
303
304
305
306
307
308
309
310
311
312

	return m_frameT.isValid() && m_frameB.isValid() && m_frameN.isValid() && m_newFrameT.isValid() && m_newFrameB.isValid() && m_newFrameN.isValid() && (m_nbCoefs > 0) ;
}

template <typename PFP>
void Approximator_HemiFuncCoefsHalfEdge<PFP>::approximate(Dart d)
{
	const Dart& dd = this->m_map.phi1(d) ;

	// get coefs of v1 and v2
	std::vector<VEC3> coefs1, coefs2 ;
	coefs1.resize(m_nbCoefs) ; coefs2.resize(m_nbCoefs) ;
	for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
313
	{
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
		coefs1[i] = m_coefs[i]->operator[](d) ;
		coefs2[i] = m_coefs[i]->operator[](dd) ;
	}

	const VEC3& T = m_newFrameT[d] ;
	const VEC3& T1 = m_frameT[d] ;
	const VEC3& T2 = m_frameT[dd] ;
	const VEC3& B1 = m_frameB[d] ;
	const VEC3& B2 = m_frameB[dd] ;
	const VEC3& N = m_newFrameN[d] ;
	const VEC3& N1 = m_frameN[d] ;
	const VEC3& N2 = m_frameN[dd] ;

	//assert(abs(T * N1) < 1e-6 || !"Approximator_FrameInterpolation<PFP>::approximate: T is not located in the first tangent plane") ;
	//assert(abs(T * N2) < 1e-6 || !"Approximator_FrameInterpolation<PFP>::approximate: T is not located in the second tangent plane") ;

	// Compute D1' and D2'
	VEC3 B1prime = N1 ^ T ;
	B1prime.normalize() ;
	VEC3 B2prime = N2 ^ T ;
	B2prime.normalize() ;

	// Rotation dans sens trigo dans le plan tangent autour de N (T1 --> T)
	const REAL gamma1 = ((B1 * T) > 0 ? 1 : -1) * acos( std::max(std::min(1.0f, T1 * T ), -1.0f)) ; // angle positif ssi
	const REAL gamma2 = ((B2 * T) > 0 ? 1 : -1) * acos( std::max(std::min(1.0f, T2 * T ), -1.0f)) ; // -PI/2 < angle(i,j1) < PI/2  ssi i*j1 > 0
	// Rotation dans le sens trigo autour de l'axe T (N1 --> N)
Kenneth Vanhoey's avatar
Kenneth Vanhoey committed
340
341
	const REAL alpha1 = (N == N1) ? 0 : ((N * B1prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N1), -1.0f) ) ; // angle positif ssi
	const REAL alpha2 = (N == N2) ? 0 : ((N * B2prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N2), -1.0f) ) ; // PI/2 < angle(j1',n) < -PI/2 ssi j1'*n < 0
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379

	assert (m_quadricHF.isValid() || !"LightfieldApproximator:approximate quadricHF is not valid") ;

//	if (T == T1)
//	{
//		// std::cout << "Fallback to p1" << std::endl ;
//
//		m_quadricHF[d] = Utils::QuadricHF<REAL>(coefs1, gamma1, 0) ;
//		m_quadricHF[d] += Utils::QuadricHF<REAL>(coefs2, gamma2, alpha2 - alpha1) ;
//
//		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
//			this->m_approx[i][d] = m_coefs[i]->operator[](d) ;
//	}
//	else if (T == T2)
//	{
//		// std::cout << "Fallback to p2" << std::endl ;
//
//		m_quadricHF[d] = Utils::QuadricHF<REAL>(coefs1, gamma1, alpha1 - alpha2) ;
//		m_quadricHF[d] += Utils::QuadricHF<REAL>(coefs2, gamma2, 0) ;
//
//		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
//			this->m_approx[i][d] = m_coefs[i]->operator[](dd) ;
//	}
//	else
	{
//		std::cout << "N1: " << N1 << std::endl ;
//		std::cout << "N2: " << N2 << std::endl ;
//		std::cout << "N: " << N << std::endl ;
//		std::cout << "T: " << T << std::endl ;
//		std::cout << "T1: " << T1 << std::endl ;
//		std::cout << "T2: " << T2 << std::endl ;
//		std::cout << "gamma1 " << gamma1 << std::endl ;
//		std::cout << "gamma2 " << gamma2 << std::endl ;
//		std::cout << "alpha1 " << alpha1 << std::endl ;
//		std::cout << "alpha2 " << alpha2 << std::endl ;

		// Add quadric for first one only (for half-edge contraction)
		m_quadricHF[d] = Utils::QuadricHF<REAL>(coefs1, gamma1, alpha1) ;
380

381
382
383
384
385
		std::vector<VEC3> coefs ;
		// Compute new function
		bool opt = m_quadricHF[d].findOptimizedCoefs(coefs) ;
		// copy result
		for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
Kenneth Vanhoey's avatar
Kenneth Vanhoey committed
386
			this->m_approx[i][d] = coefs[i] ;
387
388

		if (!opt)
Kenneth Vanhoey's avatar
Kenneth Vanhoey committed
389
390
391
392
		{
			std::cerr << "LightfieldApproximatorHalfCollapse::Optimization failed (should never happen since no optim is done)" << std::endl ;
			std::cout << alpha1 << std::endl ;
		}
393
394
395

		// Add second one for error evaluation
		m_quadricHF[d] += Utils::QuadricHF<REAL>(coefs2, gamma2, alpha2) ;
396
	}
397
398
}

399
400
401
402
403
} //namespace Decimation

} //namespace Algo

} //namespace CGoGN