lucaspalomodevelop/meshlab
1
0
Fork 0
mirror of https://github.com/lucaspalomodevelop/meshlab.git synced 2026-03-19 19:14:42 +00:00
Code Issues Packages Projects Releases Wiki Activity

Missin files and other updated files for screeend poisson

Browse Source
This commit is contained in:
Paolo Cignoni cignoni 2016-10-21 15:04:42 +00:00
parent 5f601ae2de
commit 7c81c4e0b8
15 changed files with 5105 additions and 153 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/meshlabplugins/filter_screened_poisson/Src/Array.inl
View File

@ -188,7 +188,7 @@ public:
// [WARNING] Chaning szC and szD to size_t causes some really strange behavior.
long long szC = sizeof( C );
long long szD = sizeof( D );
data = (C*)&a[0];
data = (C*)a.data;
min = ( a.minimum() * szD ) / szC;
max = ( a.maximum() * szD ) / szC;
if( min*szC!=a.minimum()*szD || max*szC!=a.maximum()*szD )
@ -306,6 +306,7 @@ public:
return (*this);
}
inline Array& operator ++ ( void ) { return (*this) += 1; }
inline Array operator++( int ){ Array< C > temp = (*this) ; (*this) +=1 ; return temp; }
Array operator - ( int idx ) const { return (*this) + (-idx); }
Array operator - ( long long idx ) const { return (*this) + (-idx); }
Array operator - ( unsigned int idx ) const { return (*this) + (-idx); }
@ -315,6 +316,7 @@ public:
Array& operator -= ( unsigned int idx ) { return (*this) += (-idx); }
Array& operator -= ( unsigned long long idx ) { return (*this) += (-idx); }
Array& operator -- ( void ) { return (*this) -= 1; }
inline Array operator--( int ){ Array< C > temp = (*this) ; (*this) -=1 ; return temp; }
long long operator - ( const Array& a ) const { return ( long long )( data - a.data ); }
void Free( void )
@ -501,6 +503,7 @@ public:
return (*this);
}
inline ConstArray& operator ++ ( void ) { return (*this) += 1; }
inline ConstArray operator++( int ){ ConstArray< C > temp = (*this) ; (*this) +=1 ; return temp; }
ConstArray operator - ( int idx ) const { return (*this) + (-idx); }
ConstArray operator - ( long long idx ) const { return (*this) + (-idx); }
ConstArray operator - ( unsigned int idx ) const { return (*this) + (-idx); }
@ -510,6 +513,7 @@ public:
ConstArray& operator -= ( unsigned int idx ) { return (*this) += (-idx); }
ConstArray& operator -= ( unsigned long long idx ) { return (*this) += (-idx); }
ConstArray& operator -- ( void ) { return (*this) -= 1; }
inline ConstArray operator--( int ){ ConstArray< C > temp = (*this) ; (*this) -=1 ; return temp; }
long long operator - ( const ConstArray& a ) const { return ( long long )( data - a.data ); }
long long operator - ( const Array< C >& a ) const { return ( long long )( data - a.pointer() ); }

60
src/meshlabplugins/filter_screened_poisson/Src/BinaryNode.h
View File

@ -29,9 +29,6 @@ DAMAGE.
#ifndef BINARY_NODE_INCLUDED
#define BINARY_NODE_INCLUDED
#define MSVC_2010_FIX 1
class BinaryNode
{
public:
@ -40,39 +37,34 @@ public:
static inline int CumulativeCenterCount( int maxDepth ) { return (1<<(maxDepth+1))-1; }
static inline int CumulativeCornerCount( int maxDepth ) { return (1<<(maxDepth+1))+maxDepth; }
static inline int CenterIndex( int depth , int offSet ) { return (1<<depth)+offSet-1; }
static inline int CornerIndex( int depth , int offSet ) { return (1<<depth)+offSet+depth; }
static inline int CornerIndex( int depth , int offSet ) { return (1<<depth)+offSet-1+depth; }
static inline int CornerIndex( int maxDepth , int depth , int offSet , int forwardCorner ){ return (offSet+forwardCorner)<<(maxDepth-depth); }
template< class Real > static inline Real CornerIndexPosition(int index,int maxDepth){ return Real(index)/(1<<maxDepth); }
template< class Real > static inline Real Width(int depth){ return Real(1.0/(1<<depth)); }
template< class Real > static inline void CenterAndWidth( int depth , int offset , Real& center , Real& width )
{
width=Real (1.0/(1<<depth) );
center=Real((0.5+offset)*width);
}
template< class Real > static inline Real Width( int depth ){ return Real(1.0/(1<<depth)); }
template< class Real > static inline void CenterAndWidth( int depth , int offset , Real& center , Real& width ){ width = Real (1.0/(1<<depth) ) , center = Real((0.5+offset)*width); }
template< class Real > static inline void CornerAndWidth( int depth , int offset , Real& corner , Real& width ){ width = Real(1.0/(1<<depth) ) , corner = Real(offset*width); }
template< class Real > static inline void CenterAndWidth( int idx , Real& center , Real& width )
{
int depth , offset;
DepthAndOffset( idx , depth , offset );
CenterAndWidth( depth , offset , center , width );
}
static inline void DepthAndOffset( int idx , int& depth , int& offset )
{
int i=idx+1;
#if MSVC_2010_FIX
depth = 0;
#else // !MSVC_2010_FIX
depth = -1;
#endif // MSVC_2010_FIX
while( i )
{
i >>= 1;
depth++;
}
#if MSVC_2010_FIX
depth--;
#endif // MSVC_2010_FIX
offset = ( idx+1 ) - (1<<depth);
}
{
int depth , offset;
CenterDepthAndOffset( idx , depth , offset );
CenterAndWidth( depth , offset , center , width );
}
template< class Real > static inline void CornerAndWidth( int idx , Real& corner , Real& width )
{
int depth , offset;
CornerDepthAndOffset( idx , depth , offset );
CornerAndWidth( depth , offset , corner , width );
}
static inline void CenterDepthAndOffset( int idx , int& depth , int& offset )
{
offset = idx , depth = 0;
while( offset>=(1<<depth) ) offset -= (1<<depth) , depth++;
}
static inline void CornerDepthAndOffset( int idx , int& depth , int& offset )
{
offset = idx , depth = 0;
while( offset>=( (1<<depth)+1 ) ) offset -= ( (1<<depth)+1 ) , depth++;
}
};
#endif // BINARY_NODE_INCLUDED

1
src/meshlabplugins/filter_screened_poisson/Src/CmdLineParser.cpp
View File

