Assignments Assignments

geometry.inl
Go to the documentation of this file.
1/*
2 Copyright (c) 2019, Michael Kazhdan
3 All rights reserved.
4
5 Redistribution and use in source and binary forms, with or without modification,
6 are permitted provided that the following conditions are met:
7
8 Redistributions of source code must retain the above copyright notice, this list of
9 conditions and the following disclaimer. Redistributions in binary form must reproduce
10 the above copyright notice, this list of conditions and the following disclaimer
11 in the documentation and/or other materials provided with the distribution.
12
13 Neither the name of the Johns Hopkins University nor the names of its contributors
14 may be used to endorse or promote products derived from this software without specific
15 prior written permission.
16
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
18 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
19 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
20 SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
21 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
22 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
23 BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
25 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
26 DAMAGE.
27*/
28
29#include <algorithm>
30#include <SVD/SVDFit.h>
31#include <SVD/MatrixMNTC.h>
32#include <Util/exceptions.h>
33
34 namespace Util
35{
37 // Point //
39 template< unsigned int Dim >
40 void Point< Dim >::_init( const double *values , unsigned int sz )
41 {
42 if ( sz==0 ) memset( _p , 0 , sizeof(_p) );
43 else if( sz==Dim ) memcpy( _p , values , sizeof(_p) );
44 else ERROR_OUT( "Should never be called" );
45 }
46
47 template< unsigned int Dim >
48 Point< Dim >::Point( void ){ memset( _p , 0 , sizeof(_p) ); }
49
50 template< unsigned int Dim >
51 Point< Dim >::Point( const Point &p ){ memcpy( _p , p._p , sizeof(_p) ); }
52
53 template< unsigned int Dim >
54 template< typename ... Doubles >
55 Point< Dim >::Point( Doubles ... values )
56 {
57 static_assert( sizeof...(values)==Dim || sizeof...(values)==0 , "[ERROR] Point< Dim >::Point: Invalid number of coefficients" );
58 const double _values[] = { static_cast< double >( values )... };
59 _init( _values , sizeof...(values) );
60 }
61
62 template< unsigned int Dim >
63 Point< Dim > Point< Dim >::operator * ( double s ) const { Point p ; for( int i=0 ; i<Dim ; i++ ) p._p[i] = _p[i] * s ; return p; }
64
65 template< unsigned int Dim >
66 Point< Dim > Point< Dim >::operator + ( const Point &p ) const { Point q ; for( int i=0 ; i<Dim ; i++ ) q._p[i] += _p[i] + p._p[i] ; return q; }
67
68 template< unsigned int Dim >
69 double Point< Dim >::dot( const Point &q ) const
70 {
71 double dot = 0;
72 for( int i=0 ; i<Dim ; i++ ) dot += _p[i] * q._p[i];
73 return dot;
74 }
75
76 template< unsigned int Dim >
77 double& Point< Dim >::operator[] ( int i ){ return _p[i]; }
78
79 template< unsigned int Dim >
80 const double &Point< Dim >::operator[] ( int i ) const { return _p[i]; }
81
82 template< unsigned int Dim >
83 Point< Dim > Point< Dim >::operator * ( const Point &q ) const
84 {
85 Point p;
86 for( int i=0 ; i<Dim ; i++ ) p[i] = _p[i]*q._p[i];
87 return p;
88 }
89
90 template< unsigned int Dim >
91 Point< Dim > Point< Dim >::operator / ( const Point &q ) const
92 {
93 Point p;
94 for( int i=0 ; i<Dim ; i++ ) p[i] = _p[i]/q._p[i];
95 return p;
96 }
97
98 template< unsigned int Dim >
99 Point< Dim > &Point< Dim >::operator *= ( const Point &q ){ return (*this) = (*this) * q; }
100
101 template< unsigned int Dim >
102 Point< Dim > &Point< Dim >::operator /= ( const Point &q ){ return (*this) = (*this) / q; }
103
104 template< unsigned int Dim >
105 template< typename ... Points >
106 Point< Dim > Point< Dim >::CrossProduct( Points ... points )
107 {
108 static_assert( sizeof ... ( points )==Dim-1 , "[ERROR] Number of points in cross-product must be one less than the dimension" );
109 const Point< Dim > _points[] = { points ... };
110 return CrossProduct( _points );
111 }
112
113 template< unsigned int Dim >
114 Point< Dim > Point< Dim >::CrossProduct( Point *points ){ return CrossProduct( (const Point *)points );}
115
116 template< unsigned int Dim >
117 Point< Dim > Point< Dim >::CrossProduct( const Point *points )
118 {
119 Matrix< Dim , Dim > M;
120 for( int d=0 ; d<Dim ; d++ ) for( int c=0 ; c<Dim-1 ; c++ ) M(d,c) = points[c][d];
121 Point p;
122 for( int d=0 ; d<Dim ; d++ ) p[d] = ( d&1 ) ? -M.subDeterminant( d , Dim-1 ) : M.subDeterminant( d , Dim-1 );
123 return p;
124 }
125
126 template< unsigned int Dim >
127 std::ostream &operator << ( std::ostream &stream , const Point< Dim > &p )
128 {
129 for( int i=0 ; i<Dim-1 ; i++ ) stream << p[i] << " ";
130 stream << p[Dim-1];
131 return stream;
132 }
133 template< unsigned int Dim >
134 std::istream &operator >> ( std::istream &stream , Point< Dim > &p )
135 {
136 for( int i=0 ; i<Dim ; i++ ) stream >> p[i];
137 return stream;
138 }
139
140 /////////////////
141 // _BaseMatrix //
143 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
144 double _BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType >::dot( const _BaseMatrix &m ) const
145 {
146 double dot = 0;
147 for( int r=0 ; r<Rows ; r++ ) for( int c=0 ; c<Cols ; c++ ) dot += operator()(r,c) * m(r,c);
148 return dot;
149 }
150
151 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
152 _BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType >::_BaseMatrix( void ){ memset( _m , 0 , sizeof(_m) ); }
153
154 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
155 double& _BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType >::operator() ( int r , int c ) { return _m[r][c]; }
156
157 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
158 const double &_BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType >::operator() ( int r , int c ) const { return _m[r][c]; }
159
160 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
161 MatrixTransposeType _BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType >::transpose( void ) const
162 {
163 MatrixTransposeType n;
164 for( int r=0 ; r<Rows ; r++ ) for( int c=0 ; c<Cols ; c++ ) n(c,r) = operator()(r,c);
165 return n;
166 }
167
168 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
169 Point< Rows > _BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType >::operator * ( const Point< Cols > &p ) const
170 {
171 Point< Rows > q;
172 for( int r=0 ; r<Rows ; r++ ) for( int c=0 ; c<Cols ; c++ ) q[r] += operator()(r,c) * p[c];
173 return q;
174 }
175
176 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
177 MatrixType _BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType >::operator * ( double s ) const { MatrixType n ; for( int r=0 ; r<Rows ; r++ ) for( int c=0 ; c<Cols ; c++ ) n(r,c) = operator()(r,c) * s ; return n; }
178
179 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
180 MatrixType _BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType >::operator + ( const MatrixType &m ) const { MatrixType n ; for( int r=0 ; r<Rows ; r++ ) for( int c=0 ; c<Cols ; c++ ) n(r,c) = operator()(r,c) + m(r,c) ; return n; }
181
183 // Matrix //
185 template< unsigned int Rows , unsigned int Cols > Matrix< Rows , Cols >::Matrix( void ) : _BaseMatrix< Rows , Cols , Matrix< Rows , Cols > , Matrix< Cols , Rows > >() {}
186
187 template< unsigned int Rows , unsigned int Cols >
188 template< unsigned int _Cols >
189 Matrix< Rows , _Cols > Matrix< Rows , Cols >::operator * ( const Matrix< Cols , _Cols > &m ) const
190 {
191 Matrix< Rows , _Cols > n;
192 for( int r=0 ; r<Rows ; r++ ) for( int c=0 ; c<_Cols ; c++ ) for( int i=0 ; i<Cols ; i++ ) n(r,c) += operator()(r,i) * m(i,c);
193 return n;
194 }
195
197 // SquareMatrix //
199 template< unsigned int Dim > SquareMatrix< Dim >::Matrix( void ) : _BaseMatrix< Dim , Dim , SquareMatrix< Dim > , SquareMatrix< Dim > >() {}
200 template< unsigned int Dim >
201 SquareMatrix< Dim >::Matrix( const SquareMatrix< Dim+1 > &n ){ for( int i=0 ; i<Dim ; i++ ) for( int j=0 ; j<Dim ; j++ ) operator()(i,j) = n(i,j); }
202
203 template< unsigned int Dim >
204 SquareMatrix< Dim >::Matrix( const SquareMatrix< Dim-1 > &n , Point< Dim-1 > p ) : Matrix()
205 {
206 for( int i=0 ; i<Dim-1 ; i++ ) for( int j=0 ; j<Dim-1 ; j++ ) operator()(i,j) = n(i,j);
207 operator()(Dim-1,Dim-1) = 1.;
208 for( int i=0 ; i<Dim-1 ; i++ ) operator()(i,Dim-1) = p[i];
209 }
210
211 template< unsigned int Dim >
212 template< unsigned int Cols >
213 Matrix< Dim , Cols > SquareMatrix< Dim >::operator * ( const Matrix< Dim , Cols > &m ) const
214 {
215 Matrix< Dim , Cols > n;
216 for( int r=0 ; r<Dim ; r++ ) for( int c=0 ; c<Cols ; c++ ) for( int i=0 ; i<Dim ; i++ ) n(r,c) += operator()(r,i) * m(i,c);
217 return n;
218 }
219
220 template< unsigned int Dim >
221 double SquareMatrix< Dim >::subDeterminant( int r , int c ) const
222 {
223 SquareMatrix< Dim-1 > m;
224 int rr[Dim-1] , cc[Dim-1];
225 for( int a=0 , _r=0 , _c=0 ; a<Dim ; a++ )
226 {
227 if( a!=r ) rr[_r++] = a;
228 if( a!=c ) cc[_c++] = a;
229 }
230 for( int _c=0 ; _c<Dim-1 ; _c++ ) for( int _r=0 ; _r<Dim-1 ; _r++ ) m(_r,_c) = operator()( rr[_r] , cc[_c] );
231 return m.determinant();
232 }
233
234 template< unsigned int Dim >
235 double SquareMatrix< Dim >::determinant( void ) const
236 {
237 double det = 0.;
238 for( int d=0 ; d<Dim ; d++ )
239 if( d&1 ) det -= operator()(0,d) * subDeterminant( 0 , d );
240 else det += operator()(0,d) * subDeterminant( 0 , d );
241 return det;
242 }
243
244 template<>
245 inline double SquareMatrix< 3 >::determinant( void ) const
246 {
247 return
248 operator()(0,0)*( operator()(1,1)*operator()(2,2) - operator()(1,2)*operator()(2,1) ) -
249 operator()(0,1)*( operator()(1,0)*operator()(2,2) - operator()(1,2)*operator()(2,0) ) +
250 operator()(0,2)*( operator()(1,0)*operator()(2,1) - operator()(1,1)*operator()(2,0) );
251 }
252
253 template<>
254 inline double SquareMatrix< 2 >::determinant( void ) const { return operator()(0,0)*operator()(1,1) - operator()(0,1)*operator()(1,0); }
255
256 template<>
257 inline double SquareMatrix< 1 >::determinant( void ) const { return operator()(0,0); }
258
259 template< unsigned int Dim >
260 double SquareMatrix< Dim >::trace( void ) const
261 {
262 double tr = 0;
263 for( int i=0 ; i<Dim ; i++ ) tr += operator()(i,i);
264 return tr;
265 }
266 template< unsigned int Dim >
267 SquareMatrix< Dim > SquareMatrix< Dim >::inverse( void ) const
268 {
269 Matrix inv;
270 if( !setInverse( inv ) ) THROW( " singular matrix" );
271 return inv;
272 }
273
274 template< unsigned int Dim >
275 bool SquareMatrix< Dim >::setInverse( Matrix &inv ) const
276 {
277 double d = determinant();
278 if( !d ) return false;
279 for( int i=0 ; i<Dim ; i++ ) for( int j=0 ; j<Dim ; j++ )
280 if( (i+j)%2==0 ) inv(i,j) = subDeterminant( j , i ) / d;
281 else inv(i,j) = -subDeterminant( j , i ) / d;
282 return true;
283 }
284
285 template<>
286 inline bool SquareMatrix< 3 >::setInverse( Matrix &inv ) const
287 {
288 double det = determinant();
289 if( !det ) return false;
290 inv(0,0) = ( operator()(1,1)*operator()(2,2) - operator()(1,2)*operator()(2,1) ) / det;
291 inv(0,1) = -( operator()(0,1)*operator()(2,2) - operator()(2,1)*operator()(0,2) ) / det;
292 inv(0,2) = ( operator()(0,1)*operator()(1,2) - operator()(0,2)*operator()(1,1) ) / det;
293 inv(1,0) = -( operator()(1,0)*operator()(2,2) - operator()(1,2)*operator()(2,0) ) / det;
294 inv(1,1) = ( operator()(0,0)*operator()(2,2) - operator()(0,2)*operator()(2,0) ) / det;
295 inv(1,2) = -( operator()(0,0)*operator()(1,2) - operator()(0,2)*operator()(1,0) ) / det;
296 inv(2,0) = ( operator()(1,0)*operator()(2,1) - operator()(1,1)*operator()(2,0) ) / det;
297 inv(2,1) = -( operator()(0,0)*operator()(2,1) - operator()(0,1)*operator()(2,0) ) / det;
298 inv(2,2) = ( operator()(0,0)*operator()(1,1) - operator()(0,1)*operator()(1,0) ) / det;
299
300 return true;
301 }
302
303 template<>
304 inline bool SquareMatrix< 2 >::setInverse( Matrix &inv ) const
305 {
306 double det = determinant();
307 if( !det ) return false;
308 inv(0,0) = operator()(1,1) / det;
309 inv(1,1) = operator()(0,0) / det;
310 inv(1,0) = -operator()(1,0) / det;
311 inv(0,1) = -operator()(0,1) / det;
312
313 return true;
314 }
315
316 template<>
317 inline bool SquareMatrix< 1 >::setInverse( Matrix &inv ) const
318 {
319 double det = operator()(0,0);
320 if( !det ) return false;
321 inv(0,0) = 1./det;
322 return true;
323 }
324
325 template< unsigned int Dim >
326 Point< Dim-1 > SquareMatrix< Dim >::operator * ( const Point< Dim-1 > &p ) const
327 {
328 Point< Dim > q;
329 for( int i=0 ; i<Dim-1 ; i++ ) q[i] = p[i];
330 q[Dim-1] = 1;
331 q = (*this) * q;
332 Point< Dim-1 > _q;
333 for( int i=0 ; i<Dim-1 ; i++ ) _q[i] = q[i] / q[Dim-1];
334 return _q;
335 }
336
337 template< unsigned int Dim >
338 SquareMatrix< Dim > SquareMatrix< Dim >::Identity( void )
339 {
340 Matrix m;
341 for( int i=0 ; i<Dim ; i++ ) m(i,i) = 1;
342 return m;
343 }
344
345 template< unsigned int Dim >
346 void SquareMatrix< Dim >::SVD( Matrix& r1 , Matrix& d , Matrix& r2 ) const
347 {
348 GXMatrixMNd M( Dim , Dim );
349 GXMatrixMNd U, W, Vt;
350
351 d = r1 = r2 = Identity();
352
353 for( int i=0 ; i<Dim ; i++ ) for( int j=0 ; j<Dim ; j++ ) M(i,j) = operator()(i,j);
354 SVDMat( M , U , W , Vt ); // M == U . DiagonalMatrix(W) . Vt
355 for( int i=0 ; i<Dim ; i++ )
356 {
357 for( int j=0 ; j<Dim ; j++ ) r1(j,i) = U(j,i) , r2(j,i) = Vt(j,i);
358 d(i,i) = W(i,0);
359 }
360 }
361
362 // Code borrowed from:
363 // Linear Combination of Transformations
364 // Marc Alexa
365 template< unsigned int Dim >
366 SquareMatrix< Dim > SquareMatrix< Dim >::SquareRoot( const Matrix& m , double eps )
367 {
368 Matrix X,Y;
369 X = m;
370 Y = Identity();
371 while( (X*X-m).squareNorm()>eps*eps )
372 {
373 Matrix iX = X.inverse();
374 Matrix iY = Y.inverse();
375 X = (X+iY) / 2;
376 Y = (Y+iX) / 2;
377 }
378 return X;
379 }
380
381 template< unsigned int Dim >
382 SquareMatrix< Dim > SquareMatrix< Dim >::Log( const Matrix& m , double eps )
383 {
384 Matrix I = Identity();
385 Matrix X , Z , A=m;
386 int k=0;
387 while( (A-I).squareNorm()>0.25 )
388 {
389 A = SquareRoot( A , eps );
390 k++;
391 }
392 A = I-A;
393 X = Z = A;
394 int i = 1;
395 while( Z.squareNorm()>eps*eps )
396 {
397 Z = Z*A;
398 i++;
399 X += Z*( 1.0/i );
400 }
401 return X * ( -pow( 2.0 , (double)k ) );
402 }
403
404 template< unsigned int Dim >
405 SquareMatrix< Dim > SquareMatrix< Dim >::symmetrize( void ) const { return ( (*this)+transpose() ) / 2; }
406
407 template< unsigned int Dim >
408 SquareMatrix< Dim > SquareMatrix< Dim >::skewSymmetrize( void ) const { return ( (*this)-transpose() ) / 2; }
409
410 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
411 std::ostream &operator << ( std::ostream &stream , const _BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType > &m )
412 {
413 for( unsigned int c=0 ; c<Cols ; c++ ) for( unsigned int r=0 ; r<Rows ; r++ )
414 {
415 stream << m(r,c);
416 if( r!=Rows-1 || c!=Cols-1 ) stream << " ";
417 }
418 return stream;
419 }
420
421 template< unsigned int Rows , unsigned int Cols , typename MatrixType , typename MatrixTransposeType >
422 std::istream &operator >> ( std::istream &stream , _BaseMatrix< Rows , Cols , MatrixType , MatrixTransposeType > &m )
423 {
424 for( unsigned int c=0 ; c<Cols ; c++ ) for( unsigned int r=0 ; r<Rows ; r++ ) stream >> m(r,c);
425 return stream;
426 }
427
428 ///////////
429 // Plane //
431 template< unsigned int Dim >
432 Plane< Dim >::Plane( void ) : distance(0) {}
433
434 template< unsigned int Dim >
435 Plane< Dim >::Plane( const Point< Dim > &n , const Point< Dim > &p )
436 {
437 normal = n.unit();
438 distance = Point< Dim >::Dot( normal , p );
439 }
440
441 template< unsigned int Dim >
442 template< typename ... Points >
443 Plane< Dim >::Plane( Points ... points )
444 {
445 static_assert( sizeof ... ( points )==Dim , "[ERROR] Number of points in plane constructor must equal the dimension" );
446 const Point< Dim > _points[] = { points ... };
447 (*this) = Plane( _points );
448 }
449
450 template< unsigned int Dim >
451 Plane< Dim >::Plane( Point< Dim > *points ) : Plane( (const Point< Dim > *)points ) {}
452
453 template< unsigned int Dim >
454 Plane< Dim >::Plane( const Point< Dim > *points )
455 {
456 Point< Dim > _points[Dim-1];
457 for( int i=1 ; i<Dim ; i++ ) _points[i-1] = points[i] - points[0];
458 normal = Point< Dim >::CrossProduct( _points ).unit();
459 distance = Point< Dim >::Dot( normal , points[0] );
460 }
461
462 template< unsigned int Dim >
463 double Plane< Dim >::operator() ( const Point< Dim > &p ) const
464 {
465 return Point< Dim >::Dot( normal , p ) - distance;
466 }
467
469 // Ray //
471 template< unsigned int Dim >
472 Ray< Dim >::Ray( void ){}
473
474 template< unsigned int Dim >
475 Ray< Dim >::Ray( const Point< Dim > &p , const Point< Dim > &d ) : position(p) , direction(d) {}
476
477 template< unsigned int Dim >
478 Point< Dim > Ray< Dim >::operator() ( double s ) const { return position+direction*s; }
479
480 template< unsigned int Dim >
481 Ray< Dim > Ray< Dim >::operator + ( const Point< Dim > &p ) const { return Ray( position+p , direction );}
482
483 template< unsigned int Dim >
484 Ray< Dim > &Ray< Dim >::operator += ( const Point< Dim > &p ){ position += p ; return *this; }
485
486 template< unsigned int Dim >
487 Ray< Dim > Ray< Dim >::operator - ( const Point< Dim > &p ) const { return Ray( position-p , direction );}
488
489 template< unsigned int Dim >
490 Ray< Dim > &Ray< Dim >::operator -= ( const Point< Dim > &p ){ position -= p ; return *this; }
491
492 template< unsigned int Dim >
493 Ray< Dim > operator * ( const Matrix< Dim+1 , Dim+1 > &m , const Ray< Dim >& r )
494 {
495 return Ray< Dim >( m * r.position , Matrix< Dim , Dim >(m) * r.direction );
496 }
497
498
500 // BoundingBox //
502 template< unsigned int Dim >
503 BoundingBox< Dim >::BoundingBox( void ){}
504
505 template< unsigned int Dim >
506 BoundingBox< Dim >::BoundingBox( const Point< Dim > &p1 , const Point< Dim > &p2 )
507 {
508 for( int d=0 ; d<Dim ; d++ ) _p[0][d] = std::min< double >( p1[d] , p2[d] ) , _p[1][d] = std::max< double >( p1[d] , p2[d] );
509 }
510
511 template< unsigned int Dim >
512 BoundingBox< Dim >::BoundingBox( const Point< Dim > *pList , int pSize )
513 {
514 if( pSize>0 )
515 {
516 _p[0] = _p[1] = pList[0];
517 for( int i=1 ; i<pSize ; i++ ) for( int j=0 ; j<Dim ; j++ ) _p[0][j] = std::min< double >( _p[0][j] , pList[i][j] ) , _p[1][j] = std::max< double >( _p[1][j] , pList[i][j] );
518 }
519 }
520
521 template< unsigned int Dim >
522 Point< Dim > &BoundingBox< Dim >::operator[] ( int idx ){ return _p[idx]; }
523
524 template< unsigned int Dim >
525 const Point< Dim > &BoundingBox< Dim >::operator[] ( int idx ) const { return _p[idx]; }
526
527 template< unsigned int Dim >
528 BoundingBox< Dim > BoundingBox< Dim >::operator + ( const BoundingBox &b ) const
529 {
530 Point< Dim > pList[4];
531 Point< Dim > q;
532
533 if( b.isEmpty() ) return *this;
534 if( isEmpty() ) return b;
535 pList[0] = _p[0];
536 pList[1] = _p[1];
537 pList[2] = b._p[0];
538 pList[3] = b._p[1];
539 return BoundingBox( pList , 4 );
540 }
541
542 template< unsigned int Dim >
543 BoundingBox< Dim > &BoundingBox< Dim >::operator += ( const BoundingBox &b ){ return (*this) = (*this) + b; }
544
545 template< unsigned int Dim >
546 BoundingBox< Dim > BoundingBox< Dim >::operator ^ ( const BoundingBox &b ) const
547 {
548
549 if( isEmpty() || b.isEmpty() ) return BoundingBox();
550 BoundingBox _b;
551 for( int j=0 ; j<Dim ; j++ ) _b._p[0][j] = std::max< double >( _p[0][j] , b._p[0][j] ) , _b._p[1][j] = std::min< double >( _p[1][j] , b._p[1][j] );
552 if( _b.isEmpty() ) _b._p[0] = _b._p[1] = ( _b._p[0] + _b._p[1] ) / 2;
553 return _b;
554 }
555
556 template< unsigned int Dim >
557 BoundingBox< Dim > &BoundingBox< Dim >::operator ^= ( const BoundingBox &b ){ return (*this) = (*this) ^ b; }
558
559 template< unsigned int Dim >
560 BoundingBox< Dim > operator * ( const Matrix< Dim+1 , Dim+1 > &m , const BoundingBox< Dim > &b )
561 {
562 Point< Dim > v[1<<Dim];
563 for( int idx=0 ; idx<(1<<Dim) ; idx++ )
564 {
565 Point< Dim > p;
566 for( int d=0 ; d<Dim ; d++ ) p[d] = b[(idx>>d)&1][d];
567 v[idx] = m * p;
568 }
569 return BoundingBox< Dim >( v , 1<<Dim );
570 }
571
572 template< unsigned int Dim >
573 bool BoundingBox< Dim >::isInside( const Point< Dim > &p ) const
574 {
575 for( int d=0 ; d<Dim ; d++ ) if( p[d]<=_p[0][d] || p[d]>=_p[1][d] ) return false;
576 return true;
577 }
578
579 template< unsigned int Dim >
580 bool BoundingBox< Dim >::isEmpty( void ) const
581 {
582 for( int d=0 ; d<Dim ; d++ ) if( _p[0][d]>=_p[1][d] ) return true;
583 return false;
584 }
585
586 template< unsigned int Dim >
587 std::ostream &operator << ( std::ostream &stream , const BoundingBox< Dim > &b )
588 {
589 stream << "[ " << b[0] << " ] [ " << b[1] << " ]";
590 return stream;
591 }
592
594 // Quadric //
596 template< unsigned int Dim >
597 Quadric< Dim >::Quadric( void ) : _C(0) {}
598
599 template< unsigned int Dim >
600 Quadric< Dim >::Quadric( const Matrix< Dim , Dim > &Q , const Point< Dim > &L , const double &C ) :_Q(Q) , _L(L) , _C(C) {}
601
602 template< unsigned int Dim >
603 Quadric< Dim >::Quadric( const Matrix< Dim+1 , Dim+1 > &Q )
604 {
605 _C = Q(Dim,Dim);
606 for( unsigned int i=0 ; i<Dim ; i++ )
607 {
608 _L[i] = ( Q(i,Dim) + Q(Dim,i) ) / 2.;
609 for( unsigned int j=0 ; j<Dim ; j++ ) _Q(i,j) = ( Q(i,j) + Q(j,i) ) / 2.;
610 }
611 }
612
613 template< unsigned int Dim >
614 SquareMatrix< Dim+1 > Quadric< Dim >::operator()( void ) const
615 {
616 SquareMatrix< Dim+1 > Q;
617 Q(Dim,Dim) = _C;
618 for( unsigned int i=0 ; i<Dim ; i++ )
619 {
620 Q(Dim,i) = Q(i,Dim) = _L[i];
621 for( unsigned int j=0 ; j<Dim ; j++ ) Q(i,j) = _Q(i,j);
622 }
623 return Q;
624 }
625
626 template< unsigned int Dim >
627 Matrix< Dim , Dim > Quadric< Dim >::getQuadratic( void ) const { return _Q; }
628
629 template< unsigned int Dim >
630 Point< Dim > Quadric< Dim >::getLinear( void ) const { return _L; }
631
632 template< unsigned int Dim >
633 double Quadric< Dim >::getConstant( void ) const { return _C; }
634
635 template< unsigned int Dim >
636 void Quadric< Dim >::setQuadratic( Matrix< Dim , Dim > Q ){ _Q = ( Q + Q.transpose() )/2.; }
637
638 template< unsigned int Dim >
639 void Quadric< Dim >::setLinear( Point< Dim > L ){ _L = L; }
640
641 template< unsigned int Dim >
642 void Quadric< Dim >::setConstant( double C ){ _C = C; }
643
644 template< unsigned int Dim >
645 double Quadric< Dim >::operator()( const Point< Dim > &p ) const { return Point< Dim >::Dot( p , _Q*p ) + 2. * Point< Dim >::Dot( p , _L ) + _C; }
646
647 template< unsigned int Dim >
648 template< unsigned int _Dim >
649 Quadric< _Dim > Quadric< Dim >::operator * ( const Matrix< Dim , _Dim > &T ) const
650 {
651 // Q(p) = _Q( T(p) )
652 // = [ T * p ]^t * _Q._Q * [ T * p ] + 2 * _L^t * T * p + _C
653 // = p^t * [ T^t *_Q._Q * T ] *p + 2 ( [T^t*_L]^t * p ) + _C
654 Quadric< _Dim > Q;
655 Q._Q = T.transpose() * _Q * T;
656 Q._L = T.transpose() * _L;
657 Q._C = _C;
658 return Q;
659 }
660
661 template< unsigned int Dim >
662 Quadric< Dim > Quadric< Dim >::operator + ( const Point< Dim > &t ) const
663 {
664 // Q(p) = _Q( p+t )
665 // = [ p + t ]^t * _Q._Q * [ p + t ] + 2 * _Q._L^t * [ p + t ] + _Q._C
666 // = p^t * _Q._Q *p + 2 ( [ _Q._Q^t * t + _Q._L ]^t * p ) + _Q._C + t^t * _Q._Q * t + 2 * _Q._L^t * t
667 Quadric< Dim > Q;
668 Q._Q = _Q;
669 Q._L = _Q*t + _L;
670 Q._C = _C + Point< Dim >::Dot( t , _Q*t ) + 2. * Point< Dim >::Dot( _L , t );
671 return Q;
672 }
673
674 template< unsigned int Dim >
675 bool Quadric< Dim >::setExtremum( Point< Dim > &extremum ) const
676 {
677 // The equation for the quardic is given by
678 // Q(p) = p^t * _Q *p + 2 * < _L , p > + c
679 // The gradient is:
680 // \nabla Q|_p = 2 * _Q * p + 2 * _L
681 // Setting the gradient to zero gives:
682 // p = - _Q^{-1} * _L
683 try{ extremum = -_Q.inverse() * _L; }
684 catch( std::exception & ){ return false; }
685 return true;
686 }
687
689 // QuadricBoundingBoxOverlap //
691 template< unsigned int Dim >
692 QuadricBoundingBoxOverlap< Dim >::QuadricBoundingBoxOverlap( void ){}
693
694 template< unsigned int Dim >
695 QuadricBoundingBoxOverlap< Dim >::QuadricBoundingBoxOverlap( const Quadric< Dim > &Q ){ _set( Q ); }
696
697 template< unsigned int Dim >
698 bool QuadricBoundingBoxOverlap< Dim >::operator()( const BoundingBox< Dim > & bBox ) const { return _intersect( _quadric , bBox ); }
699
700 template< unsigned int Dim >
701 void QuadricBoundingBoxOverlap< Dim >::_set( const Quadric< Dim > &Q )
702 {
703 _quadric = Q;
704 _Q = Q._Q;
705 try{ _Qinv = _Q.inverse(); }
706 catch(...){}
707 for( unsigned int d=0 ; d<Dim ; d++ )
708 {
709 unsigned idx = 0;
710 for( unsigned int dd=0 ; dd<Dim ; dd++ ) if( dd!=d ) _boundaryInfo[d]._T(dd,idx++) = 1.;
711 _boundaryInfo[d]._Tt = _boundaryInfo[d]._T.transpose();
712 _boundaryInfo[d]._set( Q * _boundaryInfo[d]._T );
713 _boundaryInfo[d]._Tt_Q = _boundaryInfo[d]._Tt * _Q;
714 }
715 }
716
717 template< unsigned int Dim >
718 bool QuadricBoundingBoxOverlap< Dim >::_intersect( const Quadric< Dim > &Q , const BoundingBox< Dim > &bBox ) const
719 {
720 // Intersection can be found by finding the minimizer of the quadric within the cube and testing if it evaluates to a positive
721 // or negative value.
722 // The minimizer will be:
723 // 1. The point at which the gradient vanishes, if that point is inside the cube. Or,
724 // 2. A point on the face at which the gradient vanishes
725
726 // Case 2:
727 {
728 for( unsigned int d=0 ; d<Dim ; d++ )
729 {
730 BoundingBox< Dim-1 > _bBox;
731
732 // Compute the bounding box and the transformation mapping the (Dim-1)-dimensional plane containing the face into Dim-dimensional space
733
734 // Example:
735 // We have (x,z) -> (x,c,z)
736 // Setting
737 // | 1 0 |
738 // T = | 0 0 |
739 // | 0 1 |
740 // and
741 // | 0 |
742 // v = | c |
743 // | 0 |
744 // We get:
745 // _Q(_p) = Q( T*_p + v )
746 // = _p^t * [ T^t * Q._Q * T ] * _p + 2 * < T*_p + v , Q._L > + 2 * < Q._Q*T*_p , v > + _C + < v , Q._Q * v >
747 // = _p^t * [ T^t * Q._Q * T ] * _p + 2 * < _p , T^t*Q._L > + 2 * < _p , T^t*Q._Q*v > + _C + < v , Q._Q * v > + 2 * < v , Q._L >
748
749 {
750 unsigned idx = 0;
751 for( unsigned int dd=0 ; dd<Dim ; dd++ ) if( dd!=d )
752 {
753 _bBox[0][idx] = bBox[0][dd];
754 _bBox[1][idx] = bBox[1][dd];
755 idx++;
756 }
757 }
758
759 Quadric< Dim-1 > _Q;
760 _Q._Q = _boundaryInfo[d]._Q;
761
762 Point< Dim-1 > _L = _boundaryInfo[d]._Tt * Q._L;
763 Point< Dim-1 > _v;
764 for( unsigned int i=0 ; i<Dim-1 ; i++ ) _v[i] = _boundaryInfo[d]._Tt_Q(i,d);
765
766 for( unsigned int o=0 ; o<2 ; o++ )
767 {
768 double v = bBox[o][d];
769 _Q._L = _L + _v * v;
770 _Q._C = Q._C + v * v * Q._Q(d,d) + 2. * v * Q._L[d];
771 if( _boundaryInfo[d]._intersect( _Q , _bBox ) ) return true;
772 }
773 }
774 }
775
776 // Case 1:
777 Point< Dim > extremum = -_Qinv * Q._L;
778 return bBox.isInside( extremum ) && Q( extremum )<0;
779 }
780
781 inline void QuadricBoundingBoxOverlap< 1 >::_set( const Quadric< 1 > &Q )
782 {
783 _quadric = Q;
784 _Q = Q._Q;
785 try{ _Qinv = _Q.inverse(); } catch(...){}
786 }
787
788 inline bool QuadricBoundingBoxOverlap< 1 >::_intersect( const Quadric< 1 > &Q , const BoundingBox< 1 > &bBox ) const
789 {
790 if( Q(bBox[0])<0 || Q(bBox[1])<0 ) return true;
791 Point< 1 > extremum = -_Qinv * Q._L;
792 return extremum[0]>bBox[0][0] && extremum[0]<bBox[1][0] && Q(extremum)<0;
793 }
794
796 // RotationParameter //
798 template< typename RotationParameterType , typename ParameterType >
799 RotationParameterType RotationParameter< RotationParameterType , ParameterType >::operator * ( double scale ) const
800 {
801 RotationParameterType p;
802 p.parameter = parameter * scale;
803 return p;
804 }
805
806 template< typename RotationParameterType , typename ParameterType >
807 RotationParameterType RotationParameter< RotationParameterType , ParameterType >::operator + ( const RotationParameterType &p ) const
808 {
809 RotationParameterType _p;
810 _p.parameter = parameter + p.parameter;
811 return _p;
812 }
813
815 // TransformationParameter //
817 template< typename RotationParameterType >
818 TransformationParameter< RotationParameterType > TransformationParameter< RotationParameterType >::operator * ( double scale ) const
819 {
820 TransformationParameter tp = *this;
821 tp.rotationParameter *= scale;
822 tp.translation *= scale;
823 return tp;
824 }
825
826 template< typename RotationParameterType >
827 TransformationParameter< RotationParameterType > TransformationParameter< RotationParameterType >::operator + ( const TransformationParameter &tp ) const
828 {
829 TransformationParameter _tp = *this;
830 _tp.rotationParameter += tp.rotationParameter;
831 _tp.translation += tp.translation;
832 return _tp;
833 }
834
835 template< typename RotationParameterType >
836 TransformationParameter< RotationParameterType >::TransformationParameter( void ) {}
837
838 template< typename RotationParameterType >
839 TransformationParameter< RotationParameterType >::TransformationParameter( const Matrix4D &m ) : rotationParameter( Matrix3D( m ) )
840 {
841 for( int i=0 ; i<3 ; i++ ) translation[i] = m(i,3);
842 }
843
844 template< typename RotationParameterType >
845 TransformationParameter< RotationParameterType >::TransformationParameter( const Matrix4D &m , const TransformationParameter &p ) : rotationParameter( Matrix3D( m ) , p.rotationParameter )
846 {
847 for( int i=0 ; i<3 ; i++ ) translation[i] = m(i,3);
848 }
849
850 template< typename RotationParameterType >
851 Matrix4D TransformationParameter< RotationParameterType >::operator () ( void ) const
852 {
853 return Matrix4D( rotationParameter() , translation );
854 }
855}
MatrixMNTC.h
Identity
const GXMatrixMNTC< Coord > Identity(unsigned int cOrder)
GXMatrixMNd
GXMatrixMNTC< double > GXMatrixMNd
Definition MatrixMNTC.h:150
SVDFit.h
SVDMat
void SVDMat(GXMatrixMNTC< REAL > &A, GXMatrixMNTC< REAL > &u, GXMatrixMNTC< REAL > &w, GXMatrixMNTC< REAL > &vt)
Definition SVDFit.inl:59
Util::_BaseMatrix
Definition geometry.h:126
Util::_BaseMatrix::operator*
Point< Rows > operator*(const Point< Cols > &p) const
Definition geometry.inl:169
Util::_BaseMatrix::operator()
double & operator()(int r, int c)
Definition geometry.inl:155
Util::_BaseMatrix::_m
double _m[Rows][Cols]
Definition geometry.h:128
Util::_BaseMatrix::_BaseMatrix
_BaseMatrix(void)
Definition geometry.inl:152
Util::_BaseMatrix::operator+
MatrixType operator+(const MatrixType &m) const
Definition geometry.inl:180
Util::_BaseMatrix::dot
double dot(const _BaseMatrix &p) const
Definition geometry.inl:144
Util::_BaseMatrix::transpose
MatrixTransposeType transpose(void) const
Definition geometry.inl:161
Util::_InnerProductSpace::squareNorm
double squareNorm(void) const
Definition algebra.h:117
Util::_InnerProductSpace::unit
Element unit(void) const
Definition algebra.h:123
Util::_InnerProductSpace< Point< Dim > >::Dot
static double Dot(const Point< Dim > &e1, const Point< Dim > &e2)
Definition algebra.h:102
Util::BoundingBox
Definition geometry.h:345
Util::BoundingBox::BoundingBox
BoundingBox(void)
Definition geometry.inl:503
Util::BoundingBox::isEmpty
bool isEmpty(void) const
Definition geometry.inl:580
Util::BoundingBox::operator^=
BoundingBox & operator^=(const BoundingBox &b)
Definition geometry.inl:557
Util::BoundingBox::operator+=
BoundingBox & operator+=(const BoundingBox &b)
Definition geometry.inl:543
Util::BoundingBox::isInside
bool isInside(const Point< Dim > &p) const
Definition geometry.inl:573
Util::BoundingBox::operator+
BoundingBox operator+(const BoundingBox &b) const
Definition geometry.inl:528
Util::BoundingBox::operator[]
Point< Dim > & operator[](int index)
Definition geometry.inl:522
Util::BoundingBox::operator^
BoundingBox operator^(const BoundingBox &b) const
Definition geometry.inl:546
Util::BoundingBox::_p
Point< Dim > _p[2]
Definition geometry.h:350
Util::Matrix
Definition geometry.h:168
Util::Matrix::Matrix
Matrix(void)
Definition geometry.inl:185
Util::Matrix::operator*
Matrix< Rows, _Cols > operator*(const Matrix< Cols, _Cols > &m) const
Definition geometry.inl:189
Util::Plane::distance
double distance
Definition geometry.h:274
Util::Plane::normal
Point< Dim > normal
Definition geometry.h:272
Util::Plane::Plane
Plane(void)
Definition geometry.inl:432
Util::Plane< 2 >::operator()
double operator()(const Point< Dim > &p) const
Util::Point
Definition geometry.h:52
Util::Point::CrossProduct
static Point CrossProduct(Points ... points)
Util::Point::operator*
Point operator*(const Point &p) const
Definition geometry.inl:83
Util::Point::operator[]
double & operator[](int index)
Definition geometry.inl:77
Util::Point::operator/=
Point & operator/=(const Point &p)
Definition geometry.inl:102
Util::Point::Point
Point(void)
Definition geometry.inl:48
Util::Point::operator/
Point operator/(const Point &p) const
Definition geometry.inl:91
Util::Point::_init
void _init(const double *values, unsigned int sz)
Definition geometry.inl:40
Util::Point::_p
double _p[Dim]
Definition geometry.h:54
Util::Point::operator+
Point operator+(const Point &p) const
Definition geometry.inl:66
Util::Point::operator*=
Point & operator*=(const Point &p)
Definition geometry.inl:99
Util::Point::dot
double dot(const Point &p) const
Definition geometry.inl:69
Util::Quadric
Definition geometry.h:399
Util::Quadric::getConstant
double getConstant(void) const
Definition geometry.inl:633
Util::Quadric::operator*
Quadric< _Dim > operator*(const Matrix< Dim, _Dim > &T) const
Definition geometry.inl:649
Util::Quadric::Quadric
friend class Quadric
Definition geometry.h:405
Util::Quadric::getLinear
Point< Dim > getLinear(void) const
Definition geometry.inl:630
Util::Quadric::operator()
Matrix< Dim+1, Dim+1 > operator()(void) const
Definition geometry.inl:614
Util::Quadric::setQuadratic
void setQuadratic(Matrix< Dim, Dim > Q)
Definition geometry.inl:636
Util::Quadric::setExtremum
bool setExtremum(Point< Dim > &extremum) const
Definition geometry.inl:675
Util::Quadric::setLinear
void setLinear(Point< Dim > L)
Definition geometry.inl:639
Util::Quadric::_Q
Matrix< Dim, Dim > _Q
Definition geometry.h:400
Util::Quadric::getQuadratic
Matrix< Dim, Dim > getQuadratic(void) const
Definition geometry.inl:627
Util::Quadric::setConstant
void setConstant(double C)
Definition geometry.inl:642
Util::Quadric::_L
Point< Dim > _L
Definition geometry.h:401
Util::Quadric::operator+
Quadric< Dim > operator+(const Point< Dim > &t) const
Definition geometry.inl:662
Util::Quadric::_C
double _C
Definition geometry.h:402
Util::Ray
Definition geometry.h:299
Util::Ray::operator+
Ray operator+(const Point< Dim > &p) const
Definition geometry.inl:481
Util::Ray::direction
Point< Dim > direction
Definition geometry.h:305
Util::Ray::Ray
Ray(void)
Definition geometry.inl:472
Util::Ray::position
Point< Dim > position
Definition geometry.h:302
Util::Ray< 1 >::operator()
Point< Dim > operator()(double t) const
Util::Ray::operator-
Ray operator-(const Point< Dim > &p) const
Definition geometry.inl:487
Util::Ray::operator+=
Ray & operator+=(const Point< Dim > &p)
Definition geometry.inl:484
Util::Ray::operator-=
Ray & operator-=(const Point< Dim > &p)
Definition geometry.inl:490
Util::RotationParameter::operator+
RotationParameterType operator+(const RotationParameterType &p) const
Definition geometry.inl:807
Util::RotationParameter::parameter
ParameterType parameter
Definition geometry.h:642
Util::RotationParameter::operator*
RotationParameterType operator*(double scale) const
Definition geometry.inl:799
Util::TransformationParameter::operator+
TransformationParameter operator+(const TransformationParameter &tp) const
Definition geometry.inl:827
Util::TransformationParameter::rotationParameter
RotationParameterType rotationParameter
Definition geometry.h:751
Util::TransformationParameter::TransformationParameter
TransformationParameter(void)
Definition geometry.inl:836
Util::TransformationParameter::translation
Point3D translation
Definition geometry.h:754
Util::TransformationParameter::operator*
TransformationParameter operator*(double scale) const
Definition geometry.inl:818
Util::TransformationParameter::operator()
Matrix4D operator()(void) const
Definition geometry.inl:851
exceptions.h
ERROR_OUT
#define ERROR_OUT(...)
Definition exceptions.h:154
THROW
#define THROW(...)
Definition exceptions.h:151
b
long b
Definition jpegint.h:371
Ray
Definition box.h:7
Util
Definition algebra.h:7
Util::Matrix3D
Matrix< 3, 3 > Matrix3D
Definition geometry.h:541
Util::Matrix4D
Matrix< 4, 4 > Matrix4D
Definition geometry.h:544
Util::operator<<
std::ostream & operator<<(std::ostream &stream, const Point< Dim > &p)
Definition geometry.inl:127
Util::operator>>
std::istream & operator>>(std::istream &stream, Point< Dim > &p)
Definition geometry.inl:134
Util::operator*
Ray< Dim > operator*(const Matrix< Dim+1, Dim+1 > &m, const Ray< Dim > &ray)
Definition geometry.inl:493
Util::SquareMatrix
Matrix< Dim, Dim > SquareMatrix
Definition geometry.h:532
Util::QuadricBoundingBoxOverlap< 1 >::_quadric
Quadric< 1 > _quadric
Definition geometry.h:496
Util::QuadricBoundingBoxOverlap< 1 >::_Q
Matrix< 1, 1 > _Q
Definition geometry.h:497
Util::QuadricBoundingBoxOverlap< 1 >::_Qinv
Matrix< 1, 1 > _Qinv
Definition geometry.h:497
Util::QuadricBoundingBoxOverlap::_Qinv
Matrix< Dim, Dim > _Qinv
Definition geometry.h:476
Util::QuadricBoundingBoxOverlap::_T
Matrix< Dim+1, Dim > _T
Definition geometry.h:477
Util::QuadricBoundingBoxOverlap::_Tt_Q
Matrix< Dim, Dim+1 > _Tt_Q
Definition geometry.h:478
Util::QuadricBoundingBoxOverlap::_Q
Matrix< Dim, Dim > _Q
Definition geometry.h:476
Util::QuadricBoundingBoxOverlap::_intersect
bool _intersect(const Quadric< Dim > &Q, const BoundingBox< Dim > &bBox) const
Definition geometry.inl:718
Util::QuadricBoundingBoxOverlap::QuadricBoundingBoxOverlap
friend struct QuadricBoundingBoxOverlap
Definition geometry.h:474
Util::QuadricBoundingBoxOverlap::_boundaryInfo
struct QuadricBoundingBoxOverlap< Dim-1 > _boundaryInfo[Dim]
Definition geometry.h:479
Util::QuadricBoundingBoxOverlap::_quadric
Quadric< Dim > _quadric
Definition geometry.h:475
Util::QuadricBoundingBoxOverlap::operator()
bool operator()(const BoundingBox< Dim > &bBox) const
Definition geometry.inl:698
Util::QuadricBoundingBoxOverlap::_set
void _set(const Quadric< Dim > &Q)
Definition geometry.inl:701