mirror of
https://github.com/CloudCompare/PoissonRecon.git
synced 2026-08-29 16:40:28 +08:00
mkazhdan
630 lines
27 KiB
C++
630 lines
27 KiB
C++
/*
|
|
Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
|
|
All rights reserved.
|
|
|
|
Redistribution and use in source and binary forms, with or without modification,
|
|
are permitted provided that the following conditions are met:
|
|
|
|
Redistributions of source code must retain the above copyright notice, this list of
|
|
conditions and the following disclaimer. Redistributions in binary form must reproduce
|
|
the above copyright notice, this list of conditions and the following disclaimer
|
|
in the documentation and/or other materials provided with the distribution.
|
|
|
|
Neither the name of the Johns Hopkins University nor the names of its contributors
|
|
may be used to endorse or promote products derived from this software without specific
|
|
prior written permission.
|
|
|
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
|
|
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
|
|
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
|
|
SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
|
|
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
|
|
TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
|
|
BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
|
|
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
|
|
ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
|
|
DAMAGE.
|
|
*/
|
|
|
|
#ifndef GEOMETRY_INCLUDED
|
|
#define GEOMETRY_INCLUDED
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <vector>
|
|
#include <stdlib.h>
|
|
#include <unordered_map>
|
|
#include <string.h>
|
|
#ifdef _WIN32
|
|
#include <io.h>
|
|
#endif // _WIN32
|
|
|
|
template< class Real > Real Random( void );
|
|
|
|
template< class Real , unsigned int Dim > struct XForm;
|
|
|
|
template< class Real , unsigned int Dim >
|
|
struct Point
|
|
{
|
|
void _init( unsigned int d )
|
|
{
|
|
if( !d ) memset( coords , 0 , sizeof(Real)*Dim );
|
|
else ERROR_OUT( "Should never be called" );
|
|
}
|
|
template< class _Real , class ... _Reals > void _init( unsigned int d , _Real v , _Reals ... values )
|
|
{
|
|
coords[d] = (Real)v;
|
|
if( d+1<Dim ) _init( d+1 , values... );
|
|
}
|
|
template< class ... Points >
|
|
static void _AddColumnVector( XForm< Real , Dim >& x , unsigned int c , Point point , Points ... points )
|
|
{
|
|
for( unsigned int r=0 ; r<Dim ; r++ ) x( c , r ) = point[r];
|
|
_AddColumnVector( x , c+1 , points ... );
|
|
}
|
|
static void _AddColumnVector( XForm< Real , Dim >& x , unsigned int c ){ ; }
|
|
public:
|
|
Real coords[Dim];
|
|
Point( void ) { memset( coords , 0 , sizeof(Real)*Dim ); }
|
|
Point( const Point& p ){ memcpy( coords , p.coords , sizeof(Real)*Dim ); }
|
|
template< class ... _Reals > Point( _Reals ... values ){ static_assert( sizeof...(values)==Dim || sizeof...(values)==0 , "[ERROR] Point::Point: Invalid number of coefficients" ) ; _init( 0 , values... ); }
|
|
template< class _Real > Point( const Point< _Real , Dim >& p ){ for( unsigned int d=0 ; d<Dim ; d++ ) coords[d] = (Real) p.coords[d]; }
|
|
inline Real& operator[] ( unsigned int i ) { return coords[i]; }
|
|
inline const Real& operator[] ( unsigned int i ) const { return coords[i]; }
|
|
inline Point operator - ( void ) const { Point q ; for( unsigned int d=0 ; d<Dim ; d++ ) q.coords[d] = - coords[d] ; return q; }
|
|
|
|
|
|
template< class _Real > inline Point& operator += ( Point< _Real , Dim > p ) { for( unsigned int d=0 ; d<Dim ; d++ ) coords[d] += (Real)p.coords[d] ; return *this; }
|
|
template< class _Real > inline Point operator + ( Point< _Real , Dim > p ) const { Point q ; for( unsigned int d=0 ; d<Dim ; d++ ) q.coords[d] = coords[d] + (Real)p.coords[d] ; return q; }
|
|
template< class _Real > inline Point& operator -= ( Point< _Real , Dim > p ) { return (*this)+=(-p); }
|
|
template< class _Real > inline Point operator - ( Point< _Real , Dim > p ) const { return (*this)+ (-p); }
|
|
template< class Scalar > inline Point& operator *= ( Scalar r ) { for( unsigned int d=0 ; d<Dim ; d++ ) coords[d] *= r ; return *this; }
|
|
template< class Scalar > inline Point operator * ( Scalar r ) const { Point q ; for( unsigned int d=0 ; d<Dim ; d++ ) q.coords[d] = coords[d] * r ; return q; }
|
|
template< class Scalar > inline Point& operator /= ( Scalar r ) { for( unsigned int d=0 ; d<Dim ; d++ ) coords[d] /= r ; return *this; }
|
|
template< class Scalar > inline Point operator / ( Scalar r ) const { Point q ; for( unsigned int d=0 ; d<Dim ; d++ ) q.coords[d] = coords[d] / r ; return q; }
|
|
template< class _Real > inline Point& operator *= ( Point< _Real , Dim > p ) { for( unsigned int d=0 ; d<Dim ; d++ ) coords[d] *= p.coords[d] ; return *this; }
|
|
template< class _Real > inline Point operator * ( Point< _Real , Dim > p ) const { Point q ; for( unsigned int d=0 ; d<Dim ; d++ ) q.coords[d] = coords[d] * p.coords[d] ; return q; }
|
|
template< class _Real > inline Point& operator /= ( Point< _Real , Dim > p ) { for( unsigned int d=0 ; d<Dim ; d++ ) coords[d] /= p.coords[d] ; return *this; }
|
|
template< class _Real > inline Point operator / ( Point< _Real , Dim > p ) const { Point q ; for( unsigned int d=0 ; d<Dim ; d++ ) q.coords[d] = coords[d] / p.coords[d] ; return q; }
|
|
|
|
static Real Dot( const Point& p1 , const Point& p2 ){ Real dot = {} ; for( unsigned int d=0 ; d<Dim ; d++ ) dot += p1.coords[d] * p2.coords[d] ; return dot; }
|
|
static Real SquareNorm( const Point& p ){ return Dot( p , p ); }
|
|
template< class ... Points > static Point CrossProduct( Points ... points )
|
|
{
|
|
static_assert( sizeof ... ( points )==Dim-1 , "Number of points in cross-product must be one less than the dimension" );
|
|
XForm< Real , Dim > x;
|
|
_AddColumnVector( x , 0 , points ... );
|
|
Point p;
|
|
for( unsigned int d=0 ; d<Dim ; d++ ) p[d] = ( d&1 ) ? -x.subDeterminant( Dim-1 , d ) : x.subDeterminant( Dim-1 , d );
|
|
return p;
|
|
}
|
|
static Point CrossProduct( const Point* points )
|
|
{
|
|
XForm< Real , Dim > x;
|
|
for( unsigned int d=0 ; d<Dim-1 ; d++ ) for( unsigned int c=0 ; c<Dim ; c++ ) x(d,c) = points[d][c];
|
|
Point p;
|
|
for( unsigned int d=0 ; d<Dim ; d++ ) p[d] = ( d&1 ) ? -x.subDeterminant( Dim-1 , d ) : x.subDeterminant( Dim-1 , d );
|
|
return p;
|
|
}
|
|
static Point CrossProduct( Point* points ){ return CrossProduct( ( const Point* )points ); }
|
|
};
|
|
template< class Real , unsigned int Dim > Point< Real , Dim > operator * ( Real r , Point< Real , Dim > p ){ return p*r; }
|
|
template< class Real , unsigned int Dim > Point< Real , Dim > operator / ( Real r , Point< Real , Dim > p ){ return p/r; }
|
|
|
|
template< class Real , unsigned int _Columns , unsigned int _Rows >
|
|
struct Matrix
|
|
{
|
|
static const unsigned int Columns = _Columns;
|
|
static const unsigned int Rows = _Rows;
|
|
Real coords[Columns][Rows];
|
|
Matrix( void ) { memset( coords , 0 , sizeof(coords) ); }
|
|
inline Real& operator() ( unsigned int c , unsigned int r ) { return coords[c][r]; }
|
|
inline const Real& operator() ( unsigned int c , unsigned int r ) const { return coords[c][r]; }
|
|
inline Real* operator[] ( unsigned int c ) { return coords[c] ; }
|
|
inline const Real* operator[] ( unsigned int c ) const { return coords[c] ; }
|
|
|
|
inline Matrix operator - ( void ) const { Matrix m ; for( unsigned int c=0 ; c<Columns ; c++ ) for( unsigned int r=0 ; r<Rows ; r++ ) m.coords[c][r] = - coords[c][r] ; return m; }
|
|
|
|
inline Matrix& operator += ( const Matrix& m ){ for( unsigned int c=0 ; c<Columns ; c++ ) for( unsigned int r=0 ; r<Rows ; r++ ) coords[c][r] += m.coords[c][r] ; return *this; }
|
|
inline Matrix operator + ( const Matrix& m ) const { Matrix n ; for( unsigned int c=0 ; c<Columns ; c++ ) for( unsigned int r=0 ; r<Rows ; r++ ) n.coords[c][r] = coords[c][r] + m.coords[c][r] ; return n; }
|
|
inline Matrix& operator *= ( Real s ) { for( unsigned int c=0 ; c<Columns ; c++ ) for( unsigned int r=0 ; r<Rows ; r++ ) coords[c][r] *= s ; return *this; }
|
|
inline Matrix operator * ( Real s ) const { Matrix n ; for( unsigned int c=0 ; c<Columns ; c++ ) for( unsigned int r=0 ; r<Rows ; r++ ) n.coords[c][r] = coords[c][r] * s ; return n; }
|
|
|
|
inline Matrix& operator -= ( const Matrix& m ){ return ( (*this)+=(-m) ); }
|
|
inline Matrix operator - ( const Matrix& m ) const { return (*this)+(-m); }
|
|
inline Matrix& operator /= ( Real s ){ return ( (*this)*=(Real)(1./s) ); }
|
|
inline Matrix operator / ( Real s ) const { return (*this) * ( (Real)(1./s) ); }
|
|
|
|
static Real Dot( const Matrix& m1 , const Matrix& m2 ){ Real dot = (Real)0 ; for( unsigned int c=0 ; c<Columns ; c++ ) for( unsigned int r=0 ; r<Rows ; r++ ) dot += m1.coords[c][r] * m2.coords[c][r] ; return dot; }
|
|
|
|
template< typename T >
|
|
inline Point< T , Rows > operator* ( const Point< T , Columns >& p ) const { Point< T , Rows > q ; for( unsigned int c=0 ; c<Columns ; c++ ) for( unsigned int r=0 ; r<Rows ; r++ ) q[r] += (T)( p[c] * coords[c][r] ) ; return q; }
|
|
};
|
|
|
|
template< class Real , unsigned int Dim >
|
|
struct XForm
|
|
{
|
|
Real coords[Dim][Dim];
|
|
XForm( void ) { memset( coords , 0 , sizeof(Real) * Dim * Dim ); }
|
|
static XForm Identity( void )
|
|
{
|
|
XForm xForm;
|
|
for( unsigned int d=0 ; d<Dim ; d++ ) xForm(d,d) = (Real)1.;
|
|
return xForm;
|
|
}
|
|
Real& operator() ( unsigned int i , unsigned int j ){ return coords[i][j]; }
|
|
const Real& operator() ( unsigned int i , unsigned int j ) const { return coords[i][j]; }
|
|
template< class _Real > Point< _Real , Dim-1 > operator * ( const Point< _Real , Dim-1 >& p ) const
|
|
{
|
|
Point< _Real , Dim-1 > q;
|
|
for( unsigned int i=0 ; i<Dim-1 ; i++ )
|
|
{
|
|
for( unsigned int j=0 ; j<Dim-1 ; j++ ) q[i] += (_Real)( coords[j][i] * p[j] );
|
|
q[i] += (_Real)coords[Dim-1][i];
|
|
}
|
|
return q;
|
|
}
|
|
template< class _Real > Point< _Real , Dim > operator * ( const Point< _Real , Dim >& p ) const
|
|
{
|
|
Point< _Real , Dim > q;
|
|
for( unsigned int i=0 ; i<Dim ; i++ ) for( unsigned int j=0 ; j<Dim ; j++ ) q[i] += (_Real)( coords[j][i] * p[j] );
|
|
return q;
|
|
}
|
|
XForm operator * ( const XForm& m ) const
|
|
{
|
|
XForm n;
|
|
for( unsigned int i=0 ; i<Dim ; i++ ) for( unsigned int j=0 ; j<Dim ; j++ ) for( unsigned int k=0 ; k<Dim ; k++ ) n.coords[i][j] += m.coords[i][k]*coords[k][j];
|
|
return n;
|
|
}
|
|
XForm transpose( void ) const
|
|
{
|
|
XForm xForm;
|
|
for( unsigned int i=0 ; i<Dim ; i++ ) for( unsigned int j=0 ; j<Dim ; j++ ) xForm( i , j ) = coords[j][i];
|
|
return xForm;
|
|
}
|
|
Real determinant( void ) const
|
|
{
|
|
Real det = (Real)0.;
|
|
for( unsigned int d=0 ; d<Dim ; d++ )
|
|
if( d&1 ) det -= coords[d][0] * subDeterminant( d , 0 );
|
|
else det += coords[d][0] * subDeterminant( d , 0 );
|
|
return det;
|
|
}
|
|
XForm inverse( void ) const
|
|
{
|
|
XForm xForm;
|
|
Real d = determinant();
|
|
for( unsigned int i=0 ; i<Dim ; i++ ) for( unsigned int j=0 ; j<Dim ; j++ )
|
|
if( (i+j)%2==0 ) xForm.coords[j][i] = subDeterminant( i , j ) / d;
|
|
else xForm.coords[j][i] = -subDeterminant( i , j ) / d;
|
|
return xForm;
|
|
}
|
|
Real subDeterminant( unsigned int i , unsigned int j ) const
|
|
{
|
|
XForm< Real , Dim-1 > xForm;
|
|
unsigned int ii[Dim-1] , jj[Dim-1];
|
|
for( unsigned int a=0 , _i=0 , _j=0 ; a<Dim ; a++ )
|
|
{
|
|
if( a!=i ) ii[_i++] = a;
|
|
if( a!=j ) jj[_j++] = a;
|
|
}
|
|
for( unsigned int _i=0 ; _i<Dim-1 ; _i++ ) for( unsigned int _j=0 ; _j<Dim-1 ; _j++ ) xForm( _i , _j ) = coords[ ii[_i] ][ jj[_j] ];
|
|
return xForm.determinant();
|
|
}
|
|
|
|
inline XForm operator - ( void ) const { XForm m ; for( unsigned int c=0 ; c<Dim ; c++ ) for( unsigned int r=0 ; r<Dim ; r++ ) m.coords[c][r] = - coords[c][r] ; return m; }
|
|
|
|
inline XForm& operator += ( const XForm& m ){ for( unsigned int c=0 ; c<Dim ; c++ ) for( unsigned int r=0 ; r<Dim ; r++ ) coords[c][r] += m.coords[c][r] ; return *this; }
|
|
inline XForm operator + ( const XForm& m ) const { XForm n ; for( unsigned int c=0 ; c<Dim ; c++ ) for( unsigned int r=0 ; r<Dim ; r++ ) n.coords[c][r] = coords[c][r] + m.coords[c][r] ; return n; }
|
|
inline XForm& operator *= ( Real s ) { for( unsigned int c=0 ; c<Dim ; c++ ) for( unsigned int r=0 ; r<Dim ; r++ ) coords[c][r] *= s ; return *this; }
|
|
inline XForm operator * ( Real s ) const { XForm n ; for( unsigned int c=0 ; c<Dim ; c++ ) for( unsigned int r=0 ; r<Dim ; r++ ) n.coords[c][r] = coords[c][r] * s ; return n; }
|
|
|
|
inline XForm& operator -= ( const XForm& m ){ return ( (*this)+=(-m) ); }
|
|
inline XForm operator - ( const XForm& m ) const { return (*this)+(-m); }
|
|
inline XForm& operator /= ( Real s ){ return ( (*this)*=(Real)(1./s) ); }
|
|
inline XForm operator / ( Real s ) const { return (*this) * ( (Real)(1./s) ); }
|
|
};
|
|
template<>
|
|
inline XForm< float , 1 > XForm< float , 1 >::inverse( void ) const
|
|
{
|
|
XForm< float , 1 > x;
|
|
x.coords[0][0] = (float)(1./coords[0][0] );
|
|
return x;
|
|
}
|
|
template<>
|
|
inline XForm< double , 1 > XForm< double , 1 >::inverse( void ) const
|
|
{
|
|
XForm< double , 1 > x;
|
|
x.coords[0][0] = (double)(1./coords[0][0] );
|
|
return x;
|
|
}
|
|
template<> inline float XForm< float , 1 >::determinant( void ) const { return coords[0][0]; }
|
|
template<> inline double XForm< double , 1 >::determinant( void ) const { return coords[0][0]; }
|
|
|
|
template< class Real , unsigned int Dim >
|
|
struct OrientedPoint
|
|
{
|
|
Point< Real , Dim > p , n;
|
|
OrientedPoint( Point< Real , Dim > pp = Point< Real , Dim >() , Point< Real , Dim > nn=Point< Real , Dim >() ) : p(pp) , n(nn) { ; }
|
|
template< class _Real > OrientedPoint( const OrientedPoint< _Real , Dim>& p ) : OrientedPoint( Point< Real , Dim >( p.p ) , Point< Real , Dim >( p.n ) ){ ; }
|
|
|
|
template< class _Real > inline OrientedPoint& operator += ( OrientedPoint< _Real , Dim > _p ){ p += _p.p , n += _p.n ; return *this; }
|
|
template< class _Real > inline OrientedPoint operator + ( OrientedPoint< _Real , Dim > _p ) const { return OrientedPoint< Real , Dim >( p+_p.p , n+_p.n ); }
|
|
template< class _Real > inline OrientedPoint& operator *= ( _Real r ) { p *= r , n *= r ; return *this; }
|
|
template< class _Real > inline OrientedPoint operator * ( _Real r ) const { return OrientedPoint< Real , Dim >( p*r , n*r ); }
|
|
|
|
template< class _Real > inline OrientedPoint& operator -= ( OrientedPoint< _Real , Dim > p ){ return ( (*this)+=(-p) ); }
|
|
template< class _Real > inline OrientedPoint operator - ( OrientedPoint< _Real , Dim > p ) const { return (*this)+(-p); }
|
|
template< class _Real > inline OrientedPoint& operator /= ( _Real r ){ return ( (*this)*=Real(1./r) ); }
|
|
template< class _Real > inline OrientedPoint operator / ( _Real r ) const { return (*this) * ( Real(1.)/r ); }
|
|
};
|
|
|
|
|
|
template< class Data , class Real >
|
|
struct ProjectiveData
|
|
{
|
|
Data data;
|
|
Real weight;
|
|
ProjectiveData( Data d=Data() , Real w=(Real)0 ) : data(d) , weight(w) { ; }
|
|
operator Data (){ return weight!=0 ? data/weight : data*weight; }
|
|
Data value( void ) const { return weight!=0 ? data/weight : data*weight; }
|
|
ProjectiveData& operator += ( const ProjectiveData& p ){ data += p.data , weight += p.weight ; return *this; }
|
|
ProjectiveData& operator -= ( const ProjectiveData& p ){ data -= p.data , weight -= p.weight ; return *this; }
|
|
ProjectiveData& operator *= ( Real s ){ data *= s , weight *= s ; return *this; }
|
|
ProjectiveData& operator /= ( Real s ){ data /= s , weight /= s ; return *this; }
|
|
ProjectiveData operator + ( const ProjectiveData& p ) const { return ProjectiveData( data+p.data , weight+p.weight ); }
|
|
ProjectiveData operator - ( const ProjectiveData& p ) const { return ProjectiveData( data-p.data , weight-p.weight ); }
|
|
ProjectiveData operator * ( Real s ) const { return ProjectiveData( data*s , weight*s ); }
|
|
ProjectiveData operator / ( Real s ) const { return ProjectiveData( data/s , weight/s ); }
|
|
};
|
|
|
|
template< class Real , unsigned int Dim > Point< Real , Dim > RandomBallPoint( void );
|
|
template< class Real , unsigned int Dim > Point< Real , Dim > RandomSpherePoint( void );
|
|
template< class Real , unsigned int Dim > Real Length( Point< Real , Dim > p ){ return (Real)sqrt( Point< Real , Dim >::SquareNorm( p ) ); }
|
|
template< class Real , unsigned int Dim > Real SquareLength( Point< Real , Dim > p ){ return Point< Real , Dim >::SquareNorm( p ); }
|
|
template< class Real , unsigned int Dim > Real Distance( Point< Real , Dim > p1 , Point< Real , Dim > p2 ){ return Length(p1-p2); }
|
|
template< class Real , unsigned int Dim > Real SquareDistance( Point< Real , Dim > p1 , Point< Real , Dim > p2 ){ return SquareLength( p1-p2 ); }
|
|
template< class Real > Point< Real , 3 > CrossProduct( Point< Real , 3 > p1 , Point< Real , 3 > p2 ){ return Point< Real , 3 >::CrossProduct( p1 , p2 ); }
|
|
|
|
template< class Real , unsigned int Dim > Real SquareArea( Point< Real , Dim > p1 , Point< Real , Dim > p2 , Point< Real , Dim > p3 )
|
|
{
|
|
Point< Real , Dim > v1 = p2-p1 , v2 = p3-p1;
|
|
// Area^2 = ( |v1|^2 * |v2|^2 * sin^2( < v1 ,v2 ) ) / 4
|
|
// = ( |v1|^2 * |v2|^2 * ( 1 - cos^2( < v1 ,v2 ) ) ) / 4
|
|
// = ( |v1|^2 * |v2|^2 * ( 1 - < v1 , v2 >^2 / ( |v1|^2 * |v2|^2 ) ) ) / 4
|
|
// = ( |v1|^2 * |v2|^2 - < v1 , v2 >^2 ) / 4
|
|
Real dot = Point< Real , Dim >::Dot( v1 , v2 );
|
|
Real l1 = Point< Real , Dim >::SquareNorm( v1 ) , l2 = Point< Real , Dim >::SquareNorm( v2 );
|
|
return ( l1 * l2 - dot * dot ) / 4;
|
|
}
|
|
template< class Real , unsigned int Dim > Real Area( Point< Real , Dim > p1 , Point< Real , Dim > p2 , Point< Real , Dim > p3 ){ return (Real)sqrt( SquareArea( p1 , p2 , p3 ) ); }
|
|
|
|
template< unsigned int K > struct Factorial{ static const unsigned long long Value = Factorial< K-1 >::Value * K; };
|
|
template<> struct Factorial< 0 >{ static const unsigned long long Value = 1; };
|
|
|
|
template< class Real , unsigned int Dim , unsigned int K >
|
|
struct Simplex
|
|
{
|
|
Point< Real , Dim > p[K+1];
|
|
Simplex( void ){ static_assert( K<=Dim , "[ERROR] Bad simplex dimension" ); }
|
|
Point< Real , Dim >& operator[]( unsigned int k ){ return p[k]; }
|
|
const Point< Real , Dim >& operator[]( unsigned int k ) const { return p[k]; }
|
|
Real measure( void ) const { return (Real)sqrt( squareMeasure() ); }
|
|
Real squareMeasure( void ) const
|
|
{
|
|
XForm< Real , K > mass;
|
|
for( unsigned int i=1 ; i<=K ; i++ ) for( unsigned int j=1 ; j<=K ; j++ ) mass(i-1,j-1) = Point< Real , Dim >::Dot( p[i]-p[0] , p[j]-p[0] );
|
|
return mass.determinant() / ( Factorial< K >::Value * Factorial< K >::Value );
|
|
}
|
|
Point< Real , Dim > center( void ) const
|
|
{
|
|
Point< Real , Dim > c;
|
|
for( unsigned int k=0 ; k<=K ; k++ ) c += p[k];
|
|
return c / (K+1);
|
|
}
|
|
void split( Point< Real , Dim > pNormal , Real pOffset , std::vector< Simplex >& back , std::vector< Simplex >& front ) const;
|
|
};
|
|
template< class Real , unsigned int Dim >
|
|
struct Simplex< Real , Dim , 0 >
|
|
{
|
|
Point< Real , Dim > p[1];
|
|
Point< Real , Dim >& operator[]( unsigned int k ){ return p[k]; }
|
|
const Point< Real , Dim >& operator[]( unsigned int k ) const { return p[k]; }
|
|
Real squareMeasure( void ) const { return (Real)1.; }
|
|
Real measure( void ) const { return (Real)1.; }
|
|
Point< Real , Dim > center( void ) const { return p[0]; }
|
|
void split( Point< Real , Dim > pNormal , Real pOffset , std::vector< Simplex >& back , std::vector< Simplex >& front ) const
|
|
{
|
|
if( Point< Real , Dim >::Dot( p[0] , pNormal ) < pOffset ) back.push_back( *this );
|
|
else front.push_back( *this );
|
|
}
|
|
};
|
|
template< class Real , unsigned int Dim > using Edge = Simplex< Real , Dim , 1 >;
|
|
template< class Real , unsigned int Dim > using Triangle = Simplex< Real , Dim , 2 >;
|
|
|
|
template< unsigned int K , typename Index >
|
|
struct SimplexIndex
|
|
{
|
|
Index idx[K+1];
|
|
template< class ... Ints >
|
|
SimplexIndex( Ints ... values ){ static_assert( sizeof...(values)==K+1 || sizeof...(values)==0 , "[ERROR] Invalid number of coefficients" ) ; _init( 0 , values ... ); }
|
|
Index &operator[] ( unsigned int i ) { return idx[i] ;}
|
|
const Index &operator[] ( unsigned int i ) const { return idx[i]; }
|
|
protected:
|
|
void _init( unsigned int k )
|
|
{
|
|
if( !k ) memset( idx , 0 , sizeof(idx) );
|
|
else ERROR_OUT( "Should never be called" );
|
|
}
|
|
template< class ... Ints > void _init( unsigned int k , Index v , Ints ... values )
|
|
{
|
|
idx[k] = v;
|
|
if( k<=K ) _init( k+1 , values ... );
|
|
}
|
|
};
|
|
template< typename Index > using EdgeIndex = SimplexIndex< 1 , Index >;
|
|
template< typename Index > using TriangleIndex = SimplexIndex< 2 , Index >;
|
|
|
|
template< typename Index >
|
|
class CoredPointIndex
|
|
{
|
|
public:
|
|
Index index;
|
|
char inCore;
|
|
|
|
bool operator == (const CoredPointIndex& cpi) const {return (index==cpi.index) && (inCore==cpi.inCore);};
|
|
bool operator != (const CoredPointIndex& cpi) const {return (index!=cpi.index) || (inCore!=cpi.inCore);};
|
|
};
|
|
template< typename Index > struct CoredEdgeIndex{ CoredPointIndex< Index > idx[2]; };
|
|
|
|
template< typename Index >
|
|
class TriangulationEdge
|
|
{
|
|
public:
|
|
TriangulationEdge( void ){ pIndex[0] = pIndex[1] = tIndex[0] = tIndex[1] = -1; }
|
|
Index pIndex[2] , tIndex[2];
|
|
};
|
|
|
|
template< typename Index >
|
|
class TriangulationTriangle
|
|
{
|
|
public:
|
|
TriangulationTriangle( void ){ eIndex[0] = eIndex[1] = eIndex[2] = -1; }
|
|
Index eIndex[3];
|
|
};
|
|
|
|
template< typename Index >
|
|
struct CoredVertexIndex
|
|
{
|
|
Index idx;
|
|
bool inCore;
|
|
};
|
|
|
|
template< class Vertex , typename Index >
|
|
class CoredCurveData
|
|
{
|
|
public:
|
|
std::vector< Vertex > inCorePoints;
|
|
virtual void resetIterator( void ) = 0;
|
|
|
|
virtual Index addOutOfCorePoint( const Vertex& p ) = 0;
|
|
virtual Index addOutOfCorePoint_s( unsigned int thread , const Vertex& p ) = 0;
|
|
virtual void addEdge_s( unsigned int thread , CoredVertexIndex< Index > v1 , CoredVertexIndex< Index > v2 ) = 0;
|
|
virtual void addEdge_s( unsigned int thread , Index v1 , Index v2 ) = 0;
|
|
|
|
virtual Index nextOutOfCorePoint( Vertex& p )=0;
|
|
virtual Index nextEdge( CoredVertexIndex< Index >& v1 , CoredVertexIndex< Index >& v2 ) = 0;
|
|
|
|
virtual size_t outOfCorePointCount(void)=0;
|
|
virtual size_t edgeCount( void ) = 0;
|
|
};
|
|
|
|
template< class Vertex , typename Index >
|
|
class CoredMeshData
|
|
{
|
|
public:
|
|
virtual ~CoredMeshData( void ){}
|
|
std::vector< Vertex > inCorePoints;
|
|
virtual void resetIterator( void ) = 0;
|
|
|
|
virtual Index addOutOfCorePoint( const Vertex& p ) = 0;
|
|
virtual Index addOutOfCorePoint_s( unsigned int thread , const Vertex& p ) = 0;
|
|
virtual void addPolygon_s( unsigned int thread , const std::vector< CoredVertexIndex< Index > >& vertices ) = 0;
|
|
virtual void addPolygon_s( unsigned int thread , const std::vector< Index >& vertices ) = 0;
|
|
|
|
virtual Index nextOutOfCorePoint( Vertex& p )=0;
|
|
virtual Index nextPolygon( std::vector< CoredVertexIndex< Index > >& vertices ) = 0;
|
|
|
|
virtual size_t outOfCorePointCount( void )=0;
|
|
virtual size_t polygonCount( void ) = 0;
|
|
};
|
|
|
|
template< class Vertex , typename Index >
|
|
class CoredVectorCurveData : public CoredCurveData< Vertex , Index >
|
|
{
|
|
std::vector< Vertex > oocPoints;
|
|
std::vector< std::pair< Index , Index > > edges;
|
|
unsigned int threadIndex;
|
|
Index edgeIndex;
|
|
Index oocPointIndex;
|
|
public:
|
|
CoredVectorCurveData( void );
|
|
|
|
void resetIterator( void );
|
|
|
|
Index addOutOfCorePoint( const Vertex& p );
|
|
Index addOutOfCorePoint_s( unsigned int thread , const Vertex& p );
|
|
void addEdge_s( unsigned int thread , CoredVertexIndex< Index > v1 , CoredVertexIndex< Index > v2 );
|
|
void addEdge_s( unsigned int thread , Index v1 , Index v2 );
|
|
|
|
Index nextOutOfCorePoint( Vertex& p );
|
|
Index nextEdge( CoredVertexIndex< Index > &v1 , CoredVertexIndex< Index > &v2 );
|
|
|
|
size_t outOfCorePointCount( void );
|
|
size_t edgeCount( void );
|
|
};
|
|
template< class Vertex , typename Index >
|
|
class CoredVectorMeshData : public CoredMeshData< Vertex , Index >
|
|
{
|
|
std::vector< Vertex > oocPoints;
|
|
std::vector< std::vector< std::vector< Index > > > polygons;
|
|
unsigned int threadIndex;
|
|
Index polygonIndex;
|
|
Index oocPointIndex;
|
|
public:
|
|
CoredVectorMeshData( void );
|
|
|
|
void resetIterator( void );
|
|
|
|
Index addOutOfCorePoint( const Vertex& p );
|
|
Index addOutOfCorePoint_s( unsigned int thread , const Vertex& p );
|
|
void addPolygon_s( unsigned int thread , const std::vector< CoredVertexIndex< Index > >& vertices );
|
|
void addPolygon_s( unsigned int thread , const std::vector< Index >& vertices );
|
|
|
|
Index nextOutOfCorePoint( Vertex& p );
|
|
Index nextPolygon( std::vector< CoredVertexIndex< Index > >& vertices );
|
|
|
|
size_t outOfCorePointCount( void );
|
|
size_t polygonCount( void );
|
|
};
|
|
class BufferedReadWriteFile
|
|
{
|
|
bool tempFile;
|
|
FILE* _fp;
|
|
char *_buffer , _fileName[1024];
|
|
size_t _bufferIndex , _bufferSize;
|
|
public:
|
|
BufferedReadWriteFile( const char* fileName=NULL , const char* fileHeader="" , unsigned int bufferSize=(1<<20) )
|
|
{
|
|
_bufferIndex = 0;
|
|
_bufferSize = bufferSize;
|
|
if( fileName ) strcpy( _fileName , fileName ) , tempFile = false , _fp = fopen( _fileName , "w+b" );
|
|
else
|
|
{
|
|
if( fileHeader && strlen(fileHeader) ) sprintf( _fileName , "%sXXXXXX" , fileHeader );
|
|
else strcpy( _fileName , "XXXXXX" );
|
|
#ifdef _WIN32
|
|
_mktemp( _fileName );
|
|
_fp = fopen( _fileName , "w+b" );
|
|
#else // !_WIN32
|
|
_fp = fdopen( mkstemp( _fileName ) , "w+b" );
|
|
#endif // _WIN32
|
|
tempFile = true;
|
|
}
|
|
if( !_fp ) ERROR_OUT( "Failed to open file: " , _fileName );
|
|
_buffer = (char*) malloc( _bufferSize );
|
|
}
|
|
~BufferedReadWriteFile( void )
|
|
{
|
|
free( _buffer );
|
|
fclose( _fp );
|
|
if( tempFile ) remove( _fileName );
|
|
}
|
|
bool write( const void* data , size_t size )
|
|
{
|
|
if( !size ) return true;
|
|
const char* _data = (char*) data;
|
|
size_t sz = _bufferSize - _bufferIndex;
|
|
while( sz<=size )
|
|
{
|
|
memcpy( _buffer+_bufferIndex , _data , sz );
|
|
fwrite( _buffer , 1 , _bufferSize , _fp );
|
|
_data += sz;
|
|
size -= sz;
|
|
_bufferIndex = 0;
|
|
sz = _bufferSize;
|
|
}
|
|
if( size )
|
|
{
|
|
memcpy( _buffer+_bufferIndex , _data , size );
|
|
_bufferIndex += size;
|
|
}
|
|
return true;
|
|
}
|
|
bool read( void* data , size_t size )
|
|
{
|
|
if( !size ) return true;
|
|
char *_data = (char*) data;
|
|
size_t sz = _bufferSize - _bufferIndex;
|
|
while( sz<=size )
|
|
{
|
|
if( size && !_bufferSize ) return false;
|
|
memcpy( _data , _buffer+_bufferIndex , sz );
|
|
_bufferSize = fread( _buffer , 1 , _bufferSize , _fp );
|
|
_data += sz;
|
|
size -= sz;
|
|
_bufferIndex = 0;
|
|
if( !size ) return true;
|
|
sz = _bufferSize;
|
|
}
|
|
if( size )
|
|
{
|
|
if( !_bufferSize ) return false;
|
|
memcpy( _data , _buffer+_bufferIndex , size );
|
|
_bufferIndex += size;
|
|
}
|
|
return true;
|
|
}
|
|
void reset( void )
|
|
{
|
|
if( _bufferIndex ) fwrite( _buffer , 1 , _bufferIndex , _fp );
|
|
_bufferIndex = 0;
|
|
fseek( _fp , 0 , SEEK_SET );
|
|
_bufferIndex = 0;
|
|
_bufferSize = fread( _buffer , 1 , _bufferSize , _fp );
|
|
}
|
|
};
|
|
template< class Vertex , typename Index >
|
|
class CoredFileCurveData : public CoredCurveData< Vertex , Index >
|
|
{
|
|
BufferedReadWriteFile *oocPointFile;
|
|
Index oocPoints;
|
|
std::vector< BufferedReadWriteFile* > edgeFiles;
|
|
unsigned int threadIndex;
|
|
public:
|
|
CoredFileCurveData( const char* fileHeader="" );
|
|
~CoredFileCurveData( void );
|
|
|
|
void resetIterator( void );
|
|
|
|
Index addOutOfCorePoint( const Vertex& p );
|
|
Index addOutOfCorePoint_s( unsigned int thread , const Vertex& p );
|
|
void addEdge_s( unsigned int thread , CoredVertexIndex< Index > v1 , CoredVertexIndex< Index > v2 );
|
|
void addEdge_s( unsigned int thread , Index v1 , Index v2 );
|
|
|
|
Index nextOutOfCorePoint( Vertex& p );
|
|
Index nextEdge( CoredVertexIndex< Index > &v1 , CoredVertexIndex< Index > &v2 );
|
|
|
|
size_t outOfCorePointCount( void );
|
|
size_t edgeCount( void );
|
|
};
|
|
|
|
template< class Vertex , typename Index >
|
|
class CoredFileMeshData : public CoredMeshData< Vertex , Index >
|
|
{
|
|
BufferedReadWriteFile *oocPointFile;
|
|
Index oocPoints;
|
|
std::vector< Index > polygons;
|
|
std::vector< BufferedReadWriteFile* > polygonFiles;
|
|
unsigned int threadIndex;
|
|
public:
|
|
CoredFileMeshData( const char* fileHeader="" );
|
|
~CoredFileMeshData( void );
|
|
|
|
void resetIterator( void );
|
|
|
|
Index addOutOfCorePoint( const Vertex& p );
|
|
Index addOutOfCorePoint_s( unsigned int thread , const Vertex& p );
|
|
void addPolygon_s( unsigned int thread , const std::vector< CoredVertexIndex< Index > >& vertices );
|
|
void addPolygon_s( unsigned int thread , const std::vector< Index >& vertices );
|
|
|
|
Index nextOutOfCorePoint( Vertex& p );
|
|
Index nextPolygon( std::vector< CoredVertexIndex< Index > >& vertices );
|
|
|
|
size_t outOfCorePointCount( void );
|
|
size_t polygonCount( void );
|
|
};
|
|
#include "Geometry.inl"
|
|
|
|
#endif // GEOMETRY_INCLUDED
|