@ -166,7 +166,6 @@ char* GetFileExtension(char* fileName){
char* fileNameCopy;
char* ext=NULL;
char* temp;
if(fileName==NULL) return NULL;
fileNameCopy=new char[strlen(fileName)+1];
assert(fileNameCopy);

60
src/meshlabplugins/filter_screened_poisson/Src/FunctionData.inl
View File

@ -1,33 +1,33 @@
/*
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.
*/
//////////////////
// FunctionData //
/*
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.
*/
//////////////////
// FunctionData //
//////////////////
template<int Degree,class Real>
FunctionData<Degree,Real>::FunctionData(void)

1151
src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.Evaluation.inl Normal file
View File

File diff suppressed because it is too large Load Diff

2274
src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.System.inl Normal file
View File

File diff suppressed because it is too large Load Diff

450
src/meshlabplugins/filter_screened_poisson/Src/MultiGridOctreeData.WeightedSamples.inl Normal file
View File

@ -0,0 +1,450 @@
/*
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.
*/
// evaluate the result of splatting along a plane and then evaluating at a point on the plane.
template< int Degree > double GetScaleValue( void )
{
double centerValues[Degree+1];
Polynomial< Degree >::BSplineComponentValues( 0.5 , centerValues );
double scaleValue = 0;
for( int i=0 ; i<=Degree ; i++ ) scaleValue += centerValues[i] * centerValues[i];
return 1./ scaleValue;
}
template< class Real >
template< int WeightDegree >
void Octree< Real >::_addWeightContribution( SparseNodeData< Real , WeightDegree >& densityWeights , TreeOctNode* node , Point3D< Real > position , PointSupportKey< WeightDegree >& weightKey , Real weight )
{
static const double ScaleValue = GetScaleValue< WeightDegree >();
double dx[ DIMENSION ][ PointSupportKey< WeightDegree >::Size ];
typename TreeOctNode::Neighbors< PointSupportKey< WeightDegree >::Size >& neighbors = weightKey.template getNeighbors< true >( node , _NodeInitializer );
densityWeights.reserve( NodeCount() );
Point3D< Real > start;
Real w;
_startAndWidth( node , start , w );
for( int dim=0 ; dim<DIMENSION ; dim++ ) Polynomial< WeightDegree >::BSplineComponentValues( ( position[dim]-start[dim] ) / w , dx[dim] );
weight *= (Real)ScaleValue;
for( int i=0 ; i<PointSupportKey< WeightDegree >::Size ; i++ ) for( int j=0 ; j<PointSupportKey< WeightDegree >::Size ; j++ )
{
double dxdy = dx[0][i] * dx[1][j] * weight;
TreeOctNode** _neighbors = neighbors.neighbors[i][j];
for( int k=0 ; k<PointSupportKey< WeightDegree >::Size ; k++ ) if( _neighbors[k] ) densityWeights[ _neighbors[k] ] += Real( dxdy * dx[2][k] );
}
}
template< class Real >
template< int WeightDegree , class PointSupportKey >
Real Octree< Real >::_getSamplesPerNode( const SparseNodeData< Real , WeightDegree >& densityWeights , const TreeOctNode* node , Point3D< Real > position , PointSupportKey& weightKey ) const
{
Real weight = 0;
double dx[ DIMENSION ][ PointSupportKey::Size ];
const typename PointSupportKey::template Neighbors< PointSupportKey::Size >& neighbors = weightKey.getNeighbors( node );
Point3D< Real > start;
Real w;
_startAndWidth( node , start , w );
for( int dim=0 ; dim<DIMENSION ; dim++ ) Polynomial< WeightDegree >::BSplineComponentValues( ( position[dim]-start[dim] ) / w , dx[dim] );
for( int i=0 ; i<PointSupportKey::Size ; i++ ) for( int j=0 ; j<PointSupportKey::Size ; j++ )
{
double dxdy = dx[0][i] * dx[1][j];
for( int k=0 ; k<PointSupportKey::Size ; k++ ) if( neighbors.neighbors[i][j][k] )
{
const Real* w = densityWeights( neighbors.neighbors[i][j][k] );
if( w ) weight += Real( dxdy * dx[2][k] * (*w) );
}
}
return weight;
}
template< class Real >
template< int WeightDegree , class PointSupportKey >
void Octree< Real >::_getSampleDepthAndWeight( const SparseNodeData< Real , WeightDegree >& densityWeights , const TreeOctNode* node , Point3D< Real > position , PointSupportKey& weightKey , Real& depth , Real& weight ) const
{
const TreeOctNode* temp = node;
weight = _getSamplesPerNode( densityWeights , temp , position , weightKey );
if( weight>=(Real)1. ) depth = Real( _localDepth( temp ) + log( weight ) / log(double(1<<(DIMENSION-1))) );
else
{
Real oldWeight , newWeight;
oldWeight = newWeight = weight;
while( newWeight<(Real)1. && temp->parent )
{
temp=temp->parent;
oldWeight = newWeight;
newWeight = _getSamplesPerNode( densityWeights , temp , position , weightKey );
}
if( newWeight<=0 || oldWeight<=0 )
{
fprintf( stderr , "[WARNING] Octree::_getSampleDepthAndWeight: Weights should be positive: %g %g\n" , oldWeight , newWeight );
depth = weight = 0;
return;
}
else depth = Real( _localDepth( temp ) + log( newWeight ) / log( newWeight / oldWeight ) );
}
weight = Real( pow( double(1<<(DIMENSION-1)) , -double(depth) ) );
}
template< class Real >
template< int WeightDegree , class PointSupportKey >
void Octree< Real >::_getSampleDepthAndWeight( const SparseNodeData< Real , WeightDegree >& densityWeights , Point3D< Real > position , PointSupportKey& weightKey , Real& depth , Real& weight ) const
{
TreeOctNode* temp;
Point3D< Real > myCenter( (Real)0.5 , (Real)0.5 , (Real)0.5 );
Real myWidth = Real( 1. );
// Get the finest node with depth less than or equal to the splat depth that contains the point
temp = _spaceRoot;
while( true )
{
if( !IsActiveNode( temp->children ) ) break;// fprintf( stderr , "[ERROR] Octree::GetSampleDepthAndWeight\n" ) , exit( 0 );
int cIndex = TreeOctNode::CornerIndex( myCenter , position );
if( densityWeights( temp->children + cIndex ) ) temp = temp->children + cIndex;
else break;
myWidth /= 2;
if( cIndex&1 ) myCenter[0] += myWidth/2;
else myCenter[0] -= myWidth/2;
if( cIndex&2 ) myCenter[1] += myWidth/2;
else myCenter[1] -= myWidth/2;
if( cIndex&4 ) myCenter[2] += myWidth/2;
else myCenter[2] -= myWidth/2;
}
return _getSampleDepthAndWeight( densityWeights , temp , position , weightKey , depth , weight );
}
template< class Real >
template< bool CreateNodes , int DataDegree , class V >
void Octree< Real >::_splatPointData( TreeOctNode* node , Point3D< Real > position , V v , SparseNodeData< V , DataDegree >& dataInfo , PointSupportKey< DataDegree >& dataKey )
{
double dx[ DIMENSION ][ PointSupportKey< DataDegree >::Size ];
typename TreeOctNode::Neighbors< PointSupportKey< DataDegree >::Size >& neighbors = dataKey.template getNeighbors< CreateNodes >( node , _NodeInitializer );
Point3D< Real > start;
Real w;
_startAndWidth( node , start , w );
for( int dd=0 ; dd<DIMENSION ; dd++ ) Polynomial< DataDegree >::BSplineComponentValues( ( position[dd]-start[dd] ) / w , dx[dd] );
for( int i=0 ; i<PointSupportKey< DataDegree >::Size ; i++ ) for( int j=0 ; j<PointSupportKey< DataDegree >::Size ; j++ )
{
double dxdy = dx[0][i] * dx[1][j];
for( int k=0 ; k<PointSupportKey< DataDegree >::Size ; k++ )
if( IsActiveNode( neighbors.neighbors[i][j][k] ) )
{
TreeOctNode* _node = neighbors.neighbors[i][j][k];
double dxdydz = dxdy * dx[2][k];
dataInfo[ _node ] += v * (Real)dxdydz;
}
}
}
template< class Real >
template< bool CreateNodes , int WeightDegree , int DataDegree , class V >
Real Octree< Real >::_splatPointData( const SparseNodeData< Real , WeightDegree >& densityWeights , Point3D< Real > position , V v , SparseNodeData< V , DataDegree >& dataInfo , PointSupportKey< WeightDegree >& weightKey , PointSupportKey< DataDegree >& dataKey , LocalDepth minDepth , LocalDepth maxDepth , int dim )
{
double dx;
V _v;
TreeOctNode* temp;
int cnt=0;
double width;
Point3D< Real > myCenter( (Real)0.5 , (Real)0.5 , (Real)0.5 );
Real myWidth = (Real)1.;
temp = _spaceRoot;
while( true )
{
if( !IsActiveNode( temp->children ) ) break;
int cIndex = TreeOctNode::CornerIndex( myCenter , position );
if( densityWeights( temp->children + cIndex ) ) temp = temp->children + cIndex;
else break;
myWidth /= 2;
if( cIndex&1 ) myCenter[0] += myWidth/2;
else myCenter[0] -= myWidth/2;
if( cIndex&2 ) myCenter[1] += myWidth/2;
else myCenter[1] -= myWidth/2;
if( cIndex&4 ) myCenter[2] += myWidth/2;
else myCenter[2] -= myWidth/2;
}
Real weight , depth;
_getSampleDepthAndWeight( densityWeights , temp , position , weightKey , depth , weight );
if( depth<minDepth ) depth = Real(minDepth);
if( depth>maxDepth ) depth = Real(maxDepth);
int topDepth = int(ceil(depth));
dx = 1.0-(topDepth-depth);
if ( topDepth<=minDepth ) topDepth = minDepth , dx = 1;
else if( topDepth> maxDepth ) topDepth = maxDepth , dx = 1;
while( _localDepth( temp )>topDepth ) temp=temp->parent;
while( _localDepth( temp )<topDepth )
{
if( !temp->children ) temp->initChildren( _NodeInitializer );
int cIndex = TreeOctNode::CornerIndex( myCenter , position );
temp = &temp->children[cIndex];
myWidth/=2;
if( cIndex&1 ) myCenter[0] += myWidth/2;
else myCenter[0] -= myWidth/2;
if( cIndex&2 ) myCenter[1] += myWidth/2;
else myCenter[1] -= myWidth/2;
if( cIndex&4 ) myCenter[2] += myWidth/2;
else myCenter[2] -= myWidth/2;
}
width = 1.0 / ( 1<<_localDepth( temp ) );
_v = v * weight / Real( pow( width , dim ) ) * Real( dx );
_splatPointData< CreateNodes >( temp , position , _v , dataInfo , dataKey );
if( fabs(1.0-dx) > EPSILON )
{
dx = Real(1.0-dx);
temp = temp->parent;
width = 1.0 / ( 1<<_localDepth( temp ) );
_v = v * weight / Real( pow( width , dim ) ) * Real( dx );
_splatPointData< CreateNodes >( temp , position , _v , dataInfo , dataKey );
}
return weight;
}
template< class Real >
template< bool CreateNodes , int WeightDegree , int DataDegree , class V >
Real Octree< Real >::_multiSplatPointData( const SparseNodeData< Real , WeightDegree >* densityWeights , TreeOctNode* node , Point3D< Real > position , V v , SparseNodeData< V , DataDegree >& dataInfo , PointSupportKey< WeightDegree >& weightKey , PointSupportKey< DataDegree >& dataKey , int dim )
{
Real _depth , weight;
if( densityWeights ) _getSampleDepthAndWeight( *densityWeights , position , weightKey , _depth , weight );
else weight = (Real)1.;
V _v = v * weight;
double dx[ DIMENSION ][ PointSupportKey< DataDegree >::Size ];
dataKey.template getNeighbors< CreateNodes >( node , _NodeInitializer );
for( TreeOctNode* _node=node ; _localDepth( _node )>=0 ; _node=_node->parent )
{
V __v = _v * (Real)pow( 1<<_localDepth( _node ) , dim );
Point3D< Real > start;
Real w;
_startAndWidth( _node , start , w );
for( int dd=0 ; dd<DIMENSION ; dd++ ) Polynomial< DataDegree >::BSplineComponentValues( ( position[dd]-start[dd] ) / w , dx[dd] );
typename TreeOctNode::Neighbors< PointSupportKey< DataDegree >::Size >& neighbors = dataKey.neighbors[ _localToGlobal( _localDepth( _node ) ) ];
for( int i=0 ; i<PointSupportKey< DataDegree >::Size ; i++ ) for( int j=0 ; j<PointSupportKey< DataDegree >::Size ; j++ )
{
double dxdy = dx[0][i] * dx[1][j];
for( int k=0 ; k<PointSupportKey< DataDegree >::Size ; k++ )
if( IsActiveNode( neighbors.neighbors[i][j][k] ) )
{
TreeOctNode* _node = neighbors.neighbors[i][j][k];
double dxdydz = dxdy * dx[2][k];
dataInfo[ _node ] += __v * (Real)dxdydz;
}
}
}
return weight;
}
template< class Real >
template< class V , int DataDegree , BoundaryType BType , class Coefficients >
V Octree< Real >::_evaluate( const Coefficients& coefficients , Point3D< Real > p , const BSplineData< DataDegree , BType >& bsData , const ConstPointSupportKey< DataDegree >& dataKey ) const
{
V value = V(0);
for( int d=_localToGlobal( 0 ) ; d<=dataKey.depth() ; d++ )
{
double dx[ DIMENSION ][ PointSupportKey< DataDegree >::Size ];
memset( dx , 0 , sizeof( double ) * DIMENSION * PointSupportKey< DataDegree >::Size );
{
const TreeOctNode* n = dataKey.neighbors[d].neighbors[ PointSupportKey< DataDegree >::LeftRadius ][ PointSupportKey< DataDegree >::LeftRadius ][ PointSupportKey< DataDegree >::LeftRadius ];
if( !n ) fprintf( stderr , "[ERROR] Point is not centered on a node\n" ) , exit( 0 );
int fIdx[3];
functionIndex< DataDegree , BType >( n , fIdx );
int fStart , fEnd;
BSplineData< DataDegree , BType >::FunctionSpan( _localDepth( n ) , fStart , fEnd );
for( int dd=0 ; dd<DIMENSION ; dd++ ) for( int i=-PointSupportKey< DataDegree >::LeftRadius ; i<=PointSupportKey< DataDegree >::RightRadius ; i++ )
if( fIdx[dd]+i>=fStart && fIdx[dd]+i<fEnd ) dx[dd][i] = bsData.baseBSplines[ fIdx[dd]+i ][ -i+PointSupportKey< DataDegree >::RightRadius ]( p[dd] );
}
for( int i=0 ; i<PointSupportKey< DataDegree >::Size ; i++ ) for( int j=0 ; j<PointSupportKey< DataDegree >::Size ; j++ ) for( int k=0 ; k<PointSupportKey< DataDegree >::Size ; k++ )
{
const TreeOctNode* n = dataKey.neighbors[d].neighbors[i][j][k];
if( isValidFEMNode< DataDegree , BType >( n ) )
{
const V* v = coefficients( n );
if( v ) value += (*v) * (Real) ( dx[0][i] * dx[1][j] * dx[2][k] );
}
}
}
return value;
}
template< class Real >
template< class V , int DataDegree , BoundaryType BType >
Pointer( V ) Octree< Real >::voxelEvaluate( const DenseNodeData< V , DataDegree >& coefficients , int& res , Real isoValue , LocalDepth depth , bool primal )
{
int begin , end , dim;
if( depth<=0 || depth>_maxDepth ) depth = _maxDepth;
// Initialize the coefficients at the coarsest level
Pointer( V ) _coefficients = NullPointer( V );
{
LocalDepth d = 0;
begin = _BSplineBegin< DataDegree , BType >( d ) , end = _BSplineEnd< DataDegree , BType >( d ) , dim = end - begin;
_coefficients = NewPointer< V >( dim * dim * dim );
memset( _coefficients , 0 , sizeof( V ) * dim * dim * dim );
#pragma omp parallel for num_threads( threads )
for( int i=_sNodesBegin(d) ; i<_sNodesEnd(d) ; i++ ) if( !_outOfBounds< DataDegree , BType >( _sNodes.treeNodes[i] ) )
{
LocalDepth _d ; LocalOffset _off;
_localDepthAndOffset( _sNodes.treeNodes[i] , _d , _off );
_off[0] -= begin , _off[1] -= begin , _off[2] -= begin;
_coefficients[ _off[0] + _off[1]*dim + _off[2]*dim*dim ] = coefficients[i];
}
}
// Up-sample and add in the existing coefficients
for( LocalDepth d=1 ; d<=depth ; d++ )
{
begin = _BSplineBegin< DataDegree , BType >( d ) , end = _BSplineEnd< DataDegree , BType >( d ) , dim = end - begin;
Pointer( V ) __coefficients = NewPointer< V >( dim * dim *dim );
memset( __coefficients , 0 , sizeof( V ) * dim * dim * dim );
#pragma omp parallel for num_threads( threads )
for( int i=_sNodesBegin(d) ; i<_sNodesEnd(d) ; i++ ) if( !_outOfBounds< DataDegree , BType >( _sNodes.treeNodes[i] ) )
{
LocalDepth _d ; LocalOffset _off;
_localDepthAndOffset( _sNodes.treeNodes[i] , _d , _off );
_off[0] -= begin , _off[1] -= begin , _off[2] -= begin;
__coefficients[ _off[0] + _off[1]*dim + _off[2]*dim*dim ] = coefficients[i];
}
_UpSample< V , DataDegree , BType >( d , ( ConstPointer(V) )_coefficients , __coefficients , threads );
DeletePointer( _coefficients );
_coefficients = __coefficients;
}
res = 1<<depth;
if( primal ) res++;
Pointer( V ) values = NewPointer< V >( res*res*res );
memset( values , 0 , sizeof(V)*res*res*res );
if( primal )
{
// evaluate at the cell corners
typename BSplineEvaluationData< DataDegree , BType >::CornerEvaluator::Evaluator evaluator;
BSplineEvaluationData< DataDegree , BType >::SetCornerEvaluator( evaluator , depth );
#pragma omp parallel for num_threads( threads )
for( int k=0 ; k<res ; k++ ) for( int j=0 ; j<res ; j++ ) for( int i=0 ; i<res ; i++ )
{
V value = values[ i + j*res + k*res*res ];
for( int kk=-BSplineSupportSizes< DataDegree >::CornerEnd ; kk<=-BSplineSupportSizes< DataDegree >::CornerStart ; kk++ ) if( k+kk>=begin && k+kk<end )
for( int jj=-BSplineSupportSizes< DataDegree >::CornerEnd ; jj<=-BSplineSupportSizes< DataDegree >::CornerStart ; jj++ ) if( j+jj>=begin && j+jj<end )
{
double weight = evaluator.value( k+kk , k , false ) * evaluator.value( j+jj , j , false );
int idx = (j+jj-begin)*dim + (k+kk-begin)*dim*dim;
for( int ii=-BSplineSupportSizes< DataDegree >::CornerEnd ; ii<=-BSplineSupportSizes< DataDegree >::CornerStart ; ii++ ) if( i+ii>=begin && i+ii<end )
value += _coefficients[ i + ii - begin + idx ] * Real( weight * evaluator.value( i + ii , i , false ) );
}
values[ i + j*res + k*res*res ] = value;
}
}
else
{
// evaluate at the cell centers
typename BSplineEvaluationData< DataDegree , BType >::CenterEvaluator::Evaluator evaluator;
BSplineEvaluationData< DataDegree , BType >::SetCenterEvaluator( evaluator , depth );
#pragma omp parallel for num_threads( threads )
for( int k=0 ; k<res ; k++ ) for( int j=0 ; j<res ; j++ ) for( int i=0 ; i<res ; i++ )
{
V& value = values[ i + j*res + k*res*res ];
for( int kk=-BSplineSupportSizes< DataDegree >::SupportEnd ; kk<=-BSplineSupportSizes< DataDegree >::SupportStart ; kk++ ) if( k+kk>=begin && k+kk<end )
for( int jj=-BSplineSupportSizes< DataDegree >::SupportEnd ; jj<=-BSplineSupportSizes< DataDegree >::SupportStart ; jj++ ) if( j+jj>=begin && j+jj<end )
{
double weight = evaluator.value( k+kk , k , false ) * evaluator.value( j+jj , j , false );
int idx = (j+jj-begin)*dim + (k+kk-begin)*dim*dim;
for( int ii=-BSplineSupportSizes< DataDegree >::SupportEnd ; ii<=-BSplineSupportSizes< DataDegree >::SupportStart ; ii++ ) if( i+ii>=begin && i+ii<end )
value += _coefficients[ i + ii - begin + idx ] * Real( weight * evaluator.value( i+ii , i , false ) );
}
}
}
memoryUsage();
DeletePointer( _coefficients );
for( int i=0 ; i<res*res*res ; i++ ) values[i] -= isoValue;
return values;
}
template< class Real >
template< int FEMDegree , BoundaryType BType >
SparseNodeData< Real , 0 > Octree< Real >::leafValues( const DenseNodeData< Real , FEMDegree >& coefficients ) const
{
SparseNodeData< Real , 0 > values;
DenseNodeData< Real , FEMDegree > _coefficients( _sNodesEnd(_maxDepth-1) );
memset( &_coefficients[0] , 0 , sizeof(Real)*_sNodesEnd(_maxDepth-1) );
for( int i=_sNodes.begin( _localToGlobal( 0 ) ) ; i<_sNodesEnd(_maxDepth-1) ; i++ ) _coefficients[i] = coefficients[i];
for( LocalDepth d=1 ; d<_maxDepth ; d++ ) _upSample( d , _coefficients );
for( LocalDepth d=_maxDepth ; d>=0 ; d-- )
{
_Evaluator< FEMDegree , BType > evaluator;
evaluator.set( d );
std::vector< ConstPointSupportKey< FEMDegree > > neighborKeys( std::max< int >( 1 , threads ) );
for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( _localToGlobal( d ) );
for( int i=_sNodesBegin(d) ; i<_sNodesEnd(d) ; i++ ) if( _isValidSpaceNode( _sNodes.treeNodes[i] ) )
{
ConstPointSupportKey< FEMDegree >& neighborKey = neighborKeys[ omp_get_thread_num() ];
TreeOctNode* node = _sNodes.treeNodes[i];
if( !IsActiveNode( node->children ) )
{
neighborKey.getNeighbors( node );
bool isInterior = _IsInteriorlySupported< FEMDegree >( node->parent );
values[ node ] = _getCenterValue( neighborKey , node , coefficients , _coefficients , evaluator , isInterior );
}
}
}
return values;
}
template< class Real >
template< int FEMDegree , BoundaryType BType >
SparseNodeData< Point3D< Real > , 0 > Octree< Real >::leafGradients( const DenseNodeData< Real , FEMDegree >& coefficients ) const
{
SparseNodeData< Point3D< Real > , 0 > gradients;
DenseNodeData< Real , FEMDegree > _coefficients( _sNodesEnd(_maxDepth-1 ) );
memset( &_coefficients[0] , 0 , sizeof(Real)*_sNodesEnd(_maxDepth-1) );
for( int i=_sNodesBegin(0) ; i<_sNodesEnd(_maxDepth-1) ; i++ ) _coefficients[i] = coefficients[i];
for( LocalDepth d=1 ; d<_maxDepth ; d++ ) _upSample( d , _coefficients );
for( LocalDepth d=_maxDepth ; d>=0 ; d-- )
{
_Evaluator< FEMDegree , BType > evaluator;
evaluator.set( d );
std::vector< ConstPointSupportKey< FEMDegree > > neighborKeys( std::max< int >( 1 , threads ) );
for( size_t i=0 ; i<neighborKeys.size() ; i++ ) neighborKeys[i].set( _localToGlobal( d ) );
for( int i=_sNodesBegin(d) ; i<_sNodesEnd(d) ; i++ ) if( _isValidSpaceNode( _sNodes.treeNodes[i] ) )
{
ConstPointSupportKey< FEMDegree >& neighborKey = neighborKeys[ omp_get_thread_num() ];
TreeOctNode* node = _sNodes.treeNodes[i];
if( !IsActiveNode( node->children ) )
{
neighborKey.getNeighbors( node );
bool isInterior = _IsInteriorlySupported< FEMDegree >( node->parent );
gradients[ node ] = _getCenterValueAndGradient( neighborKey , node , coefficients , _coefficients , evaluator , isInterior ).second;
}
}
}
return gradients;
}

51
src/meshlabplugins/filter_screened_poisson/Src/MyTime.h Normal file
View File

@ -0,0 +1,51 @@
/*
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 MY_TIME_INCLUDED
#define MY_TIME_INCLUDED
#include <string.h>
#include <sys/timeb.h>
#ifndef WIN32
#include <sys/time.h>
#endif // WIN32
inline double Time( void )
{
#ifdef WIN32
struct _timeb t;
_ftime( &t );
return double( t.time ) + double( t.millitm ) / 1000.0;
#else // WIN32
struct timeval t;
gettimeofday( &t , NULL );
return t.tv_sec + double( t.tv_usec ) / 1000000;
#endif // WIN32
}
#endif // MY_TIME_INCLUDED

56
src/meshlabplugins/filter_screened_poisson/Src/PPolynomial.h
View File

@ -30,57 +30,56 @@ DAMAGE.
#define P_POLYNOMIAL_INCLUDED
#include <vector>
#include "Polynomial.h"
#include "Array.h"
template<int Degree>
class StartingPolynomial{
template< int Degree >
class StartingPolynomial
{
public:
Polynomial<Degree> p;
Polynomial< Degree > p;
double start;
template<int Degree2>
StartingPolynomial<Degree+Degree2> operator * (const StartingPolynomial<Degree2>& p) const;
StartingPolynomial scale(double s) const;
StartingPolynomial shift(double t) const;
int operator < (const StartingPolynomial& sp) const;
static int Compare(const void* v1,const void* v2);
template< int Degree2 >
StartingPolynomial< Degree+Degree2 > operator * ( const StartingPolynomial< Degree2 >& p ) const;
StartingPolynomial scale( double s ) const;
StartingPolynomial shift( double t ) const;
int operator < ( const StartingPolynomial& sp ) const;
static int Compare( const void* v1 , const void* v2 );
};
template<int Degree>
template< int Degree >
class PPolynomial
{
public:
size_t polyCount;
StartingPolynomial<Degree>* polys;
Pointer( StartingPolynomial< Degree > ) polys;
PPolynomial(void);
PPolynomial(const PPolynomial<Degree>& p);
~PPolynomial(void);
PPolynomial( void );
PPolynomial( const PPolynomial<Degree>& p );
~PPolynomial( void );
PPolynomial& operator = (const PPolynomial& p);
PPolynomial& operator = ( const PPolynomial& p );
int size(void) const;
int size( void ) const;
void set( size_t size );
// Note: this method will sort the elements in sps
void set( StartingPolynomial<Degree>* sps , int count );
void set( Pointer( StartingPolynomial<Degree> ) sps , int count );
void reset( size_t newSize );
PPolynomial& compress( double delta=0. );
double operator()( double t ) const;
double integral( double tMin , double tMax ) const;
double Integral( void ) const;
template<int Degree2>
PPolynomial<Degree>& operator = (const PPolynomial<Degree2>& p);
template< int Degree2 > PPolynomial< Degree >& operator = ( const PPolynomial< Degree2 >& p );
PPolynomial operator + (const PPolynomial& p) const;
PPolynomial operator - (const PPolynomial& p) const;
PPolynomial operator + ( const PPolynomial& p ) const;
PPolynomial operator - ( const PPolynomial& p ) const;
template<int Degree2>
PPolynomial<Degree+Degree2> operator * (const Polynomial<Degree2>& p) const;
template<int Degree2>
PPolynomial<Degree+Degree2> operator * (const PPolynomial<Degree2>& p) const;
template< int Degree2 > PPolynomial< Degree+Degree2 > operator * ( const Polynomial< Degree2 >& p ) const;
template< int Degree2 > PPolynomial< Degree+Degree2 > operator * ( const PPolynomial< Degree2 >& p) const;
PPolynomial& operator += ( double s );
@ -92,15 +91,16 @@ public:
PPolynomial operator * ( double s ) const;
PPolynomial operator / ( double s ) const;
PPolynomial& addScaled(const PPolynomial& poly,double scale);
PPolynomial& addScaled( const PPolynomial& poly , double scale );
PPolynomial scale( double s ) const;
PPolynomial shift( double t ) const;
PPolynomial reflect( double r=0. ) const;
PPolynomial< Degree-1 > derivative(void) const;
PPolynomial< Degree+1 > integral(void) const;
void getSolutions(double c,std::vector<double>& roots,double EPS,double min=-DBL_MAX,double max=DBL_MAX) const;
void getSolutions( double c , std::vector< double >& roots , double EPS , double min=-DBL_MAX , double max=DBL_MAX ) const;
void printnl( void ) const;

283
src/meshlabplugins/filter_screened_poisson/Src/Ply.h
View File

@ -33,6 +33,8 @@
#ifndef __PLY_H__
#define __PLY_H__
#define USE_PLY_WRAPPER 1
#ifndef WIN32
#define _strdup strdup
#endif
@ -239,12 +241,22 @@ static PlyProperty face_props[] =
{ _strdup( "vertex_indices" ) , PLY_INT , PLY_INT , offsetof( PlyFace , vertices ) , 1 , PLY_UCHAR, PLY_UCHAR , offsetof(PlyFace,nr_vertices) },
};
///////////////////
// PlyVertexType //
///////////////////
// The "Wrapper" class indicates the class to cast to/from in order to support linear operations.
template< class Real >
class PlyVertex
{
public:
const static int Components=3;
static PlyProperty Properties[];
typedef PlyVertex Wrapper;
const static int ReadComponents=3;
const static int WriteComponents=3;
static PlyProperty ReadProperties[];
static PlyProperty WriteProperties[];
Point3D< Real > point;
@ -260,7 +272,13 @@ public:
template< class _Real > PlyVertex& operator /= ( _Real s ) { point /= s ; return *this; }
};
template< class Real , class _Real > PlyVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyVertex< Real > v ) { return PlyVertex< Real >( xForm * v.point ); }
template< class Real > PlyProperty PlyVertex< Real >::Properties[]=
template< class Real > PlyProperty PlyVertex< Real >::ReadProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 }
};
template< class Real > PlyProperty PlyVertex< Real >::WriteProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
@ -270,8 +288,12 @@ template< class Real >
class PlyValueVertex
{
public:
const static int Components=4;
static PlyProperty Properties[];
typedef PlyValueVertex Wrapper;
const static int ReadComponents=4;
const static int WriteComponents=4;
static PlyProperty ReadProperties[];
static PlyProperty WriteProperties[];
Point3D<Real> point;
Real value;
@ -288,8 +310,14 @@ public:
template< class _Real > PlyValueVertex& operator /= ( _Real s ) { point /= s , value /= Real(s) ; return *this; }
};
template< class Real , class _Real > PlyValueVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyValueVertex< Real > v ) { return PlyValueVertex< Real >( xForm * v.point , v.value ); }
template< class Real >
PlyProperty PlyValueVertex< Real >::Properties[]=
template< class Real > PlyProperty PlyValueVertex< Real >::ReadProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "value" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , value ) ) , 0 , 0 , 0 , 0 }
};
template< class Real > PlyProperty PlyValueVertex< Real >::WriteProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
@ -300,8 +328,12 @@ template< class Real >
class PlyOrientedVertex
{
public:
const static int Components=6;
static PlyProperty Properties[];
typedef PlyOrientedVertex Wrapper;
const static int ReadComponents=6;
const static int WriteComponents=6;
static PlyProperty ReadProperties[];
static PlyProperty WriteProperties[];
Point3D<Real> point , normal;
@ -317,8 +349,16 @@ public:
template< class _Real > PlyOrientedVertex& operator /= ( _Real s ) { point /= s , normal /= s ; return *this; }
};
template< class Real , class _Real > PlyOrientedVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyOrientedVertex< Real > v ) { return PlyOrientedVertex< Real >( xForm * v.point , xForm.inverse().transpose() * v.normal ); }
template< class Real >
PlyProperty PlyOrientedVertex< Real >::Properties[]=
template< class Real > PlyProperty PlyOrientedVertex< Real >::ReadProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "nx" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , normal.coords[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "ny" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , normal.coords[1] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "nz" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , normal.coords[2] ) ) , 0 , 0 , 0 , 0 }
};
template< class Real > PlyProperty PlyOrientedVertex< Real >::WriteProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyOrientedVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
@ -331,19 +371,61 @@ template< class Real >
class PlyColorVertex
{
public:
const static int Components=6;
static PlyProperty Properties[];
struct _PlyColorVertex
{
Point3D< Real > point , color;
_PlyColorVertex( void ) { ; }
_PlyColorVertex( Point3D< Real > p , Point3D< Real > c ) : point(p) , color(c) { ; }
_PlyColorVertex( PlyColorVertex< Real > p ){ point = p.point ; for( int c=0 ; c<3 ; c++ ) color[c] = (Real) p.color[c]; }
operator PlyColorVertex< Real > ()
{
PlyColorVertex< Real > p;
p.point = point;
for( int c=0 ; c<3 ; c++ ) p.color[c] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( color[c]+0.5 ) ) );
return p;
}
Point3D<Real> point;
_PlyColorVertex operator + ( _PlyColorVertex p ) const { return _PlyColorVertex( point+p.point , color+p.color ); }
_PlyColorVertex operator - ( _PlyColorVertex p ) const { return _PlyColorVertex( point-p.value , color-p.color ); }
template< class _Real > _PlyColorVertex operator * ( _Real s ) const { return _PlyColorVertex( point*s , color*s ); }
template< class _Real > _PlyColorVertex operator / ( _Real s ) const { return _PlyColorVertex( point/s , color/s ); }
_PlyColorVertex& operator += ( _PlyColorVertex p ) { point += p.point , color += p.color ; return *this; }
_PlyColorVertex& operator -= ( _PlyColorVertex p ) { point -= p.point , color -= p.color ; return *this; }
template< class _Real > _PlyColorVertex& operator *= ( _Real s ) { point *= s , color *= s ; return *this; }
template< class _Real > _PlyColorVertex& operator /= ( _Real s ) { point /= s , color /= s ; return *this; }
};
typedef _PlyColorVertex Wrapper;
const static int ReadComponents=9;
const static int WriteComponents=6;
static PlyProperty ReadProperties[];
static PlyProperty WriteProperties[];
Point3D< Real > point;
unsigned char color[3];
operator Point3D<Real>& () {return point;}
operator const Point3D<Real>& () const {return point;}
PlyColorVertex(void) {point.coords[0]=point.coords[1]=point.coords[2]=0,color[0]=color[1]=color[2]=0;}
PlyColorVertex(const Point3D<Real>& p) {point=p;}
operator Point3D< Real >& (){ return point; }
operator const Point3D< Real >& () const { return point; }
PlyColorVertex( void ) { point.coords[0] = point.coords[1] = point.coords[2] = 0 , color[0] = color[1] = color[2] = 0; }
PlyColorVertex( const Point3D<Real>& p ) { point=p; }
PlyColorVertex( const Point3D< Real >& p , const unsigned char c[3] ) { point = p , color[0] = c[0] , color[1] = c[1] , color[2] = c[2]; }
};
template< class Real >
PlyProperty PlyColorVertex< Real >::Properties[]=
template< class Real , class _Real > PlyColorVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyColorVertex< Real > v ) { return PlyColorVertex< Real >( xForm * v.point , v.color ); }
template< class Real > PlyProperty PlyColorVertex< Real >::ReadProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "z" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "red" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[1] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "blue" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[2] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "r" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "g" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[1] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "b" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[2] ) ) , 0 , 0 , 0 , 0 }
};
template< class Real > PlyProperty PlyColorVertex< Real >::WriteProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >(), int( offsetof( PlyColorVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 },
@ -352,6 +434,77 @@ PlyProperty PlyColorVertex< Real >::Properties[]=
{ _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[1] ) ) , 0 , 0 , 0 , 0 },
{ _strdup( "blue" ) , PLYType< unsigned char >() , PLYType< unsigned char >(), int( offsetof( PlyColorVertex , color[2] ) ) , 0 , 0 , 0 , 0 }
};
template< class Real >
class PlyColorAndValueVertex
{
public:
struct _PlyColorAndValueVertex
{
Point3D< Real > point , color;
Real value;
_PlyColorAndValueVertex( void ) : value(0) { ; }
_PlyColorAndValueVertex( Point3D< Real > p , Point3D< Real > c , Real v ) : point(p) , color(c) , value(v) { ; }
_PlyColorAndValueVertex( PlyColorAndValueVertex< Real > p ){ point = p.point ; for( int c=0 ; c<3 ; c++ ) color[c] = (Real) p.color[c] ; value = p.value; }
operator PlyColorAndValueVertex< Real > ()
{
PlyColorAndValueVertex< Real > p;
p.point = point;
for( int c=0 ; c<3 ; c++ ) p.color[c] = (unsigned char)std::max< int >( 0 , std::min< int >( 255 , (int)( color[c]+0.5 ) ) );
p.value = value;
return p;
}
_PlyColorAndValueVertex operator + ( _PlyColorAndValueVertex p ) const { return _PlyColorAndValueVertex( point+p.point , color+p.color , value+p.value ); }
_PlyColorAndValueVertex operator - ( _PlyColorAndValueVertex p ) const { return _PlyColorAndValueVertex( point-p.value , color-p.color , value+p.value ); }
template< class _Real > _PlyColorAndValueVertex operator * ( _Real s ) const { return _PlyColorAndValueVertex( point*s , color*s , value*s ); }
template< class _Real > _PlyColorAndValueVertex operator / ( _Real s ) const { return _PlyColorAndValueVertex( point/s , color/s , value/s ); }
_PlyColorAndValueVertex& operator += ( _PlyColorAndValueVertex p ) { point += p.point , color += p.color , value += p.value ; return *this; }
_PlyColorAndValueVertex& operator -= ( _PlyColorAndValueVertex p ) { point -= p.point , color -= p.color , value -= p.value ; return *this; }
template< class _Real > _PlyColorAndValueVertex& operator *= ( _Real s ) { point *= s , color *= s , value *= (Real)s ; return *this; }
template< class _Real > _PlyColorAndValueVertex& operator /= ( _Real s ) { point /= s , color /= s , value /= (Real)s ; return *this; }
};
typedef _PlyColorAndValueVertex Wrapper;
const static int ReadComponents=10;
const static int WriteComponents=7;
static PlyProperty ReadProperties[];
static PlyProperty WriteProperties[];
Point3D< Real > point;
unsigned char color[3];
Real value;
operator Point3D< Real >& (){ return point; }
operator const Point3D< Real >& () const { return point; }
PlyColorAndValueVertex( void ) { point.coords[0] = point.coords[1] = point.coords[2] = (Real)0 , color[0] = color[1] = color[2] = 0 , value = (Real)0; }
PlyColorAndValueVertex( const Point3D< Real >& p ) { point=p; }
PlyColorAndValueVertex( const Point3D< Real >& p , const unsigned char c[3] , Real v) { point = p , color[0] = c[0] , color[1] = c[1] , color[2] = c[2] , value = v; }
};
template< class Real , class _Real > PlyColorAndValueVertex< Real > operator * ( XForm4x4< _Real > xForm , PlyColorAndValueVertex< Real > v ) { return PlyColorAndValueVertex< Real >( xForm * v.point , v.color , v.value ); }
template< class Real > PlyProperty PlyColorAndValueVertex< Real >::ReadProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "value" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , value ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "red" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[0] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[1] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "blue" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[2] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "r" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[0] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "g" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[1] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "b" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[2] ) ) , 0 , 0 , 0 , 0 }
};
template< class Real > PlyProperty PlyColorAndValueVertex< Real >::WriteProperties[]=
{
{ _strdup( "x" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[0] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "y" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[1] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "z" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , point.coords[2] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "value" ) , PLYType< Real >() , PLYType< Real >() , int( offsetof( PlyColorAndValueVertex , value ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "red" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[0] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "green" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[1] ) ) , 0 , 0 , 0 , 0 } ,
{ _strdup( "blue" ) , PLYType< unsigned char >() , PLYType< unsigned char >() , int( offsetof( PlyColorAndValueVertex , color[2] ) ) , 0 , 0 , 0 , 0 }
};
template< class Vertex , class Real >
int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_type , const Point3D< float >& translate , float scale , char** comments=NULL , int commentNum=0 , XForm4x4< Real > xForm=XForm4x4< Real >::Identity() );
@ -359,6 +512,94 @@ int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file
template< class Vertex , class Real >
int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_type , char** comments=NULL , int commentNum=0 , XForm4x4< Real > xForm=XForm4x4< Real >::Identity() );
inline bool PlyReadHeader( char* fileName , PlyProperty* properties , int propertyNum , bool* readFlags , int& file_type )
{
int nr_elems;
char **elist;
float version;
PlyFile* ply;
char* elem_name;
int num_elems;
int nr_props;
PlyProperty** plist;
ply = ply_open_for_reading( fileName , &nr_elems , &elist , &file_type , &version );
if( !ply ) return false;
for( int i=0 ; i<nr_elems ; i++ )
{
elem_name = elist[i];
plist = ply_get_element_description( ply , elem_name , &num_elems , &nr_props );
if( !plist )
{
for( int i=0 ; i<nr_elems ; i++ )
{
free( ply->elems[i]->name );
free( ply->elems[i]->store_prop );
for( int j=0 ; j<ply->elems[i]->nprops ; j++ )
{
free( ply->elems[i]->props[j]->name );
free( ply->elems[i]->props[j] );
}
free( ply->elems[i]->props );
}
for( int i=0 ; i<nr_elems ; i++ ) free( ply->elems[i] );
free( ply->elems );
for( int i=0 ; i<ply->num_comments ; i++ ) free( ply->comments[i] );
free( ply->comments );
for( int i=0 ; i<ply->num_obj_info ; i++ ) free( ply->obj_info[i] );
free( ply->obj_info );
ply_free_other_elements( ply->other_elems );
for( int i=0 ; i<nr_elems ; i++ ) free( elist[i] );
free( elist );
ply_close( ply );
return 0;
}
if( equal_strings( "vertex" , elem_name ) )
for( int i=0 ; i<propertyNum ; i++ )
if( readFlags ) readFlags[i] = ply_get_property( ply , elem_name , &properties[i] )!=0;
for( int j=0 ; j<nr_props ; j++ )
{
free( plist[j]->name );
free( plist[j] );
}
free( plist );
} // for each type of element
for( int i=0 ; i<nr_elems ; i++ )
{
free( ply->elems[i]->name );
free( ply->elems[i]->store_prop );
for( int j=0 ; j<ply->elems[i]->nprops ; j++ )
{
free( ply->elems[i]->props[j]->name );
free( ply->elems[i]->props[j] );
}
if( ply->elems[i]->props && ply->elems[i]->nprops ) free(ply->elems[i]->props);
}
for( int i=0 ; i<nr_elems ; i++ ) free(ply->elems[i]);
free( ply->elems) ;
for( int i=0 ; i<ply->num_comments ; i++ ) free( ply->comments[i] );
free( ply->comments );
for( int i=0 ; i<ply->num_obj_info ; i++ ) free( ply->obj_info[i] );
free( ply->obj_info );
ply_free_other_elements(ply->other_elems);
for( int i=0 ; i<nr_elems ; i++ ) free( elist[i] );
free( elist );
ply_close( ply );
return true;
}
inline bool PlyReadHeader( char* fileName , PlyProperty* properties , int propertyNum , bool* readFlags )
{
int file_type;
return PlyReadHeader( fileName , properties , propertyNum , readFlags , file_type );
}
template<class Vertex>
int PlyReadPolygons(char* fileName,
std::vector<Vertex>& vertices,std::vector<std::vector<int> >& polygons,
@ -626,7 +867,7 @@ int PlyWritePolygons( char* fileName , CoredMeshData< Vertex >* mesh , int file_
// describe vertex and face properties
//
ply_element_count( ply , "vertex" , nr_vertices );
for( int i=0 ; i<Vertex::Components ; i++ ) ply_describe_property( ply , "vertex" , &Vertex::Properties[i] );
for( int i=0 ; i<Vertex::WriteComponents ; i++ ) ply_describe_property( ply , "vertex" , &Vertex::WriteProperties[i] );
ply_element_count( ply , "face" , nr_faces );
ply_describe_property( ply , "face" , &face_props[0] );

4
src/meshlabplugins/filter_screened_poisson/Src/PlyFile.cpp
View File

@ -259,6 +259,7 @@ PlyFile *ply_open_for_writing(
/* open the file for writing */
fp = fopen (name, "wb");
free(name);
if (fp == NULL) {
return (NULL);
}
@ -878,6 +879,7 @@ Open a polygon file for reading.
/* open the file for reading */
fp = fopen (name, "rb");
free(name);
if (fp == NULL)
return (NULL);
@ -892,7 +894,6 @@ Open a polygon file for reading.
/* return a pointer to the file's information */
free(name);
return (plyfile);
}
@ -1873,6 +1874,7 @@ Read an element from a binary file.
/* read in a line */
result = fgets (str, BIG_STRING, fp);
if (result == NULL) {
free(words);
*nwords = 0;
*orig_line = NULL;
return (NULL);

9
src/meshlabplugins/filter_screened_poisson/Src/Polynomial.h
View File

@ -29,10 +29,13 @@ DAMAGE.
#ifndef POLYNOMIAL_INCLUDED
#define POLYNOMIAL_INCLUDED
#define NEW_POLYNOMIAL_CODE 1
#include <vector>
template< int Degree >
class Polynomial{
class Polynomial
{
public:
double coefficients[Degree+1];
@ -86,7 +89,11 @@ public:
void getSolutions(double c,std::vector<double>& roots,double EPS) const;
int getSolutions( double c , double* roots , double EPS ) const;
// [NOTE] Both of these methods define the indexing according to DeBoor's algorithm, so that
// Polynomial< Degree >BSplineComponent( 0 )( 1.0 )=0 for all Degree>0.
static Polynomial BSplineComponent( int i );
static void BSplineComponentValues( double x , double* values );
static void BinomialCoefficients( int bCoefficients[Degree+1] );
};
#include "Polynomial.inl"

51
src/meshlabplugins/filter_screened_poisson/Src/Polynomial.inl
View File

@ -306,6 +306,7 @@ int Polynomial<Degree>::getSolutions( double c , double* roots , double EPS ) co
for( int i=0 ; i<_rCount ; i++ ) if( fabs(_roots[i][1])<=EPS ) roots[rCount++] = _roots[i][0];
return rCount;
}
// The 0-th order B-spline
template< >
Polynomial< 0 > Polynomial< 0 >::BSplineComponent( int i )
{
@ -313,20 +314,56 @@ Polynomial< 0 > Polynomial< 0 >::BSplineComponent( int i )
p.coefficients[0] = 1.;
return p;
}
// The Degree-th order B-spline
template< int Degree >
Polynomial< Degree > Polynomial< Degree >::BSplineComponent( int i )
{
// B_d^i(x) = \int_x^1 B_{d-1}^{i}(y) dy + \int_0^x B_{d-1}^{i-1} y dy
// = \int_0^1 B_{d-1}^{i}(y) dy - \int_0^x B_{d-1}^{i}(y) dy + \int_0^x B_{d-1}^{i-1} y dy
Polynomial p;
if( i>0 )
{
Polynomial< Degree > _p = Polynomial< Degree-1 >::BSplineComponent( i-1 ).integral();
p -= _p;
p.coefficients[0] += _p(1);
}
if( i<Degree )
{
Polynomial< Degree > _p = Polynomial< Degree-1 >::BSplineComponent( i ).integral();
p -= _p;
p.coefficients[0] += _p(1);
}
if( i>0 )
{
Polynomial< Degree > _p = Polynomial< Degree-1 >::BSplineComponent( i-1 ).integral();
p += _p;
}
return p;
}
}
// The 0-th order B-spline values
template< > void Polynomial< 0 >::BSplineComponentValues( double x , double* values ){ values[0] = 1.; }
// The Degree-th order B-spline
template< int Degree > void Polynomial< Degree >::BSplineComponentValues( double x , double* values )
{
const double Scale = 1./Degree;
Polynomial< Degree-1 >::BSplineComponentValues( x , values+1 );
values[0] = values[1] * (1.-x) * Scale;
for( int i=1 ; i<Degree ; i++ )
{
double x1 = (x-i+Degree) , x2 = (-x+i+1);
values[i] = ( values[i]*x1 + values[i+1]*x2 ) * Scale;
}
values[Degree] *= x * Scale;
}
// Using the recurrence formulation for Pascal's triangle
template< > void Polynomial< 0 >::BinomialCoefficients( int bCoefficients[1] ){ bCoefficients[0] = 1; }
template< int Degree > void Polynomial< Degree >::BinomialCoefficients( int bCoefficients[Degree+1] )
{
Polynomial< Degree-1 >::BinomialCoefficients( bCoefficients );
int leftValue = 0;
for( int i=0 ; i<Degree ; i++ )
{
int temp = bCoefficients[i];
bCoefficients[i] += leftValue;
leftValue = temp;
}
bCoefficients[Degree] = 1;
}

760
src/meshlabplugins/filter_screened_poisson/Src/SSDRecon.cpp Normal file
View File

@ -0,0 +1,760 @@
/*
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.
*/
#undef FAST_COMPILE
#undef ARRAY_DEBUG
#define BRUNO_LEVY_FIX
#define FOR_RELEASE
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#if defined( _WIN32 ) || defined( _WIN64 )
#include <Windows.h>
#include <Psapi.h>
#endif // _WIN32 || _WIN64
#include "MyTime.h"
#include "MarchingCubes.h"
#include "Octree.h"
#include "SparseMatrix.h"
#include "CmdLineParser.h"
#include "PPolynomial.h"
#include "Ply.h"
#include "MemoryUsage.h"
#ifdef _OPENMP
#include "omp.h"
#endif // _OPENMP
void DumpOutput( const char* format , ... );
void DumpOutput2( std::vector< char* >& comments , const char* format , ... );
#include "MultiGridOctreeData.h"
#define DEFAULT_FULL_DEPTH 5
#define XSTR(x) STR(x)
#define STR(x) #x
#if DEFAULT_FULL_DEPTH
#pragma message ( "[WARNING] Setting default full depth to " XSTR(DEFAULT_FULL_DEPTH) )
#endif // DEFAULT_FULL_DEPTH
#include <stdarg.h>
char* outputFile=NULL;
int echoStdout=0;
void DumpOutput( const char* format , ... )
{
if( outputFile )
{
FILE* fp = fopen( outputFile , "a" );
va_list args;
va_start( args , format );
vfprintf( fp , format , args );
fclose( fp );
va_end( args );
}
if( echoStdout )
{
va_list args;
va_start( args , format );
vprintf( format , args );
va_end( args );
}
}
void DumpOutput2( std::vector< char* >& comments , const char* format , ... )
{
if( outputFile )
{
FILE* fp = fopen( outputFile , "a" );
va_list args;
va_start( args , format );
vfprintf( fp , format , args );
fclose( fp );
va_end( args );
}
if( echoStdout )
{
va_list args;
va_start( args , format );
vprintf( format , args );
va_end( args );
}
comments.push_back( new char[1024] );
char* str = comments.back();
va_list args;
va_start( args , format );
vsprintf( str , format , args );
va_end( args );
if( str[strlen(str)-1]=='\n' ) str[strlen(str)-1] = 0;
}
cmdLineString
In( "in" ) ,
Out( "out" ) ,
VoxelGrid( "voxel" ) ,
XForm( "xForm" );
cmdLineReadable
#if defined( _WIN32 ) || defined( _WIN64 )
Performance( "performance" ) ,
#endif // _WIN32 || _WIN64
ShowResidual( "showResidual" ) ,
NoComments( "noComments" ) ,
PolygonMesh( "polygonMesh" ) ,
Confidence( "confidence" ) ,
NormalWeights( "nWeights" ) ,
NonManifold( "nonManifold" ) ,
ASCII( "ascii" ) ,
Density( "density" ) ,
NonLinearFit( "nonLinearFit" ) ,
PrimalVoxel( "primalVoxel" ) ,
#ifndef FAST_COMPILE
FreeBoundary( "freeBoundary" ) ,
Double( "double" ) ,
#endif // !FAST_COMPILE
Verbose( "verbose" );
cmdLineInt
#ifndef FAST_COMPILE
Degree( "degree" , 2 ) ,
#endif // !FAST_COMPILE
Depth( "depth" , 8 ) ,
CGDepth( "cgDepth" , 0 ) ,
KernelDepth( "kernelDepth" ) ,
AdaptiveExponent( "adaptiveExp" , 1 ) ,
Iters( "iters" , 8 ) ,
VoxelDepth( "voxelDepth" , -1 ) ,
FullDepth( "fullDepth" , DEFAULT_FULL_DEPTH ) ,
MaxSolveDepth( "maxSolveDepth" ) ,
Threads( "threads" , omp_get_num_procs() );
cmdLineFloat
Color( "color" , 16.f ) ,
SamplesPerNode( "samplesPerNode" , 1.5f ) ,
Scale( "scale" , 1.1f ) ,
CGSolverAccuracy( "cgAccuracy" , 1e-3f ) ,
LowResIterMultiplier( "iterMultiplier" , 1.5f ) ,
ValueWeight ( "valueWeight" , 4e-0f ) ,
GradientWeight( "gradientWeight" , 1e-3f ) ,
BiLapWeight ( "biLapWeight" , 1e-5f );
cmdLineReadable* params[] =
{
#ifndef FAST_COMPILE
&Degree , &Double , &FreeBoundary ,
#endif // !FAST_COMPILE
&In , &Depth , &Out , &XForm ,
&Scale , &Verbose , &CGSolverAccuracy , &NoComments , &LowResIterMultiplier ,
&KernelDepth , &SamplesPerNode , &Confidence , &NormalWeights , &NonManifold , &PolygonMesh , &ASCII , &ShowResidual , &VoxelDepth ,
&BiLapWeight ,
&ValueWeight , &GradientWeight , &VoxelGrid , &Threads , &MaxSolveDepth ,
&AdaptiveExponent ,
&Density ,
&FullDepth ,
&CGDepth , &Iters ,
&Color ,
&NonLinearFit ,
&PrimalVoxel ,
#if defined( _WIN32 ) || defined( _WIN64 )
&Performance ,
#endif // _WIN32 || _WIN64
};
void ShowUsage( char* ex )
{
printf( "Usage: %s\n" , ex );
printf( "\t --%s <input points>\n" , In.name );
printf( "\t[--%s <ouput triangle mesh>]\n" , Out.name );
printf( "\t[--%s <ouput voxel grid>]\n" , VoxelGrid.name );
#ifndef FAST_COMPILE
printf( "\t[--%s <b-spline degree>=%d]\n" , Degree.name , Degree.value );
#ifndef FOR_RELEASE
printf( "\t[--%s]\n" , FreeBoundary.name );
#endif // !FOR_RELEASE
#endif // !FAST_COMPILE
printf( "\t[--%s <maximum reconstruction depth>=%d]\n" , Depth.name , Depth.value );
printf( "\t[--%s <scale factor>=%f]\n" , Scale.name , Scale.value );
printf( "\t[--%s <minimum number of samples per node>=%f]\n" , SamplesPerNode.name, SamplesPerNode.value );
printf( "\t[--%s <zero-crossing weight>=%.3e]\n" , ValueWeight.name , ValueWeight.value );
printf( "\t[--%s <gradient weight>=%.3e]\n" , GradientWeight.name , GradientWeight.value );
printf( "\t[--%s <bi-laplacian weight>=%.3e]\n" , BiLapWeight.name , BiLapWeight.value );
printf( "\t[--%s]\n" , Confidence.name );
printf( "\t[--%s]\n" , NormalWeights.name );
#ifndef FOR_RELEASE
printf( "\t[--%s <adaptive weighting exponent>=%d]\n", AdaptiveExponent.name , AdaptiveExponent.value );
#endif // !FOR_RELEASE
printf( "\t[--%s <iterations>=%d]\n" , Iters.name , Iters.value );
#ifndef FOR_RELEASE
printf( "\t[--%s <low-resolution iteration multiplier>=%f]\n" , LowResIterMultiplier.name , LowResIterMultiplier.value );
#endif // FOR_RELEASE
printf( "\t[--%s <conjugate-gradients depth>=%d]\n" , CGDepth.name , CGDepth.value );
#ifndef FOR_RELEASE
printf( "\t[--%s <conjugate-gradients solver accuracy>=%g]\n" , CGSolverAccuracy.name , CGSolverAccuracy.value );
#endif // !FOR_RELEASE
printf( "\t[--%s <full depth>=%d]\n" , FullDepth.name , FullDepth.value );
printf( "\t[--%s <depth at which to extract the voxel grid>=<%s>]\n" , VoxelDepth.name , Depth.name );
printf( "\t[--%s]\n" , PrimalVoxel.name );
printf( "\t[--%s <pull factor>]\n" , Color.name );
printf( "\t[--%s]\n" , Density.name );
printf( "\t[--%s]\n" , NonLinearFit.name );
printf( "\t[--%s]\n" , PolygonMesh.name);
#ifndef FOR_RELEASE
printf( "\t[--%s]\n" , NonManifold.name );
#endif // !FOR_RELEASE
#ifdef _OPENMP
printf( "\t[--%s <num threads>=%d]\n" , Threads.name , Threads.value );
#endif // _OPENMP
printf( "\t[--%s]\n" , Verbose.name );
#ifndef FOR_RELEASE
#if defined( _WIN32 ) || defined( _WIN64 )
printf( "\t[--%s]\n" , Performance.name );
#endif // _WIN32 || _WIN64
#endif // !FOR_RELEASE
#ifndef FOR_RELEASE
printf( "\t[--%s]\n" , ASCII.name );
printf( "\t[--%s]\n" , NoComments.name );
#endif // !FOR_RELEASE
#ifndef FAST_COMPILE
printf( "\t[--%s]\n" , Double.name );
#endif // !FAST_COMPILE
}
template< class Real >
struct ColorInfo
{
static Point3D< Real > ReadASCII( FILE* fp )
{
Point3D< unsigned char > c;
if( fscanf( fp , " %c %c %c " , &c[0] , &c[1] , &c[2] )!=3 ) fprintf( stderr , "[ERROR] Failed to read color\n" ) , exit( 0 );
return Point3D< Real >( (Real)c[0] , (Real)c[1] , (Real)c[2] );
};
static bool ValidPlyProperties( const bool* props ){ return ( props[0] || props[3] ) && ( props[1] || props[4] ) && ( props[2] || props[5] ); }
const static PlyProperty PlyProperties[];
};
template<>
const PlyProperty ColorInfo< float >::PlyProperties[] =
{
{ "r" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[0] ) ) , 0 , 0 , 0 , 0 } ,
{ "g" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[1] ) ) , 0 , 0 , 0 , 0 } ,
{ "b" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[2] ) ) , 0 , 0 , 0 , 0 } ,
{ "red" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[0] ) ) , 0 , 0 , 0 , 0 } ,
{ "green" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[1] ) ) , 0 , 0 , 0 , 0 } ,
{ "blue" , PLY_UCHAR , PLY_FLOAT , int( offsetof( Point3D< float > , coords[2] ) ) , 0 , 0 , 0 , 0 }
};
template<>
const PlyProperty ColorInfo< double >::PlyProperties[] =
{
{ "r" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[0] ) ) , 0 , 0 , 0 , 0 } ,
{ "g" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[1] ) ) , 0 , 0 , 0 , 0 } ,
{ "b" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[2] ) ) , 0 , 0 , 0 , 0 } ,
{ "red" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[0] ) ) , 0 , 0 , 0 , 0 } ,
{ "green" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[1] ) ) , 0 , 0 , 0 , 0 } ,
{ "blue" , PLY_UCHAR , PLY_DOUBLE , int( offsetof( Point3D< double > , coords[2] ) ) , 0 , 0 , 0 , 0 }
};
bool ValidPlyColorProperties( const bool* props ){ return ( props[0] || props[3] ) && ( props[1] || props[4] ) && ( props[2] || props[5] ); }
double Weight( double v , double start , double end )
{
v = ( v - start ) / ( end - start );
if ( v<0 ) return 1.;
else if( v>1 ) return 0.;
else
{
// P(x) = a x^3 + b x^2 + c x + d
// P (0) = 1 , P (1) = 0 , P'(0) = 0 , P'(1) = 0
// => d = 1 , a + b + c + d = 0 , c = 0 , 3a + 2b + c = 0
// => c = 0 , d = 1 , a + b = -1 , 3a + 2b = 0
// => a = 2 , b = -3 , c = 0 , d = 1
// => P(x) = 2 x^3 - 3 x^2 + 1
return 2. * v * v * v - 3. * v * v + 1.;
}
}
#if defined( _WIN32 ) || defined( _WIN64 )
double PeakMemoryUsageMB( void )
{
HANDLE h = GetCurrentProcess();
PROCESS_MEMORY_COUNTERS pmc;
return GetProcessMemoryInfo( h , &pmc , sizeof(pmc) ) ? ( (double)pmc.PeakWorkingSetSize )/(1<<20) : 0;
}
#endif // _WIN32 || _WIN64
template< class Real >
struct OctreeProfiler
{
Octree< Real >& tree;
double t;
OctreeProfiler( Octree< Real >& t ) : tree(t) { ; }
void start( void ){ t = Time() , tree.resetLocalMemoryUsage(); }
void print( const char* header ) const
{
tree.memoryUsage();
#if defined( _WIN32 ) || defined( _WIN64 )
if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
else printf( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
#else // !_WIN32 && !_WIN64
if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
else printf( "%9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
#endif // _WIN32 || _WIN64
}
void dumpOutput( const char* header ) const
{
tree.memoryUsage();
#if defined( _WIN32 ) || defined( _WIN64 )
if( header ) DumpOutput( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
else DumpOutput( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
#else // !_WIN32 && !_WIN64
if( header ) DumpOutput( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
else DumpOutput( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
#endif // _WIN32 || _WIN64
}
void dumpOutput2( std::vector< char* >& comments , const char* header ) const
{
tree.memoryUsage();
#if defined( _WIN32 ) || defined( _WIN64 )
if( header ) DumpOutput2( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
else DumpOutput2( comments , "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() , PeakMemoryUsageMB() );
#else // !_WIN32 && !_WIN64
if( header ) DumpOutput2( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
else DumpOutput2( comments , "%9.1f (s), %9.1f (MB) / %9.1f (MB)\n" , Time()-t , tree.localMemoryUsage() , tree.maxMemoryUsage() );
#endif // _WIN32 || _WIN64
}
};
template< class Real >
XForm4x4< Real > GetPointXForm( OrientedPointStream< Real >& stream , Real scaleFactor )
{
Point3D< Real > min , max;
stream.boundingBox( min , max );
Point3D< Real > center = ( max + min ) / 2;
Real scale = std::max< Real >( max[0]-min[0] , std::max< Real >( max[1]-min[1] , max[2]-min[2] ) );
scale *= scaleFactor;
for( int i=0 ; i<3 ; i++ ) center[i] -= scale/2;
XForm4x4< Real > tXForm = XForm4x4< Real >::Identity() , sXForm = XForm4x4< Real >::Identity();
for( int i=0 ; i<3 ; i++ ) sXForm(i,i) = (Real)(1./scale ) , tXForm(3,i) = -center[i];
return sXForm * tXForm;
}
template< class Real , int Degree , BoundaryType BType , class Vertex >
int _Execute( int argc , char* argv[] )
{
typedef typename Octree< Real >::template InterpolationInfo< true > InterpolationInfo;
typedef OrientedPointStream< Real > PointStream;
typedef OrientedPointStreamWithData< Real , Point3D< Real > > PointStreamWithData;
typedef TransformedOrientedPointStream< Real > XPointStream;
typedef TransformedOrientedPointStreamWithData< Real , Point3D< Real > > XPointStreamWithData;
Reset< Real >();
int paramNum = sizeof(params)/sizeof(cmdLineReadable*);
std::vector< char* > comments;
if( Verbose.set ) echoStdout=1;
XForm4x4< Real > xForm , iXForm;
if( XForm.set )
{
FILE* fp = fopen( XForm.value , "r" );
if( !fp )
{
fprintf( stderr , "[WARNING] Could not read x-form from: %s\n" , XForm.value );
xForm = XForm4x4< Real >::Identity();
}
else
{
for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ )
{
float f;
if( fscanf( fp , " %f " , &f )!=1 ) fprintf( stderr , "[ERROR] Execute: Failed to read xform\n" ) , exit( 0 );
xForm(i,j) = (Real)f;
}
fclose( fp );
}
}
else xForm = XForm4x4< Real >::Identity();
DumpOutput2( comments , "Running SSD Reconstruction (Version 9.0)\n" );
char str[1024];
for( int i=0 ; i<paramNum ; i++ )
if( params[i]->set )
{
params[i]->writeValue( str );
if( strlen( str ) ) DumpOutput2( comments , "\t--%s %s\n" , params[i]->name , str );
else DumpOutput2( comments , "\t--%s\n" , params[i]->name );
}
double startTime = Time();
Real isoValue = 0;
Octree< Real > tree;
OctreeProfiler< Real > profiler( tree );
tree.threads = Threads.value;
if( !In.set )
{
ShowUsage( argv[0] );
return 0;
}
if( !MaxSolveDepth.set ) MaxSolveDepth.value = Depth.value;
OctNode< TreeNodeData >::SetAllocator( MEMORY_ALLOCATOR_BLOCK_SIZE );
int kernelDepth = KernelDepth.set ? KernelDepth.value : Depth.value-2;
if( kernelDepth>Depth.value )
{
fprintf( stderr,"[WARNING] %s can't be greater than %s: %d <= %d\n" , KernelDepth.name , Depth.name , KernelDepth.value , Depth.value );
kernelDepth = Depth.value;
}
int pointCount;
Real pointWeightSum;
std::vector< typename Octree< Real >::PointSample >* samples = new std::vector< typename Octree< Real >::PointSample >();
std::vector< ProjectiveData< Point3D< Real > , Real > >* sampleData = NULL;
SparseNodeData< Real , WEIGHT_DEGREE >* density = NULL;
SparseNodeData< Point3D< Real > , NORMAL_DEGREE >* normalInfo = NULL;
Real targetValue = (Real)0.;
// Read in the samples (and color data)
{
profiler.start();
PointStream* pointStream;
char* ext = GetFileExtension( In.value );
if( Color.set && Color.value>0 )
{
sampleData = new std::vector< ProjectiveData< Point3D< Real > , Real > >();
if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStreamWithData< Real , Point3D< Real > , float , Point3D< unsigned char > >( In.value );
else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYOrientedPointStreamWithData< Real , Point3D< Real > >( In.value , ColorInfo< Real >::PlyProperties , 6 , ColorInfo< Real >::ValidPlyProperties );
else pointStream = new ASCIIOrientedPointStreamWithData< Real , Point3D< Real > >( In.value , ColorInfo< Real >::ReadASCII );
}
else
{
if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryOrientedPointStream< Real , float >( In.value );
else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYOrientedPointStream< Real >( In.value );
else pointStream = new ASCIIOrientedPointStream< Real >( In.value );
}
delete[] ext;
XPointStream _pointStream( xForm , *pointStream );
xForm = GetPointXForm( _pointStream , (Real)Scale.value ) * xForm;
if( sampleData )
{
XPointStreamWithData _pointStream( xForm , ( PointStreamWithData& )*pointStream );
pointCount = tree.template init< Point3D< Real > >( _pointStream , Depth.value , Confidence.set , *samples , sampleData );
}
else
{
XPointStream _pointStream( xForm , *pointStream );
pointCount = tree.template init< Point3D< Real > >( _pointStream , Depth.value , Confidence.set , *samples , sampleData );
}
iXForm = xForm.inverse();
delete pointStream;
#pragma omp parallel for num_threads( Threads.value )
for( int i=0 ; i<(int)samples->size() ; i++ ) (*samples)[i].sample.data.n *= (Real)-1;
DumpOutput( "Input Points / Samples: %d / %d\n" , pointCount , samples->size() );
profiler.dumpOutput2( comments , "# Read input into tree:" );
}
DenseNodeData< Real , Degree > solution;
// Solve
{
DenseNodeData< Real , Degree > constraints;
InterpolationInfo* iInfo = NULL;
int solveDepth = MaxSolveDepth.value;
tree.resetNodeIndices();
// Get the kernel density estimator
{
profiler.start();
density = new SparseNodeData< Real , WEIGHT_DEGREE >();
*density = tree.template setDensityEstimator< WEIGHT_DEGREE >( *samples , kernelDepth , SamplesPerNode.value );
profiler.dumpOutput2( comments , "# Got kernel density:" );
}
// Transform the Hermite samples into a vector field
{
profiler.start();
normalInfo = new SparseNodeData< Point3D< Real > , NORMAL_DEGREE >();
*normalInfo = tree.template setNormalField< NORMAL_DEGREE >( *samples , *density , pointWeightSum , BType==BOUNDARY_NEUMANN );
profiler.dumpOutput2( comments , "# Got normal field:" );
}
if( !Density.set ) delete density , density = NULL;
// Trim the tree and prepare for multigrid
{
profiler.start();
std::vector< int > indexMap;
constexpr int MAX_DEGREE = NORMAL_DEGREE > Degree ? NORMAL_DEGREE : Degree;
tree.template inalizeForBroodedMultigrid< MAX_DEGREE , Degree , BType >( FullDepth.value , typename Octree< Real >::template HasNormalDataFunctor< NORMAL_DEGREE >( *normalInfo ) , &indexMap );
if( normalInfo ) normalInfo->remapIndices( indexMap );
if( density ) density->remapIndices( indexMap );
profiler.dumpOutput2( comments , "# Finalized tree:" );
}
// Free up the normal info
if( normalInfo ) delete normalInfo , normalInfo = NULL;
// Add the interpolation constraints
if( ValueWeight.value>0 || GradientWeight.value>0 )
{
profiler.start();
iInfo = new InterpolationInfo( tree , *samples , targetValue , AdaptiveExponent.value , (Real)ValueWeight.value * pointWeightSum , (Real)GradientWeight.value * pointWeightSum );
constraints = tree.template initDenseNodeData< Degree >( );
tree.template addInterpolationConstraints< Degree , BType >( *iInfo , constraints , solveDepth );
profiler.dumpOutput2( comments , "#Set point constraints:" );
}
DumpOutput( "Leaf Nodes / Active Nodes / Ghost Nodes: %d / %d / %d\n" , (int)tree.leaves() , (int)tree.nodes() , (int)tree.ghostNodes() );
DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage())/(1<<20) );
// Solve the linear system
{
profiler.start();
typename Octree< Real >::SolverInfo solverInfo;
solverInfo.cgDepth = CGDepth.value , solverInfo.iters = Iters.value , solverInfo.cgAccuracy = CGSolverAccuracy.value , solverInfo.verbose = Verbose.set , solverInfo.showResidual = ShowResidual.set , solverInfo.lowResIterMultiplier = std::max< double >( 1. , LowResIterMultiplier.value );
solution = tree.template solveSystem< Degree , BType >( FEMSystemFunctor< Degree , BType >( 0 , 0 , BiLapWeight.value ) , iInfo , constraints , solveDepth , solverInfo );
profiler.dumpOutput2( comments , "# Linear system solved:" );
DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) );
if( iInfo ) delete iInfo , iInfo = NULL;
}
}
CoredFileMeshData< Vertex > mesh;
{
profiler.start();
double valueSum = 0 , weightSum = 0;
typename Octree< Real >::template MultiThreadedEvaluator< Degree , BType > evaluator( &tree , solution , Threads.value );
#pragma omp parallel for num_threads( Threads.value ) reduction( + : valueSum , weightSum )
for( int j=0 ; j<samples->size() ; j++ )
{
ProjectiveData< OrientedPoint3D< Real > , Real >& sample = (*samples)[j].sample;
Real w = sample.weight;
if( w>0 ) weightSum += w , valueSum += evaluator.value( sample.data.p / sample.weight , omp_get_thread_num() , (*samples)[j].node ) * w;
}
isoValue = (Real)( valueSum / weightSum );
if( !( Color.set && Color.value>0 ) && samples ) delete samples , samples = NULL;
profiler.dumpOutput( "Got average:" );
DumpOutput( "Iso-Value: %e\n" , isoValue );
}
if( VoxelGrid.set )
{
profiler.start();
FILE* fp = fopen( VoxelGrid.value , "wb" );
if( !fp ) fprintf( stderr , "Failed to open voxel file for writing: %s\n" , VoxelGrid.value );
else
{
int res = 0;
Pointer( Real ) values = tree.template voxelEvaluate< Real , Degree , BType >( solution , res , isoValue , VoxelDepth.value , PrimalVoxel.set );
fwrite( &res , sizeof(int) , 1 , fp );
if( sizeof(Real)==sizeof(float) ) fwrite( values , sizeof(float) , res*res*res , fp );
else
{
float *fValues = new float[res*res*res];
for( int i=0 ; i<res*res*res ; i++ ) fValues[i] = float( values[i] );
fwrite( fValues , sizeof(float) , res*res*res , fp );
delete[] fValues;
}
fclose( fp );
DeletePointer( values );
}
profiler.dumpOutput( "Got voxel grid:" );
}
if( Out.set )
{
profiler.start();
SparseNodeData< ProjectiveData< Point3D< Real > , Real > , DATA_DEGREE >* colorData = NULL;
if( sampleData )
{
colorData = new SparseNodeData< ProjectiveData< Point3D< Real > , Real > , DATA_DEGREE >();
*colorData = tree.template setDataField< DATA_DEGREE , false >( *samples , *sampleData , (SparseNodeData< Real , WEIGHT_DEGREE >*)NULL );
delete sampleData , sampleData = NULL;
for( const OctNode< TreeNodeData >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) )
{
ProjectiveData< Point3D< Real > , Real >* clr = (*colorData)( n );
if( clr ) (*clr) *= (Real)pow( Color.value , tree.depth( n ) );
}
}
tree.template getMCIsoSurface< Degree , BType , WEIGHT_DEGREE , DATA_DEGREE >( density , colorData , solution , isoValue , mesh , NonLinearFit.set , !NonManifold.set , PolygonMesh.set );
DumpOutput( "Vertices / Polygons: %d / %d\n" , mesh.outOfCorePointCount()+mesh.inCorePoints.size() , mesh.polygonCount() );
if( PolygonMesh.set ) profiler.dumpOutput2( comments , "# Got polygons:" );
else profiler.dumpOutput2( comments , "# Got triangles:" );
if( colorData ) delete colorData , colorData = NULL;
if( NoComments.set )
{
if( ASCII.set ) PlyWritePolygons( Out.value , &mesh , PLY_ASCII , NULL , 0 , iXForm );
else PlyWritePolygons( Out.value , &mesh , PLY_BINARY_NATIVE , NULL , 0 , iXForm );
}
else
{
if( ASCII.set ) PlyWritePolygons( Out.value , &mesh , PLY_ASCII , &comments[0] , (int)comments.size() , iXForm );
else PlyWritePolygons( Out.value , &mesh , PLY_BINARY_NATIVE , &comments[0] , (int)comments.size() , iXForm );
}
}
if( density ) delete density , density = NULL;
DumpOutput2( comments , "# Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-startTime , tree.maxMemoryUsage() );
return 1;
}
#if defined( _WIN32 ) || defined( _WIN64 )
inline double to_seconds( const FILETIME& ft )
{
const double low_to_sec=100e-9; // 100 nanoseconds
const double high_to_sec=low_to_sec*4294967296.0;
return ft.dwLowDateTime*low_to_sec+ft.dwHighDateTime*high_to_sec;
}
#endif // _WIN32 || _WIN64
#ifndef FAST_COMPILE
template< class Real , class Vertex >
int Execute( int argc , char* argv[] )
{
if( FreeBoundary.set )
switch( Degree.value )
{
case 2: return _Execute< Real , 2 , BOUNDARY_FREE , Vertex >( argc , argv );
case 3: return _Execute< Real , 3 , BOUNDARY_FREE , Vertex >( argc , argv );
case 4: return _Execute< Real , 4 , BOUNDARY_FREE , Vertex >( argc , argv );
default: fprintf( stderr , "[ERROR] Only B-Splines of degree 2 - 4 are supported" ) ; return EXIT_FAILURE;
}
else
switch( Degree.value )
{
case 2: return _Execute< Real , 2 , BOUNDARY_NEUMANN , Vertex >( argc , argv );
case 3: return _Execute< Real , 3 , BOUNDARY_NEUMANN , Vertex >( argc , argv );
case 4: return _Execute< Real , 4 , BOUNDARY_NEUMANN , Vertex >( argc , argv );
default: fprintf( stderr , "[ERROR] Only B-Splines of degree 2 - 4 are supported" ) ; return EXIT_FAILURE;
}
}
#endif // !FAST_COMPILE
int main( int argc , char* argv[] )
{
#if defined(WIN32) && defined(MAX_MEMORY_GB)
if( MAX_MEMORY_GB>0 )
{
SIZE_T peakMemory = 1;
peakMemory <<= 30;
peakMemory *= MAX_MEMORY_GB;
printf( "Limiting memory usage to %.2f GB\n" , float( peakMemory>>30 ) );
HANDLE h = CreateJobObject( NULL , NULL );
AssignProcessToJobObject( h , GetCurrentProcess() );
JOBOBJECT_EXTENDED_LIMIT_INFORMATION jeli = { 0 };
jeli.BasicLimitInformation.LimitFlags = JOB_OBJECT_LIMIT_JOB_MEMORY;
jeli.JobMemoryLimit = peakMemory;
if( !SetInformationJobObject( h , JobObjectExtendedLimitInformation , &jeli , sizeof( jeli ) ) )
fprintf( stderr , "Failed to set memory limit\n" );
}
#endif // defined(WIN32) && defined(MAX_MEMORY_GB)
double t = Time();
cmdLineParse( argc-1 , &argv[1] , sizeof(params)/sizeof(cmdLineReadable*) , params , 1 );
if( GradientWeight.value<=0 ) fprintf( stderr , "[ERROR] Gradient weight must be positive: %g>=0\n" , GradientWeight.value ) , exit( 0 );
if( BiLapWeight.value<=0 ) fprintf( stderr , "[ERROR] Bi-Laplacian weight must be positive: %g>=0\n" , BiLapWeight.value ) , exit( 0 );
#ifdef FAST_COMPILE
static const int Degree = 2;
static const BoundaryType BType = BOUNDARY_NEUMANN;
fprintf( stderr , "[WARNING] Compiling for degree-%d, boundary-%s, single-precision _only_\n" , Degree , BoundaryNames[ BType ] );
if( Density.set )
if( Color.set && Color.value>0 ) _Execute< float , Degree , BType , PlyColorAndValueVertex< float > >( argc , argv );
else _Execute< float , Degree , BType , PlyValueVertex< float > >( argc , argv );
else
if( Color.set && Color.value>0 ) _Execute< float , Degree , BType , PlyColorVertex< float > >( argc , argv );
else _Execute< float , Degree , BType , PlyVertex< float > >( argc , argv );
#else // !FAST_COMPILE
if( Density.set )
if( Color.set && Color.value>0 )
if( Double.set ) Execute< double , PlyColorAndValueVertex< float > >( argc , argv );
else Execute< float , PlyColorAndValueVertex< float > >( argc , argv );
else
if( Double.set ) Execute< double , PlyValueVertex< float > >( argc , argv );
else Execute< float , PlyValueVertex< float > >( argc , argv );
else
if( Color.set && Color.value>0 )
if( Double.set ) Execute< double , PlyColorVertex< float > >( argc , argv );
else Execute< float , PlyColorVertex< float > >( argc , argv );
else
if( Double.set ) Execute< double , PlyVertex< float > >( argc , argv );
else Execute< float , PlyVertex< float > >( argc , argv );
#endif // FAST_COMPILE
#if defined( _WIN32 ) || defined( _WIN64 )
if( Performance.set )
{
HANDLE cur_thread=GetCurrentThread();
FILETIME tcreat, texit, tkernel, tuser;
if( GetThreadTimes( cur_thread , &tcreat , &texit , &tkernel , &tuser ) )
printf( "Time (Wall/User/Kernel): %.2f / %.2f / %.2f\n" , Time()-t , to_seconds( tuser ) , to_seconds( tkernel ) );
else printf( "Time: %.2f\n" , Time()-t );
HANDLE h = GetCurrentProcess();
PROCESS_MEMORY_COUNTERS pmc;
if( GetProcessMemoryInfo( h , &pmc , sizeof(pmc) ) ) printf( "Peak Memory (MB): %d\n" , (int)(pmc.PeakWorkingSetSize>>20) );
}
#endif // _WIN32 || _WIN64
return EXIT_SUCCESS;
}

42
src/meshlabplugins/filter_screened_poisson/Src/SparseMatrix.h
View File

@ -29,10 +29,10 @@ DAMAGE.
#ifndef __SPARSEMATRIX_HPP
#define __SPARSEMATRIX_HPP
#define NEW_SPARSE_MATRIX 1
#define ZERO_TESTING_JACOBI 1
#include "Vector.h"
#include "Array.h"
template <class T>
@ -76,43 +76,26 @@ public:
SparseMatrix<T> operator * (const T& V) const;
SparseMatrix<T>& operator *= (const T& V);
template<class T2>
Vector<T2> operator * (const Vector<T2>& V) const;
template<class T2>
Vector<T2> Multiply( const Vector<T2>& V ) const;
template<class T2>
void Multiply( const Vector<T2>& In , Vector<T2>& Out , int threads=1 ) const;
static int Solve (const SparseMatrix<T>& M,const Vector<T>& b, int iters,Vector<T>& solution,const T eps=1e-8);
template<class T2>
static int SolveSymmetric( const SparseMatrix<T>& M , const Vector<T2>& b , int iters , Vector<T2>& solution , const T2 eps=1e-8 , int reset=1 , int threads=1 );
template< class T2 > void Multiply( ConstPointer( T2 ) in , Pointer( T2 ) out , int threads=1 ) const;
template< class T2 > void MultiplyAndAddAverage( ConstPointer( T2 ) in , Pointer( T2 ) out , int threads=1 ) const;
bool write( FILE* fp ) const;
bool write( const char* fileName ) const;
bool read( FILE* fp );
bool read( const char* fileName );
template< class T2 >
static int SolveJacobi( const SparseMatrix<T>& M , const Vector<T2>& diagonal , const Vector<T2>& b , Vector<T2>& x , Vector<T2>& Mx , T2 sor , int threads=1 , int offset=0 );
template< class T2 >
static int SolveGS( const SparseMatrix<T>& M , const Vector<T2>& diagonal , const Vector<T2>& b , Vector<T2>& x , bool forward , int offset=0 );
template< class T2 >
static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix<T>& M , const Vector<T2>& diagonal , const Vector<T2>& b , Vector<T2>& x , bool forward , int threads=1 , int offset=0 );
template< class T2 >
static int SolveJacobi( const SparseMatrix<T>& M , const Vector<T2>& b , Vector<T2>& x , Vector<T2>& Mx , T2 sor , int threads=1 , int offset=0 );
template< class T2 >
static int SolveGS( const SparseMatrix<T>& M , const Vector<T2>& b , Vector<T2>& x , bool forward , int offset=0 );
template< class T2 >
static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix<T>& M , const Vector<T2>& b , Vector<T2>& x , bool forward , int threads=1 , int offset=0 );
template< class T2 >
void getDiagonal( Vector< T2 >& diagonal , int threads=1 , int offset=0 ) const;
template< class T2 > void getDiagonal( Pointer( T2 ) diagonal , int threads=1 ) const;
template< class T2 > static int SolveJacobi( const SparseMatrix<T>& M , ConstPointer( T2 ) b , Pointer( T2 ) x , Pointer( T2 ) Mx , T2 sor , int threads=1 );
template< class T2 > static int SolveJacobi( const SparseMatrix<T>& M , ConstPointer( T2 ) diagonal , ConstPointer( T2 ) b , Pointer( T2 ) x , Pointer( T2 ) Mx , T2 sor , int threads=1 );
template< class T2 > static int SolveGS( const SparseMatrix<T>& M , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward );
template< class T2 > static int SolveGS( const SparseMatrix<T>& M , ConstPointer( T2 ) diagonal , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward );
template< class T2 > static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix<T>& M , ConstPointer( T2 ) diagonal , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward , int threads=1 );
template< class T2 > static int SolveGS( const std::vector< std::vector< int > >& mcIndices , const SparseMatrix<T>& M , ConstPointer( T2 ) b , Pointer( T2 ) x , bool forward , int threads=1 );
template< class T2 > static int SolveCG( const SparseMatrix<T>& M , ConstPointer( T2 ) b , int iters , Pointer( T2 ) x , T2 eps=1e-8 , int reset=1 , bool addDCTerm=false , bool solveNormal=false , int threads=1 );
};
#if !NEW_SPARSE_MATRIX
template< class T2 >
struct MapReduceVector
{
@ -203,6 +186,7 @@ public:
template< class T2 >
void getDiagonal( Vector< T2 >& diagonal , int threads=1 ) const;
};
#endif // !NEW_SPARSE_MATRIX
#include "SparseMatrix.inl"