+
+
+PointInterpolant:
+Fits a function to a set of sample values (and gradients)
+
+
+
--in <input sample positions and values>
+
This string is the name of the file from which the positions and values will be read.
+The file should be an ascii file with groups of Dim+1 or 2*Dim+1 (depending on whether gradients are provided or note)
+white space delimited, numbers: the coordinates of the point's position, followed by the value at that point (and gradient).
+No information about the number of samples should be specified.
+
+
[--dim <dimension of the samples>]
+
This integerl value is the dimension of the samples.
+The default value is 2.
+
+
[--tree <output octree and coefficients>]
+
This string is the name of the file to which the the octree and function coefficients are to be written.
+
+
[--grid <output grid>]
+
This string is the name of the file to which the sampled implicit function will be written.
+The file is wrtten out in binary, with the first 4 bytes corresponding to the (integer) sampling resolution, 2^d,
+and the next 4 x 2^d x 2^d x ... bytes corresponding to the (single precision) floating point values
+of the implicit function.
+
+
[--degree <B-spline degree>]
+
This integer specifies the degree of the B-spline that is to be used to define the finite elements system.
+Larger degrees support higher order approximations, but come at the cost of denser system matrices (incurring a cost in both space and time).
+The default value for this parameter is 2.
+
+
[--bType <boundary type>]
+
This integer specifies the boundary type for the finite elements. Valid values are:
+
+
1: Free boundary constraints
+
2: Dirichlet boundary constraints
+
3: Neumann boundary constraints
+
+The default value for this parameter is 1 (free).
+
+
[--depth <reconstruction depth>]
+
This integer is the maximum depth of the tree that will be used for surface reconstruction.
+Running at depth d corresponds to solving on a grid whose resolution is no larger than
+2^d x 2^d x ... Note that since the reconstructor adapts the octree to the
+sampling density, the specified reconstruction depth is only an upper bound.
+The default value for this parameter is 8.
+
+
[--width <finest cell width>]
+
This floating point value specifies the target width of the finest level octree cells.
+This parameter is ignored if the --depth is also specified.
+
+
[--scale <scale factor>]
+
This floating point value specifies the ratio between the diameter of the cube used for reconstruction
+and the diameter of the samples' bounding cube.
+The default value is 1.1.
+
+
[--valueWeight <value interpolation weight>]
+
This floating point value specifies the importance that interpolation of the samples' values
+is given in the fitting of the function.
+The default value for this parameter is 1000.
This floating point value specifies the importance that interpolation of the samples' gradients
+is given in the fitting of the function.
+The default value for this parameter is 1.
+This value is ignored unless gradient interpolation is specified.
+
+
[--lapWeight <Laplacian weight>]
+
This floating point value specifies the importance that Laplacian regularization
+is given in the fitting of the function.
+The default value for this parameter is 0.
+
+
[--biLapWeight <bi-Laplacian weight>]
+
This floating point value specifies the importance that bi-Laplacian regularization
+is given in the fitting of the function.
+The default value for this parameter is 1.
+
+
[--iters <GS iters>]
+
This integer value specifies the number of Gauss-Seidel relaxations to be performed at each level of the hiearchy.
+The default value for this parameter is 8.
+
+
[--useGradients]
+
Enabling this flag indicates that the input file contains gradients as well as sample values.
+
+
[--performance]
+
Enabling this flag provides running time and peak memory usage at the end of the execution.
+
+
[--verbose]
+
Enabling this flag provides a more verbose description of the running times and memory usages of
+individual components of the surface reconstructor.
+
+
+
+
+
+
@@ -504,9 +602,15 @@ The default value for this parameter is equal to the numer of (virtual) processo
Extracts iso-surfaces and a sampling on a regular grid from an implicit function represented over an adapted tree
--in <input tree and coefficients>
-
This string is the name of the file from which the tree and implicit functions coefficients are to be read.
+
This string is the name of the file from which the tree and implicit functions coefficients are to be read.
-
[--grid <output value grid>]
+
[--samples <input sample positions>]
+
This string is the name of the file from which sampling positions are to be read.
+The file should be an ascii file with groups of Dim white space delimited, numbers giving the coordinates of the sampling points' position.
+No information about the number of samples should be specified.
+
+
+
[--grid <output value grid>]
This string is the name of the file to which the sampling of the implicit along a regular grid will be written.
The file is written out in binary, with the first 4 bytes corresponding to the (integer) sampling resolution, R,
and the next 4 x R^D bytes corresponding to the (single precision) floating point values of the implicit function. (Here, D is the dimension.)
@@ -629,6 +733,36 @@ which partitions the reconstruction into 11 pieces.
+
+
+
+
+PointInterpolant / AdaptiveTreeVisualization
+
+For testing purposes, a pair of point-sets is provided:
+
+
+
fitting samples:
+A set of 1000 random 2D samples from within the square [-1,1,]x[-1,1] along with the evaluation of the quadratic f(x,y)=x*x+y*y at each sample point (represented in ASCII format).
+
evaluation samples:
+A set of 4 2D positions at which the fit function is to be evaluated (represented in ASCII format).
+
+
+The function fitting the input samples can be by calling the point interpolant:
+
+Note that because the last evaluation position, (2,2), is outside the bounding box of the fitting samples, the function cannot be evaluated at this point and a value of "nan" is output.
+
+
+
@@ -871,4 +1005,3 @@ Similarly, to reduce compilation times, support for specific degrees can be remo
SUPPORT
This work genersouly supported by NSF grants #0746039 and #1422325.
-
diff --git a/Src/AdaptiveTreeVisualization.cpp b/Src/AdaptiveTreeVisualization.cpp
index eabcafa..ac6a32d 100644
--- a/Src/AdaptiveTreeVisualization.cpp
+++ b/Src/AdaptiveTreeVisualization.cpp
@@ -42,6 +42,7 @@ DAMAGE.
cmdLineParameter< char* >
In( "in" ) ,
+ Samples( "samples" ) ,
OutMesh( "mesh" ) ,
OutGrid( "grid" );
@@ -66,6 +67,7 @@ cmdLineParameter< float >
cmdLineReadable* params[] =
{
&In ,
+ &Samples ,
&OutMesh , &NonManifold , &PolygonMesh , &FlipOrientation , &ASCII , &NonLinearFit , &IsoValue ,
&OutGrid , &PrimalGrid ,
&Threads ,
@@ -81,6 +83,7 @@ void ShowUsage( char* ex )
{
printf( "Usage: %s\n" , ex );
printf( "\t --%s \n" , In.name );
+ printf( "\t[--%s sample positions>]\n" , Samples.name );
printf( "\t[--%s ]\n" , OutMesh.name );
printf( "\t[--%s ]\n" , OutGrid.name );
printf( "\t[--%s =%d]\n" , Threads.name , Threads.value );
@@ -158,7 +161,7 @@ void WriteGrid( ConstPointer( Real ) values , int res , const char *fileName )
}
template< unsigned int Dim , class Real , unsigned int FEMSig >
-void _Execute( const FEMTree< Dim , Real >* tree , FILE* fp )
+void _Execute( const FEMTree< Dim , Real >* tree , XForm< Real , Dim+1 > xForm , FILE* fp )
{
ThreadPool::Init( (ThreadPool::ParallelType)ParallelType.value , Threads.value );
static const unsigned int Degree = FEMSignature< FEMSig >::Degree;
@@ -166,6 +169,37 @@ void _Execute( const FEMTree< Dim , Real >* tree , FILE* fp )
coefficients.read( fp );
+ // Evaluate at the sample positions
+ if( Samples.set )
+ {
+ InputPointStream< Real , Dim > *pointStream;
+ char* ext = GetFileExtension( Samples.value );
+ if ( !strcasecmp( ext , "bpts" ) ) pointStream = new BinaryInputPointStream< Real , Dim >( Samples.value );
+ else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYInputPointStream< Real , Dim >( Samples.value );
+ else pointStream = new ASCIIInputPointStream< Real , Dim >( Samples.value );
+ delete[] ext;
+ typename FEMTree< Dim , Real >::template MultiThreadedEvaluator< IsotropicUIntPack< Dim , FEMSig > , 0 > evaluator( tree , coefficients );
+ static const unsigned int CHUNK_SIZE = 1024;
+ Point< Real , Dim > points[ CHUNK_SIZE ];
+ Real values[ CHUNK_SIZE ];
+ size_t pointsRead;
+ while( ( pointsRead=pointStream->nextPoints( points , CHUNK_SIZE ) ) )
+ {
+ ThreadPool::Parallel_for( 0 , pointsRead , [&]( unsigned int thread , size_t j )
+ {
+ Point< Real , Dim > p = xForm * points[j];
+ bool inBounds = true;
+ for( int d=0 ; d1 ) inBounds = false;
+ if( inBounds ) values[j] = evaluator.values( xForm * points[j] , thread )[0];
+ else values[j] = (Real)nan( "" );
+ }
+ );
+ for( int j=0 ; j* tree , FILE* fp )
if( Verbose.set ) printf( "Vertices / Polygons: %llu / %llu\n" , (unsigned long long)( mesh.outOfCorePointCount()+mesh.inCorePoints.size() ) , (unsigned long long)mesh.polygonCount() );
std::vector< std::string > comments;
- if( !PlyWritePolygons< Vertex , node_index_type , Real , Dim >( OutMesh.value , &mesh , ASCII.set ? PLY_ASCII : PLY_BINARY_NATIVE , comments , XForm< Real , Dim+1 >::Identity() ) )
+ if( !PlyWritePolygons< Vertex , node_index_type , Real , Dim >( OutMesh.value , &mesh , ASCII.set ? PLY_ASCII : PLY_BINARY_NATIVE , comments , xForm.inverse() ) )
ERROR_OUT( "Could not write mesh to: " , OutMesh.value );
}
}
@@ -205,7 +239,8 @@ void _Execute( const FEMTree< Dim , Real >* tree , FILE* fp )
template< unsigned int Dim , class Real >
void Execute( FILE* fp , int degree , BoundaryType bType )
{
- FEMTree< Dim , Real > tree( fp , MEMORY_ALLOCATOR_BLOCK_SIZE );
+ XForm< Real , Dim+1 > xForm;
+ FEMTree< Dim , Real > tree( fp , xForm , MEMORY_ALLOCATOR_BLOCK_SIZE );
if( Verbose.set ) printf( "Leaf Nodes / Active Nodes / Ghost Nodes: %llu / %llu / %llu\n" , (unsigned long long)tree.leaves() , (unsigned long long)tree.nodes() , (unsigned long long)tree.ghostNodes() );
@@ -215,10 +250,10 @@ void Execute( FILE* fp , int degree , BoundaryType bType )
{
switch( degree )
{
- case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_FREE >::Signature >( &tree , fp ) ; break;
- case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_FREE >::Signature >( &tree , fp ) ; break;
- case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_FREE >::Signature >( &tree , fp ) ; break;
- case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_FREE >::Signature >( &tree , fp ) ; break;
+ case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_FREE >::Signature >( &tree , xForm , fp ) ; break;
+ case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_FREE >::Signature >( &tree , xForm , fp ) ; break;
+ case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_FREE >::Signature >( &tree , xForm , fp ) ; break;
+ case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_FREE >::Signature >( &tree , xForm , fp ) ; break;
default: ERROR_OUT( "Only B-Splines of degree 1 - 4 are supported" );
}
}
@@ -227,10 +262,10 @@ void Execute( FILE* fp , int degree , BoundaryType bType )
{
switch( degree )
{
- case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_NEUMANN >::Signature >( &tree , fp ) ; break;
- case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_NEUMANN >::Signature >( &tree , fp ) ; break;
- case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_NEUMANN >::Signature >( &tree , fp ) ; break;
- case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_NEUMANN >::Signature >( &tree , fp ) ; break;
+ case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_NEUMANN >::Signature >( &tree , xForm , fp ) ; break;
+ case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_NEUMANN >::Signature >( &tree , xForm , fp ) ; break;
+ case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_NEUMANN >::Signature >( &tree , xForm , fp ) ; break;
+ case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_NEUMANN >::Signature >( &tree , xForm , fp ) ; break;
default: ERROR_OUT( "Only B-Splines of degree 1 - 4 are supported" );
}
}
@@ -239,10 +274,10 @@ void Execute( FILE* fp , int degree , BoundaryType bType )
{
switch( degree )
{
- case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_DIRICHLET >::Signature >( &tree , fp ) ; break;
- case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_DIRICHLET >::Signature >( &tree , fp ) ; break;
- case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_DIRICHLET >::Signature >( &tree , fp ) ; break;
- case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_DIRICHLET >::Signature >( &tree , fp ) ; break;
+ case 1: _Execute< Dim , Real , FEMDegreeAndBType< 1 , BOUNDARY_DIRICHLET >::Signature >( &tree , xForm , fp ) ; break;
+ case 2: _Execute< Dim , Real , FEMDegreeAndBType< 2 , BOUNDARY_DIRICHLET >::Signature >( &tree , xForm , fp ) ; break;
+ case 3: _Execute< Dim , Real , FEMDegreeAndBType< 3 , BOUNDARY_DIRICHLET >::Signature >( &tree , xForm , fp ) ; break;
+ case 4: _Execute< Dim , Real , FEMDegreeAndBType< 4 , BOUNDARY_DIRICHLET >::Signature >( &tree , xForm , fp ) ; break;
default: ERROR_OUT( "Only B-Splines of degree 1 - 4 are supported" );
}
}
diff --git a/Src/EDTInHeat.cpp b/Src/EDTInHeat.cpp
index b5a2a73..40cb341 100644
--- a/Src/EDTInHeat.cpp
+++ b/Src/EDTInHeat.cpp
@@ -545,7 +545,7 @@ void _Execute( int argc , char* argv[] )
if( !fp ) ERROR_OUT( "Failed to open file for writing: " , Out.value );
FEMTree< Dim , Real >::WriteParameter( fp );
DenseNodeData< Real , IsotropicUIntPack< Dim , FEMSig > >::WriteSignatures( fp );
- tree.write( fp );
+ tree.write( fp , xForm );
edtSolution.write( fp );
fclose( fp );
}
diff --git a/Src/FEMTree.WeightedSamples.inl b/Src/FEMTree.WeightedSamples.inl
index a3e82b2..b263dda 100644
--- a/Src/FEMTree.WeightedSamples.inl
+++ b/Src/FEMTree.WeightedSamples.inl
@@ -116,7 +116,7 @@ void FEMTree< Dim , Real >::_getSampleDepthAndWeight( const DensityEstimator< We
{
Real oldWeight , newWeight;
oldWeight = newWeight = weight;
- while( newWeight<(Real)1. && temp->parent )
+ while( newWeight<(Real)1. && _localDepth(temp) )
{
temp=temp->parent;
oldWeight = newWeight;
diff --git a/Src/FEMTree.h b/Src/FEMTree.h
index bc8d110..71702c8 100644
--- a/Src/FEMTree.h
+++ b/Src/FEMTree.h
@@ -2129,13 +2129,13 @@ public:
static void ResetLocalMemoryUsage( void ){ _LocalMemoryUsage = 0; }
static double MemoryUsage( void );
FEMTree( size_t blockSize );
- FEMTree( FILE* fp , size_t blockSize );
+ FEMTree( FILE* fp , XForm< Real , Dim+1 > &xForm , size_t blockSize );
~FEMTree( void )
{
if( _tree ) for( int c=0 ; c<(1< xForm ) const;
static void WriteParameter( FILE* fp )
{
FEMTreeRealType realType;
@@ -2153,17 +2153,19 @@ public:
template< unsigned int Radius , class IsThickenNode , class ... DenseOrSparseNodeData > void thicken( IsThickenNode F , DenseOrSparseNodeData* ... data ){ thicken< Radius , Radius >( F , data ... ); }
template< unsigned int DensityDegree >
typename FEMTree::template DensityEstimator< DensityDegree >* setDensityEstimator( const std::vector< PointSample >& samples , LocalDepth splatDepth , Real samplesPerNode , int coDimension );
- template< unsigned int ... NormalSigs , unsigned int DensityDegree , class Data >
+ template< unsigned int ... DataSigs , unsigned int DensityDegree , class InData , class OutData >
+ SparseNodeData< OutData , UIntPack< DataSigs ... > > setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData & , Real & ) > ConversionAndBiasFunction );
+ template< unsigned int ... DataSigs , unsigned int DensityDegree , class InData , class OutData >
#if defined(_WIN32) || defined(_WIN64)
- SparseNodeData< Point< Real , Dim > , UIntPack< NormalSigs ... > > setNormalField( UIntPack< NormalSigs ... > , const std::vector< PointSample >& samples , const std::vector< Data >& normalData , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< Real ( Real ) > BiasFunction = []( Real ){ return 0.f; } );
+ SparseNodeData< OutData , UIntPack< DataSigs ... > > setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData& ) > ConversionFunction , std::function< Real ( InData ) > BiasFunction = []( InData ){ return 0.f; } );
#else // !_WIN32 && !_WIN64
- SparseNodeData< Point< Real , Dim > , UIntPack< NormalSigs ... > > setNormalField( UIntPack< NormalSigs ... > , const std::vector< PointSample >& samples , const std::vector< Data >& normalData , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< Real ( Real ) > BiasFunction = []( Real ){ return (Real)0; } );
+ SparseNodeData< OutData , UIntPack< DataSigs ... > > setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData& ) > ConversionFunction , std::function< Real ( InData ) > BiasFunction = []( InData ){ return (Real)0; } );
#endif // _WIN32 || _WIN64
template< unsigned int DataSig , bool CreateNodes , unsigned int DensityDegree , class Data >
SparseNodeData< Data , IsotropicUIntPack< Dim , DataSig > > setSingleDepthDataField( const std::vector< PointSample >& samples , const std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density );
template< unsigned int DataSig , bool CreateNodes , unsigned int DensityDegree , class Data >
- SparseNodeData< ProjectiveData< Data , Real > , IsotropicUIntPack< Dim , DataSig > > setDataField( const std::vector< PointSample >& samples , std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density , bool nearest=false );
+ SparseNodeData< ProjectiveData< Data , Real > , IsotropicUIntPack< Dim , DataSig > > setMultiDepthDataField( const std::vector< PointSample >& samples , std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density , bool nearest=false );
template< unsigned int MaxDegree , class HasDataFunctor , class ... DenseOrSparseNodeData > void finalizeForMultigrid( LocalDepth fullDepth , const HasDataFunctor F , DenseOrSparseNodeData* ... data );
template< unsigned int ... FEMSigs > DenseNodeData< Real , UIntPack< FEMSigs ... > > initDenseNodeData( UIntPack< FEMSigs ... > ) const;
diff --git a/Src/FEMTree.inl b/Src/FEMTree.inl
index e6cd90a..e4c2310 100644
--- a/Src/FEMTree.inl
+++ b/Src/FEMTree.inl
@@ -73,7 +73,7 @@ template< unsigned int Dim , class Real > FEMTree< Dim , Real >::FEMTree( size_t
memset( _refinableSigs , -1 , sizeof( _refinableSigs ) );
}
template< unsigned int Dim , class Real >
-FEMTree< Dim , Real >::FEMTree( FILE* fp , size_t blockSize ) : _nodeInitializer( *this )
+FEMTree< Dim , Real >::FEMTree( FILE* fp , XForm< Real , Dim+1 > &xForm , size_t blockSize ) : _nodeInitializer( *this )
{
if( blockSize )
{
@@ -87,6 +87,7 @@ FEMTree< Dim , Real >::FEMTree( FILE* fp , size_t blockSize ) : _nodeInitializer
Allocator< FEMTreeNode > *nodeAllocator = nodeAllocators.size() ? nodeAllocators[0] : NULL;
if( fp )
{
+ if( fread( xForm.coords , sizeof( Real ) , (Dim+1)*(Dim+1) , fp )!=(Dim+1)*(Dim+1) ) ERROR_OUT( "Failed to read transform" );
if( fread( &_depthOffset , sizeof( int ) , 1 , fp )!=1 ) ERROR_OUT( "Failed to read depth offset" );
_tree = FEMTreeNode::NewBrood( nodeAllocator , _nodeInitializer );
_tree->read( fp , nodeAllocator , _nodeInitializer );
@@ -113,8 +114,9 @@ FEMTree< Dim , Real >::FEMTree( FILE* fp , size_t blockSize ) : _nodeInitializer
_depthOffset = 0;
}
}
-template< unsigned int Dim , class Real > void FEMTree< Dim , Real >::write( FILE* fp ) const
+template< unsigned int Dim , class Real > void FEMTree< Dim , Real >::write( FILE* fp , XForm< Real , Dim+1 > xForm ) const
{
+ fwrite( xForm.coords , sizeof( Real ) , (Dim+1)*(Dim+1) , fp );
fwrite( &_depthOffset , sizeof( int ) , 1 , fp );
_tree->write( fp );
}
@@ -321,68 +323,81 @@ typename FEMTree< Dim , Real >::template DensityEstimator< DensityDegree >* FEMT
MemoryUsage();
return _density;
}
+
template< unsigned int Dim , class Real >
-template< unsigned int ... NormalSigs , unsigned int DensityDegree , class Data >
-SparseNodeData< Point< Real , Dim > , UIntPack< NormalSigs ... > > FEMTree< Dim , Real >::setNormalField( UIntPack< NormalSigs ... > , const std::vector< PointSample >& samples , const std::vector< Data >& normalData , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< Real ( Real ) > BiasFunction )
+template< unsigned int ... DataSigs , unsigned int DensityDegree , class InData , class OutData >
+SparseNodeData< OutData , UIntPack< DataSigs ... > > FEMTree< Dim , Real >::setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData& ) > ConversionFunction , std::function< Real ( InData ) > BiasFunction )
+{
+ std::function< bool ( InData , OutData & , Real & ) > ConversionAndBiasFunction = [&]( InData in , OutData &out , Real &bias )
+ {
+ if( ConversionFunction( in , out ) )
+ {
+ bias = BiasFunction( in );
+ return true;
+ }
+ else return false;
+ };
+ return setDataField( UIntPack< DataSigs ... >() , samples , data , density , pointWeightSum , ConversionAndBiasFunction );
+}
+template< unsigned int Dim , class Real >
+template< unsigned int ... DataSigs , unsigned int DensityDegree , class InData , class OutData >
+SparseNodeData< OutData , UIntPack< DataSigs ... > > FEMTree< Dim , Real >::setDataField( UIntPack< DataSigs ... > , const std::vector< PointSample >& samples , const std::vector< InData >& data , const DensityEstimator< DensityDegree >* density , Real& pointWeightSum , std::function< bool ( InData , OutData & , Real & ) > ConversionAndBiasFunction )
{
LocalDepth maxDepth = _spaceRoot->maxDepth();
typedef PointSupportKey< IsotropicUIntPack< Dim , DensityDegree > > DensityKey;
- typedef UIntPack< FEMSignature< NormalSigs >::Degree ... > NormalDegrees;
- typedef PointSupportKey< UIntPack< FEMSignature< NormalSigs >::Degree ... > > NormalKey;
+ typedef UIntPack< FEMSignature< DataSigs >::Degree ... > DataDegrees;
+ typedef PointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > > DataKey;
std::vector< DensityKey > densityKeys( ThreadPool::NumThreads() );
- std::vector< NormalKey > normalKeys( ThreadPool::NumThreads() );
- bool oneKey = DensityDegree==NormalDegrees::Min() && DensityDegree==NormalDegrees::Max();
+ std::vector< DataKey > dataKeys( ThreadPool::NumThreads() );
+ bool oneKey = DensityDegree==DataDegrees::Min() && DensityDegree==DataDegrees::Max();
for( size_t i=0 ; i , UIntPack< NormalSigs ... > > normalField;
+ SparseNodeData< OutData , UIntPack< DataSigs ... > > dataField;
Real _pointWeightSum = 0;
ThreadPool::Parallel_for( 0 , samples.size() , [&]( unsigned int thread , size_t i )
{
DensityKey& densityKey = densityKeys[ thread ];
- NormalKey& normalKey = normalKeys[ thread ];
+ DataKey& dataKey = dataKeys[ thread ];
const ProjectiveData< Point< Real , Dim > , Real >& sample = samples[i].sample;
if( sample.weight>0 )
{
- Point< Real , Dim > p = sample.data / sample.weight , n = normalData[i].template data<0>();
- Real l = (Real)Length( n );
- // It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
- if( !l ) return;
- Real confidence = l / sample.weight;
- n *= sample.weight / l;
- Real depthBias = BiasFunction( confidence );
- AddAtomic( weightSum , sample.weight );
- if( !_InBounds(p) )
- {
- WARN( "Point sample is out of bounds" );
- return;
- }
- Allocator< FEMTreeNode > *nodeAllocator = nodeAllocators.size() ? nodeAllocators[ thread ] : NULL;
+ Point< Real , Dim > p = sample.data / sample.weight;
+ InData in = data[i] / sample.weight;
+ OutData out;
-#if defined( __GNUC__ ) && __GNUC__ < 5
-#warning "you've got me gcc version<5"
- if( density ) AddAtomic( _pointWeightSum , _splatPointData< true , true , DensityDegree , Point< Real , Dim > >( nodeAllocator , *density , p , n , normalField , densityKey , oneKey ? *( (NormalKey*)&densityKey ) : normalKey , 0 , maxDepth , Dim , depthBias ) * sample.weight );
-#else // !__GNUC__ || __GNUC__ >=5
- if( density ) AddAtomic( _pointWeightSum , _splatPointData< true , true , DensityDegree , Point< Real , Dim > , NormalSigs ... >( nodeAllocator , *density , p , n , normalField , densityKey , oneKey ? *( (NormalKey*)&densityKey ) : normalKey , 0 , maxDepth , Dim , depthBias ) * sample.weight );
-#endif // __GNUC__ || __GNUC__ < 4
- else
+ Real depthBias;
+ if( !_InBounds(p) ) WARN( "Point sample is out of bounds" );
+ else if( ConversionAndBiasFunction( in , out , depthBias ) )
{
- Real width = (Real)( 1.0 / ( 1< *nodeAllocator = nodeAllocators.size() ? nodeAllocators[ thread ] : NULL;
#if defined( __GNUC__ ) && __GNUC__ < 5
-#warning "you've got me gcc version<5"
- _splatPointData< true , true , Point< Real , Dim > >( nodeAllocator , _leaf< true >( nodeAllocator , p , maxDepth ) , p , n / (Real)pow( width , Dim ) , normalField , oneKey ? *( (NormalKey*)&densityKey ) : normalKey );
+ #warning "you've got me gcc version<5"
+ if( density ) AddAtomic( _pointWeightSum , _splatPointData< true , true , DensityDegree , OutData >( nodeAllocator , *density , p , out , dataField , densityKey , oneKey ? *( (DataKey*)&densityKey ) : dataKey , 0 , maxDepth , Dim , depthBias ) * sample.weight );
#else // !__GNUC__ || __GNUC__ >=5
- _splatPointData< true , true , Point< Real , Dim > , NormalSigs ... >( nodeAllocator , _leaf< true >( nodeAllocator , p , maxDepth ) , p , n / (Real)pow( width , Dim ) , normalField , oneKey ? *( (NormalKey*)&densityKey ) : normalKey );
+ if( density ) AddAtomic( _pointWeightSum , _splatPointData< true , true , DensityDegree , OutData , DataSigs ... >( nodeAllocator , *density , p , out , dataField , densityKey , oneKey ? *( (DataKey*)&densityKey ) : dataKey , 0 , maxDepth , Dim , depthBias ) * sample.weight );
#endif // __GNUC__ || __GNUC__ < 4
- AddAtomic( _pointWeightSum , sample.weight );
+ else
+ {
+ Real width = (Real)( 1.0 / ( 1<( nodeAllocator , _leaf< true >( nodeAllocator , p , maxDepth ) , p , out / (Real)pow( width , Dim ) , dataField , oneKey ? *( (DataKey*)&densityKey ) : dataKey );
+#else // !__GNUC__ || __GNUC__ >=5
+ _splatPointData< true , true , OutData , DataSigs ... >( nodeAllocator , _leaf< true >( nodeAllocator , p , maxDepth ) , p , out / (Real)pow( width , Dim ) , dataField , oneKey ? *( (DataKey*)&densityKey ) : dataKey );
+#endif // __GNUC__ || __GNUC__ < 4
+ AddAtomic( _pointWeightSum , sample.weight );
+ }
}
}
}
);
pointWeightSum = _pointWeightSum / weightSum;
MemoryUsage();
- return normalField;
+ return dataField;
}
template< unsigned int Dim , class Real >
template< unsigned int DataSig , bool CreateNodes , unsigned int DensityDegree , class Data >
@@ -412,7 +427,7 @@ SparseNodeData< Data , IsotropicUIntPack< Dim , DataSig > > FEMTree< Dim , Real
}
template< unsigned int Dim , class Real >
template< unsigned int DataSig , bool CreateNodes , unsigned int DensityDegree , class Data >
-SparseNodeData< ProjectiveData< Data , Real > , IsotropicUIntPack< Dim , DataSig > > FEMTree< Dim , Real >::setDataField( const std::vector< PointSample >& samples , std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density , bool nearest )
+SparseNodeData< ProjectiveData< Data , Real > , IsotropicUIntPack< Dim , DataSig > > FEMTree< Dim , Real >::setMultiDepthDataField( const std::vector< PointSample >& samples , std::vector< Data >& sampleData , const DensityEstimator< DensityDegree >* density , bool nearest )
{
Allocator< FEMTreeNode > *nodeAllocator = nodeAllocators.size() ? nodeAllocators[0] : NULL;
LocalDepth maxDepth = _spaceRoot->maxDepth();
diff --git a/Src/Image.h b/Src/Image.h
index 563ff47..067147f 100644
--- a/Src/Image.h
+++ b/Src/Image.h
@@ -22,7 +22,7 @@ struct ImageReader
{
unsigned int channels;
ImageReader* reader = Get( fileName );
- width = reader->width() , height = reader->height();
+ width = reader->width() , height = reader->height() , channels = reader->channels();
if( channels!=1 && channels!=3 ) ERROR_OUT( "Requres one- or three-channel input" );
unsigned char* pixels = new unsigned char[ width*height*3 ];
unsigned char* pixelRow = new unsigned char[ width*channels];
diff --git a/Src/Ply.h b/Src/Ply.h
index 1255977..c04fa1d 100644
--- a/Src/Ply.h
+++ b/Src/Ply.h
@@ -30,8 +30,8 @@ WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
*/
-#ifndef __PLY_H__
-#define __PLY_H__
+#ifndef PLY_INCLUDED
+#define PLY_INCLUDED
#include
#include
@@ -569,4 +569,4 @@ int PlyWritePolygons( const char* fileName , CoredMeshData< Vertex , Index >* me
}
inline int PlyDefaultFileType( void ){ return PLY_ASCII; }
-#endif /* !__PLY_H__ */
+#endif // PLY_INCLUDED
diff --git a/Src/PlyFile.h b/Src/PlyFile.h
index c851c65..42f8c35 100644
--- a/Src/PlyFile.h
+++ b/Src/PlyFile.h
@@ -29,8 +29,8 @@ EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
*/
-#ifndef __PLY_FILE_H__
-#define __PLY_FILE_H__
+#ifndef PLY_FILE_INCLUDED
+#define PLY_FILE_INCLUDED
#include
#include
@@ -198,4 +198,4 @@ protected:
};
#include "PlyFile.inl"
-#endif /* !__PLY_FILE_H__ */
+#endif // PLY_FILE_INCLUDED
diff --git a/Src/PointInterpolant.cpp b/Src/PointInterpolant.cpp
new file mode 100644
index 0000000..bece832
--- /dev/null
+++ b/Src/PointInterpolant.cpp
@@ -0,0 +1,820 @@
+/*
+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.
+*/
+
+#include "PreProcessor.h"
+
+#undef USE_DOUBLE // If enabled, double-precesion is used
+#define WEIGHT_DEGREE 2 // The order of the B-Spline used to splat in the weights for density estimation
+#define DEFAULT_FEM_DEGREE 2 // The default finite-element degree
+#define DEFAULT_FEM_BOUNDARY BOUNDARY_FREE // The default finite-element boundary type
+#define DEFAULT_DIMENSION 2 // The dimension of the system
+
+#include
+#include
+#include
+#include
+#include "MyMiscellany.h"
+#include "CmdLineParser.h"
+#include "PPolynomial.h"
+#include "FEMTree.h"
+#include "Ply.h"
+#include "PointStreamData.h"
+#include "Image.h"
+
+MessageWriter messageWriter;
+
+cmdLineParameter< char* >
+ In( "in" ) ,
+ Out( "out" ) ,
+ TempDir( "tempDir" ) ,
+ Grid( "grid" ) ,
+ Tree( "tree" ) ,
+ Transform( "xForm" );
+
+cmdLineReadable
+ Performance( "performance" ) ,
+ ShowResidual( "showResidual" ) ,
+ PrimalGrid( "primalGrid" ) ,
+ ExactInterpolation( "exact" ) ,
+ InCore( "inCore" ) ,
+ NoValueConstraints( "noValues" ) ,
+ UseGradientConstraints( "useGradients" ) ,
+ NoComments( "noComments" ) ,
+ PolygonMesh( "polygonMesh" ) ,
+ NonManifold( "nonManifold" ) ,
+ NonLinearFit( "nonLinearFit" ) ,
+ ASCII( "ascii" ) ,
+ Verbose( "verbose" );
+
+cmdLineParameter< int >
+#ifndef FAST_COMPILE
+ Degree( "degree" , DEFAULT_FEM_DEGREE ) ,
+#endif // !FAST_COMPILE
+ Depth( "depth" , 8 ) ,
+ Iters( "iters" , 8 ) ,
+ FullDepth( "fullDepth" , 5 ) ,
+ BaseDepth( "baseDepth" , 5 ) ,
+ BaseVCycles( "baseVCycles" , 4 ) ,
+#ifndef FAST_COMPILE
+ BType( "bType" , DEFAULT_FEM_BOUNDARY+1 ) ,
+ Dimension( "dim" , DEFAULT_DIMENSION ) ,
+#endif // !FAST_COMPILE
+ MaxMemoryGB( "maxMemory" , 0 ) ,
+ ParallelType( "parallel" , (int)ThreadPool::OPEN_MP ) ,
+ ScheduleType( "schedule" , (int)ThreadPool::DefaultSchedule ) ,
+ ThreadChunkSize( "chunkSize" , (int)ThreadPool::DefaultChunkSize ) ,
+ Threads( "threads" , (int)std::thread::hardware_concurrency() );
+
+cmdLineParameter< float >
+ Scale( "scale" , 1.1f ) ,
+ Width( "width" , 0.f ) ,
+ CGSolverAccuracy( "cgAccuracy" , 1e-3f ) ,
+ IsoValue( "iso" , 0.f ) ,
+ ValueWeight ( "valueWeight" , 1000.f ) ,
+ GradientWeight( "gradientWeight" , 1.f ) ,
+ LapWeight ( "lapWeight" , 0.f ) ,
+ BiLapWeight ( "biLapWeight" , 1.f );
+
+cmdLineReadable* params[] =
+{
+#ifndef FAST_COMPILE
+ &Degree , &BType , &Dimension ,
+#endif // !FAST_COMPILE
+ &In , &Out , &Depth , &Transform ,
+ &Width ,
+ &Scale , &Verbose , &CGSolverAccuracy , &NoComments ,
+ &NonManifold , &PolygonMesh , &ASCII , &ShowResidual ,
+ &ValueWeight , &GradientWeight ,
+ &LapWeight , &BiLapWeight ,
+ &Grid , &Threads ,
+ &Tree ,
+ &FullDepth ,
+ &BaseDepth , &BaseVCycles ,
+ &Iters ,
+ &IsoValue ,
+ &PrimalGrid ,
+ &ExactInterpolation ,
+ &Performance ,
+ &MaxMemoryGB ,
+ &InCore ,
+ &ParallelType ,
+ &ScheduleType ,
+ &ThreadChunkSize ,
+ &NoValueConstraints ,
+ &UseGradientConstraints ,
+ &NonLinearFit ,
+ NULL
+};
+
+void ShowUsage(char* ex)
+{
+ printf( "Usage: %s\n" , ex );
+ printf( "\t --%s \n" , In.name );
+ printf( "\t[--%s ]\n" , Out.name );
+ printf( "\t[--%s ]\n" , Grid.name );
+ printf( "\t[--%s ]\n" , Tree.name );
+#ifndef FAST_COMPILE
+ printf( "\t[--%s =%d]\n" , Dimension.name , Dimension.value );
+ printf( "\t[--%s =%d]\n" , Degree.name , Degree.value );
+ printf( "\t[--%s =%d]\n" , BType.name , BType.value );
+ for( int i=0 ; i=%d]\n" , Depth.name , Depth.value );
+ printf( "\t[--%s ]\n" , Width.name );
+ printf( "\t[--%s =%d]\n" , FullDepth.name , FullDepth.value );
+ printf( "\t[--%s =%d]\n" , BaseDepth.name , BaseDepth.value );
+ printf( "\t[--%s =%d]\n" , BaseVCycles.name , BaseVCycles.value );
+ printf( "\t[--%s =%f]\n" , Scale.name , Scale.value );
+ printf( "\t[--%s =%.3e]\n" , ValueWeight.name , ValueWeight.value );
+ printf( "\t[--%s =%.3e]\n" , GradientWeight.name , GradientWeight.value );
+ printf( "\t[--%s =%.3e]\n" , LapWeight.name , LapWeight.value );
+ printf( "\t[--%s =%.3e]\n" , BiLapWeight.name , BiLapWeight.value );
+ printf( "\t[--%s =%d]\n" , Iters.name , Iters.value );
+ printf( "\t[--%s]\n" , ExactInterpolation.name );
+ printf( "\t[--%s =%d]\n" , Threads.name , Threads.value );
+ printf( "\t[--%s =%d]\n" , ParallelType.name , ParallelType.value );
+ for( size_t i=0 ; i=%d]\n" , ScheduleType.name , ScheduleType.value );
+ for( size_t i=0 ; i=%d]\n" , ThreadChunkSize.name , ThreadChunkSize.value );
+ printf( "\t[--%s =%g]\n" , CGSolverAccuracy.name , CGSolverAccuracy.value );
+ printf( "\t[--%s =%d]\n" , MaxMemoryGB.name , MaxMemoryGB.value );
+ printf( "\t[--%s =%f]\n" , IsoValue.name , IsoValue.value );
+ printf( "\t[--%s]\n" , NoValueConstraints.name );
+ printf( "\t[--%s]\n" , UseGradientConstraints.name );
+ printf( "\t[--%s]\n" , Performance.name );
+ printf( "\t[--%s]\n" , PrimalGrid.name );
+ printf( "\t[--%s]\n" , NoComments.name );
+ printf( "\t[--%s]\n" , PolygonMesh.name );
+ printf( "\t[--%s]\n" , NonManifold.name );
+ printf( "\t[--%s]\n" , NonLinearFit.name );
+ printf( "\t[--%s]\n" , ASCII.name );
+ printf( "\t[--%s]\n" , InCore.name );
+ printf( "\t[--%s]\n" , Verbose.name );
+}
+
+template< unsigned int Dim , class Real >
+struct FEMTreeProfiler
+{
+ FEMTree< Dim , Real >& tree;
+ double t;
+
+ FEMTreeProfiler( FEMTree< Dim , Real >& t ) : tree(t) { ; }
+ void start( void ){ t = Time() , FEMTree< Dim , Real >::ResetLocalMemoryUsage(); }
+ void print( const char* header ) const
+ {
+ FEMTree< Dim , Real >::MemoryUsage();
+ if( header ) printf( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+ else printf( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+ }
+ void dumpOutput( const char* header ) const
+ {
+ FEMTree< Dim , Real >::MemoryUsage();
+ if( header ) messageWriter( "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+ else messageWriter( "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+ }
+ void dumpOutput2( std::vector< std::string >& comments , const char* header ) const
+ {
+ FEMTree< Dim , Real >::MemoryUsage();
+ if( header ) messageWriter( comments , "%s %9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , header , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+ else messageWriter( comments , "%9.1f (s), %9.1f (MB) / %9.1f (MB) / %9.1f (MB)\n" , Time()-t , FEMTree< Dim , Real >::LocalMemoryUsage() , FEMTree< Dim , Real >::MaxMemoryUsage() , MemoryInfo::PeakMemoryUsageMB() );
+ }
+};
+
+template< class Real , unsigned int Dim >
+XForm< Real , Dim+1 > GetBoundingBoxXForm( Point< Real , Dim > min , Point< Real , Dim > max , Real scaleFactor )
+{
+ Point< Real , Dim > center = ( max + min ) / 2;
+ Real scale = max[0] - min[0];
+ for( int d=1 ; d( scale , max[d]-min[d] );
+ scale *= scaleFactor;
+ for( int i=0 ; i tXForm = XForm< Real , Dim+1 >::Identity() , sXForm = XForm< Real , Dim+1 >::Identity();
+ for( int i=0 ; i
+XForm< Real , Dim+1 > GetBoundingBoxXForm( Point< Real , Dim > min , Point< Real , Dim > max , Real width , Real scaleFactor , int& depth )
+{
+ // Get the target resolution (along the largest dimension)
+ Real resolution = ( max[0]-min[0] ) / width;
+ for( int d=1 ; d( resolution , ( max[d]-min[d] ) / width );
+ resolution *= scaleFactor;
+ depth = 0;
+ while( (1< center = ( max + min ) / 2;
+ Real scale = (1< tXForm = XForm< Real , Dim+1 >::Identity() , sXForm = XForm< Real , Dim+1 >::Identity();
+ for( int i=0 ; i
+XForm< Real , Dim+1 > GetPointXForm( InputPointStream< Real , Dim >& stream , Real width , Real scaleFactor , int& depth )
+{
+ Point< Real , Dim > min , max;
+ stream.boundingBox( min , max );
+ return GetBoundingBoxXForm( min , max , width , scaleFactor , depth );
+}
+template< class Real , unsigned int Dim >
+XForm< Real , Dim+1 > GetPointXForm( InputPointStream< Real , Dim >& stream , Real scaleFactor )
+{
+ Point< Real , Dim > min , max;
+ stream.boundingBox( min , max );
+ return GetBoundingBoxXForm( min , max , scaleFactor );
+}
+
+template< unsigned int Dim , typename Real , typename TotalPointSampleData > struct ValueAndGradientFromSample;
+
+template< unsigned int Dim , typename Real >
+struct ValueAndGradientFromSample< Dim , Real , MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > >
+{
+ typedef MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > TotalPointSampleData;
+ std::pair< Real , Point< Real , Dim > > operator()( TotalPointSampleData d ) const { return std::pair< Real , Point< Real , Dim > >( d.template data<0>() , d.template data<1>() ); }
+};
+
+template< unsigned int Dim , typename Real >
+struct ValueAndGradientFromSample< Dim , Real , MultiPointStreamData< Real , PointStreamValue< Real > > >
+{
+ typedef MultiPointStreamData< Real , PointStreamValue< Real > > TotalPointSampleData;
+ std::pair< Real , Point< Real , Dim > > operator()( TotalPointSampleData d ) const { return std::pair< Real , Point< Real , Dim > >( d.template data<0>() , Point< Real , Dim >() ); }
+};
+
+template< unsigned int Dim , typename Real >
+struct ValueAndGradientFromSample< Dim , Real , MultiPointStreamData< Real , PointStreamNormal< Real , Dim > > >
+{
+ typedef MultiPointStreamData< Real , PointStreamNormal< Real , Dim > > TotalPointSampleData;
+ std::pair< Real , Point< Real , Dim > > operator()( TotalPointSampleData d ) const { return std::pair< Real , Point< Real , Dim > >( (Real)0 , d.template data<0>() ); }
+};
+
+
+template< unsigned int Dim , typename Real , typename TotalPointSampleData > struct ConstraintDual;
+
+template< unsigned int Dim , typename Real >
+struct ConstraintDual< Dim , Real , MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > >
+{
+ typedef MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > TotalPointSampleData;
+ Real vWeight , gWeight;
+ ConstraintDual( Real v , Real g ) : vWeight(v) , gWeight(g) { }
+ CumulativeDerivativeValues< Real , Dim , 1 > operator()( const Point< Real , Dim >& p , const TotalPointSampleData& data ) const
+ {
+ Real value = data.template data<0>();
+ Point< Real , Dim > gradient = data.template data<1>();
+ CumulativeDerivativeValues< Real , Dim , 1 > cdv;
+ cdv[0] = value*vWeight;
+ for( int d=0 ; d
+struct ConstraintDual< Dim , Real , MultiPointStreamData< Real , PointStreamValue< Real > > >
+{
+ typedef MultiPointStreamData< Real , PointStreamValue< Real > > TotalPointSampleData;
+ Real vWeight , gWeight;
+ ConstraintDual( Real v , Real g ) : vWeight(v) , gWeight(g) { }
+ CumulativeDerivativeValues< Real , Dim , 1 > operator()( const Point< Real , Dim >& p , const TotalPointSampleData& data ) const
+ {
+ Real value = data.template data<0>();
+ CumulativeDerivativeValues< Real , Dim , 1 > cdv;
+ cdv[0] = value*vWeight;
+ return cdv;
+ }
+};
+template< unsigned int Dim , typename Real >
+struct ConstraintDual< Dim , Real , MultiPointStreamData< Real , PointStreamNormal< Real , Dim > > >
+{
+ typedef MultiPointStreamData< Real , PointStreamNormal< Real , Dim > > TotalPointSampleData;
+ Real vWeight , gWeight;
+ ConstraintDual( Real v , Real g ) : vWeight(v) , gWeight(g) { }
+ CumulativeDerivativeValues< Real , Dim , 1 > operator()( const Point< Real , Dim >& p , const TotalPointSampleData& data ) const
+ {
+ Point< Real , Dim > gradient = data.template data<0>();
+ CumulativeDerivativeValues< Real , Dim , 1 > cdv;
+ for( int d=0 ; d
+struct SystemDual
+{
+ CumulativeDerivativeValues< Real , Dim , 1 > weight;
+ SystemDual( Real v , Real g )
+ {
+ weight[0] = v;
+ for( int d=0 ; d operator()( Point< Real , Dim > p , const TotalPointSampleData& data , const CumulativeDerivativeValues< Real , Dim , 1 >& dValues ) const
+ {
+ return dValues * weight;
+ }
+ CumulativeDerivativeValues< double , Dim , 1 > operator()( Point< Real , Dim > p , const TotalPointSampleData& data , const CumulativeDerivativeValues< double , Dim , 1 >& dValues ) const
+ {
+ return dValues * weight;
+ };
+};
+template< unsigned int Dim , class TotalPointSampleData >
+struct SystemDual< Dim , double , TotalPointSampleData >
+{
+ typedef double Real;
+ CumulativeDerivativeValues< Real , Dim , 1 > weight;
+ SystemDual( Real v , Real g ) : weight( v , g , g , g ) { }
+ CumulativeDerivativeValues< Real , Dim , 1 > operator()( Point< Real , Dim > p , const TotalPointSampleData& data , const CumulativeDerivativeValues< Real , Dim , 1 >& dValues ) const
+ {
+ return dValues * weight;
+ }
+};
+
+template< typename Vertex , typename Real , unsigned int ... FEMSigs , typename TotalPointSampleData >
+void ExtractMesh( UIntPack< FEMSigs ... > , FEMTree< sizeof ... ( FEMSigs ) , Real >& tree , const DenseNodeData< Real , UIntPack< FEMSigs ... > >& solution , Real isoValue , const std::vector< typename FEMTree< sizeof ... ( FEMSigs ) , Real >::PointSample >* samples , std::function< void ( Vertex& , Point< Real , sizeof ... ( FEMSigs ) > , Real , TotalPointSampleData ) > SetVertex , std::vector< std::string > &comments , XForm< Real , sizeof...(FEMSigs)+1 > iXForm )
+{
+ static const int Dim = sizeof ... ( FEMSigs );
+ typedef UIntPack< FEMSigs ... > Sigs;
+ static const unsigned int DataSig = FEMDegreeAndBType< WEIGHT_DEGREE , BOUNDARY_FREE >::Signature;
+
+ FEMTreeProfiler< Dim , Real > profiler( tree );
+
+ char tempHeader[1024];
+ {
+ char tempPath[1024];
+ tempPath[0] = 0;
+ if( TempDir.set ) strcpy( tempPath , TempDir.value );
+ else SetTempDirectory( tempPath , sizeof(tempPath) );
+ if( strlen(tempPath)==0 ) sprintf( tempPath , ".%c" , FileSeparator );
+ if( tempPath[ strlen( tempPath )-1 ]==FileSeparator ) sprintf( tempHeader , "%sPR_" , tempPath );
+ else sprintf( tempHeader , "%s%cPR_" , tempPath , FileSeparator );
+ }
+
+ CoredMeshData< Vertex , node_index_type > *mesh;
+ if( InCore.set ) mesh = new CoredVectorMeshData< Vertex , node_index_type >();
+ else mesh = new CoredFileMeshData< Vertex , node_index_type >( tempHeader );
+ profiler.start();
+ typename IsoSurfaceExtractor< Dim , Real , Vertex >::IsoStats isoStats;
+#if defined( __GNUC__ ) && __GNUC__ < 5
+#warning "you've got me gcc version<5"
+ isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< TotalPointSampleData >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , (typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE >*)NULL , (SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > *)NULL , solution , isoValue , *mesh , SetVertex , NonLinearFit.set , !NonManifold.set , PolygonMesh.set , false );
+#else // !__GNUC__ || __GNUC__ >=5
+ isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< TotalPointSampleData >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , (typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE >*)NULL , NULL , solution , isoValue , *mesh , SetVertex , NonLinearFit.set , !NonManifold.set , PolygonMesh.set , false );
+#endif // __GNUC__ || __GNUC__ < 4
+ messageWriter( "Vertices / Polygons: %llu / %llu\n" , (unsigned long long)( mesh->outOfCorePointCount()+mesh->inCorePoints.size() ) , (unsigned long long)mesh->polygonCount() );
+ std::string isoStatsString = isoStats.toString() + std::string( "\n" );
+ messageWriter( isoStatsString.c_str() );
+ if( PolygonMesh.set ) profiler.dumpOutput2( comments , "# Got polygons:" );
+ else profiler.dumpOutput2( comments , "# Got triangles:" );
+
+ std::vector< std::string > noComments;
+ if( !PlyWritePolygons< Vertex , node_index_type , Real , Dim >( Out.value , mesh , ASCII.set ? PLY_ASCII : PLY_BINARY_NATIVE , NoComments.set ? noComments : comments , iXForm ) )
+ ERROR_OUT( "Could not write mesh to: " , Out.value );
+ delete mesh;
+}
+
+template< typename Real , unsigned int Dim >
+void WriteGrid( ConstPointer( Real ) values , int res , const char *fileName )
+{
+ int resolution = 1;
+ for( int d=0 ; d avgs( ThreadPool::NumThreads() , 0 );
+ ThreadPool::Parallel_for( 0 , resolution , [&]( unsigned int thread , size_t i ){ avgs[thread] += values[i]; } );
+ for( unsigned int t=0 ; t stds( ThreadPool::NumThreads() , 0 );
+ ThreadPool::Parallel_for( 0 , resolution , [&]( unsigned int thread , size_t i ){ stds[thread] += ( values[i] - avg ) * ( values[i] - avg ); } );
+ for( unsigned int t=0 ; t [0,255]\n" , avg - 2*std , avg + 2*std );
+
+ unsigned char *pixels = new unsigned char[ resolution*3 ];
+ ThreadPool::Parallel_for( 0 , resolution , [&]( unsigned int , size_t i )
+ {
+ Real v = (Real)std::min< Real >( (Real)1. , std::max< Real >( (Real)-1. , ( values[i] - avg ) / (2*std ) ) );
+ v = (Real)( ( v + 1. ) / 2. * 256. );
+ unsigned char color = (unsigned char )std::min< Real >( (Real)255. , std::max< Real >( (Real)0. , v ) );
+ for( int c=0 ; c<3 ; c++ ) pixels[i*3+c ] = color;
+ }
+ );
+ ImageWriter::Write( fileName , pixels , res , res , 3 );
+ delete[] pixels;
+ }
+ else
+ {
+
+ FILE *fp = fopen( fileName , "wb" );
+ if( !fp ) ERROR_OUT( "Failed to open grid file for writing: " , fileName );
+ else
+ {
+ fwrite( &res , sizeof(int) , 1 , fp );
+ if( typeid(Real)==typeid(float) ) fwrite( values , sizeof(float) , resolution , fp );
+ else
+ {
+ float *fValues = new float[resolution];
+ for( int i=0 ; i
+void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
+{
+ static const int Dim = sizeof ... ( FEMSigs );
+ typedef UIntPack< FEMSigs ... > Sigs;
+ typedef UIntPack< FEMSignature< FEMSigs >::Degree ... > Degrees;
+ typedef UIntPack< FEMDegreeAndBType< WEIGHT_DEGREE , DerivativeBoundary< FEMSignature< FEMSigs >::BType , 1 >::BType >::Signature ... > DataSigs;
+ typedef typename FEMTree< Dim , Real >::template DensityEstimator< WEIGHT_DEGREE > DensityEstimator;
+ typedef typename FEMTree< Dim , Real >::template InterpolationInfo< Real , 1 > InterpolationInfo;
+ typedef InputPointStreamWithData< Real , Dim , TotalPointSampleData > InputPointStream;
+ typedef TransformedInputPointStreamWithData< Real , Dim , TotalPointSampleData > XInputPointStream;
+ std::vector< std::string > comments;
+ messageWriter( comments , "***********************************************\n" );
+ messageWriter( comments , "***********************************************\n" );
+ messageWriter( comments , "** Running Point Interpolant (Version %s) **\n" , VERSION );
+ messageWriter( comments , "***********************************************\n" );
+ messageWriter( comments , "***********************************************\n" );
+ if( !Threads.set ) messageWriter( comments , "Running with %d threads\n" , Threads.value );
+
+ ThreadPool::Init( (ThreadPool::ParallelType)ParallelType.value , Threads.value );
+
+ XForm< Real , Dim+1 > xForm , iXForm;
+ if( Transform.set )
+ {
+ FILE* fp = fopen( Transform.value , "r" );
+ if( !fp )
+ {
+ WARN( "Could not read x-form from: " , Transform.value );
+ xForm = XForm< Real , Dim+1 >::Identity();
+ }
+ else
+ {
+ for( int i=0 ; i::Identity();
+
+ char str[1024];
+ for( int i=0 ; params[i] ; i++ )
+ if( params[i]->set )
+ {
+ params[i]->writeValue( str );
+ if( strlen( str ) ) messageWriter( comments , "\t--%s %s\n" , params[i]->name , str );
+ else messageWriter( comments , "\t--%s\n" , params[i]->name );
+ }
+
+ double startTime = Time();
+
+ FEMTree< Dim , Real > tree( MEMORY_ALLOCATOR_BLOCK_SIZE );
+ FEMTreeProfiler< Dim , Real > profiler( tree );
+
+ if( Depth.set && Width.value>0 )
+ {
+ WARN( "Both --" , Depth.name , " and --" , Width.name , " set, ignoring --" , Width.name );
+ Width.value = 0;
+ }
+
+ size_t pointCount;
+
+ std::vector< typename FEMTree< Dim , Real >::PointSample >* samples = new std::vector< typename FEMTree< Dim , Real >::PointSample >();
+ std::vector< TotalPointSampleData >* sampleData = NULL;
+
+ // Read in the samples
+ {
+ profiler.start();
+ InputPointStream* pointStream;
+ char* ext = GetFileExtension( In.value );
+ sampleData = new std::vector< TotalPointSampleData >();
+ std::vector< std::pair< Point< Real , Dim > , TotalPointSampleData > > inCorePoints;
+ if( InCore.set )
+ {
+ InputPointStream *_pointStream;
+ if ( !strcasecmp( ext , "bnpts" ) ) _pointStream = new BinaryInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::ReadBinary );
+ else if( !strcasecmp( ext , "ply" ) ) _pointStream = new PLYInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::PlyReadProperties() , TotalPointSampleData::PlyReadNum , TotalPointSampleData::ValidPlyReadProperties );
+ else _pointStream = new ASCIIInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::ReadASCII );
+ Point< Real , Dim > p;
+ TotalPointSampleData d;
+ while( _pointStream->nextPoint( p , d ) ) inCorePoints.push_back( std::pair< Point< Real , Dim > , TotalPointSampleData >( p , d ) );
+ delete _pointStream;
+
+ pointStream = new MemoryInputPointStreamWithData< Real , Dim , TotalPointSampleData >( inCorePoints.size() , &inCorePoints[0] );
+ }
+ else
+ {
+ if ( !strcasecmp( ext , "bnpts" ) ) pointStream = new BinaryInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::ReadBinary );
+ else if( !strcasecmp( ext , "ply" ) ) pointStream = new PLYInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::PlyReadProperties() , TotalPointSampleData::PlyReadNum , TotalPointSampleData::ValidPlyReadProperties );
+ else pointStream = new ASCIIInputPointStreamWithData< Real , Dim , TotalPointSampleData >( In.value , TotalPointSampleData::ReadASCII );
+ }
+ delete[] ext;
+ typename TotalPointSampleData::Transform _xForm( xForm );
+ XInputPointStream _pointStream( [&]( Point< Real , Dim >& p , TotalPointSampleData& d ){ p = xForm*p , d = _xForm(d); } , *pointStream );
+ if( Width.value>0 ) xForm = GetPointXForm< Real , Dim >( _pointStream , Width.value , (Real)( Scale.value>0 ? Scale.value : 1. ) , Depth.value ) * xForm;
+ else xForm = Scale.value>0 ? GetPointXForm< Real , Dim >( _pointStream , (Real)Scale.value ) * xForm : xForm;
+ {
+ typename TotalPointSampleData::Transform _xForm( xForm );
+ XInputPointStream _pointStream( [&]( Point< Real , Dim >& p , TotalPointSampleData& d ){ p = xForm*p , d = _xForm(d); } , *pointStream );
+ auto ProcessData = []( const Point< Real , Dim >& p , TotalPointSampleData& d ){ return (Real)1.; };
+ pointCount = FEMTreeInitializer< Dim , Real >::template Initialize< TotalPointSampleData >( tree.spaceRoot() , _pointStream , Depth.value , *samples , *sampleData , true , tree.nodeAllocators.size() ? tree.nodeAllocators[0] : NULL , tree.initializer() , ProcessData );
+ }
+ iXForm = xForm.inverse();
+ delete pointStream;
+
+ messageWriter( "Input Points / Samples: %llu / %llu\n" , pointCount , (unsigned long long)samples->size() );
+ profiler.dumpOutput2( comments , "# Read input into tree:" );
+ }
+
+ DenseNodeData< Real , Sigs > solution;
+ {
+ DenseNodeData< Real , Sigs > constraints;
+ InterpolationInfo* iInfo = NULL;
+ int solveDepth = Depth.value;
+
+ tree.resetNodeIndices();
+
+ // Prepare for multigrid
+ {
+ profiler.start();
+ tree.template finalizeForMultigrid< Degrees::Max() >( FullDepth.value , []( const typename FEMTree< Dim , Real >::FEMTreeNode * ){ return true; } );
+ profiler.dumpOutput2( comments , "# Finalized tree:" );
+ }
+
+ // Add the interpolation constraints
+ {
+ profiler.start();
+ if( ExactInterpolation.set ) iInfo = FEMTree< Dim , Real >::template InitializeExactPointAndDataInterpolationInfo< Real , TotalPointSampleData , 1 >( tree , *samples , GetPointer( *sampleData ) , ConstraintDual< Dim , Real , TotalPointSampleData >( (Real)ValueWeight.value , (Real)GradientWeight.value ) , SystemDual< Dim , Real , TotalPointSampleData >( (Real)ValueWeight.value , (Real)GradientWeight.value ) , true , false );
+ else iInfo = FEMTree< Dim , Real >::template InitializeApproximatePointAndDataInterpolationInfo< Real , TotalPointSampleData , 1 >( tree , *samples , GetPointer( *sampleData ) , ConstraintDual< Dim , Real , TotalPointSampleData >( (Real)ValueWeight.value , (Real)GradientWeight.value ) , SystemDual< Dim , Real , TotalPointSampleData >( (Real)ValueWeight.value , (Real)GradientWeight.value ) , true , 1 );
+ constraints = tree.initDenseNodeData( Sigs() );
+ tree.addInterpolationConstraints( constraints , solveDepth , *iInfo );
+ profiler.dumpOutput2( comments , "#Set point constraints:" );
+ }
+
+ messageWriter( "Leaf Nodes / Active Nodes / Ghost Nodes: %llu / %llu / %llu\n" , (unsigned long long)tree.leaves() , (unsigned long long)tree.nodes() , (unsigned long long)tree.ghostNodes() );
+ messageWriter( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage())/(1<<20) );
+
+ // Solve the linear system
+ {
+ profiler.start();
+ typename FEMTree< Dim , Real >::SolverInfo sInfo;
+ sInfo.cgDepth = 0 , sInfo.cascadic = true , sInfo.vCycles = 1 , sInfo.iters = Iters.value , sInfo.cgAccuracy = CGSolverAccuracy.value , sInfo.verbose = Verbose.set , sInfo.showResidual = ShowResidual.set , sInfo.showGlobalResidual = SHOW_GLOBAL_RESIDUAL_NONE , sInfo.sliceBlockSize = 1;
+ sInfo.baseDepth = BaseDepth.value , sInfo.baseVCycles = BaseVCycles.value;
+ typename FEMIntegrator::template System< Sigs , IsotropicUIntPack< Dim , 2 > > F( { 0. , (double)LapWeight.value , (double)BiLapWeight.value } );
+ solution = tree.solveSystem( Sigs() , F , constraints , solveDepth , sInfo , iInfo );
+ profiler.dumpOutput2( comments , "# Linear system solved:" );
+ if( iInfo ) delete iInfo , iInfo = NULL;
+ }
+ }
+
+ if( Verbose.set )
+ {
+ typename FEMTree< Dim , Real >::template MultiThreadedEvaluator< Sigs , 1 > evaluator( &tree , solution );
+ std::pair< double , double > valueStat(0,0) , gradientStat(0,0);
+ std::vector< std::pair< double , double > > valueStats( ThreadPool::NumThreads() , std::pair< double , double >(0,0) ) , gradientStats( ThreadPool::NumThreads() , std::pair< double , double >(0,0) );
+ ValueAndGradientFromSample< Dim , Real , TotalPointSampleData > valueAndGradientFromSample;
+ ThreadPool::Parallel_for( 0 , samples->size() , [&]( unsigned int thread , size_t j )
+ {
+ ProjectiveData< Point< Real , Dim > , Real >& sample = (*samples)[j].sample;
+ Real w = sample.weight;
+ if( w>0 )
+ {
+ CumulativeDerivativeValues< Real , Dim , 1 > values = evaluator.values( sample.data / sample.weight , thread , (*samples)[j].node );
+ Real value = values[0];
+ Point< Real , Dim > gradient;
+ for( int d=0 ; d > valueAndGradient = valueAndGradientFromSample( (*sampleData)[j] / w );
+ valueStats[ thread ].first += ( value - valueAndGradient.first ) * ( value - valueAndGradient.first ) * w;
+ valueStats[ thread ].second += ( value * value + valueAndGradient.first * valueAndGradient.first ) * w;
+ gradientStats[ thread ].first += Point< Real , Dim >::SquareNorm( gradient - valueAndGradient.second ) * w;
+ gradientStats[ thread ].second += ( Point< Real , Dim >::SquareNorm( gradient ) + Point< Real , Dim >::SquareNorm( valueAndGradient.second ) ) * w;
+ }
+ }
+ );
+ for( unsigned int t=0 ; t0 && GradientWeight.value>0 ) messageWriter( "Value / Gradient Error: %g / %g\n" , (Real)sqrt( valueStat.first / valueStat.second ) , (Real)sqrt( gradientStat.first / gradientStat.second ) );
+ else if( ValueWeight.value>0 ) messageWriter( "Value Error: %g\n" , (Real)sqrt( valueStat.first / valueStat.second ) );
+ else if( GradientWeight.value>0 ) messageWriter( "Gradient Error: %g\n" , (Real)sqrt( gradientStat.first / gradientStat.second ) );
+ }
+
+ delete samples , samples = NULL;
+ delete sampleData , sampleData = NULL;
+
+
+ if( Tree.set )
+ {
+ FILE* fp = fopen( Tree.value , "wb" );
+ if( !fp ) ERROR_OUT( "Failed to open file for writing: " , Tree.value );
+ FEMTree< Dim , Real >::WriteParameter( fp );
+ DenseNodeData< Real , Sigs >::WriteSignatures( fp );
+ tree.write( fp , xForm );
+ solution.write( fp );
+ fclose( fp );
+ }
+
+ if( Grid.set )
+ {
+ int res = 0;
+ profiler.start();
+ Pointer( Real ) values = tree.template regularGridEvaluate< true >( solution , res , -1 , PrimalGrid.set );
+ size_t resolution = 1;
+ for( int d=0 ; d( values , res , Grid.value );
+ DeletePointer( values );
+ if( Verbose.set )
+ {
+ printf( "Transform:\n" );
+ for( int i=0 ; i Vertex;
+ std::function< void ( Vertex& , Point< Real , Dim > , Real , TotalPointSampleData ) > SetVertex = []( Vertex& v , Point< Real , Dim > p , Real , TotalPointSampleData ){ v.point = p; };
+ ExtractMesh< Vertex >( UIntPack< FEMSigs ... >() , tree , solution , IsoValue.value , samples , SetVertex , comments , iXForm );
+ }
+
+ messageWriter( comments , "# Total Solve: %9.1f (s), %9.1f (MB)\n" , Time()-startTime , FEMTree< Dim , Real >::MaxMemoryUsage() );
+}
+
+template< class Real , unsigned int ... FEMSigs >
+void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
+{
+ static const int Dim = sizeof ... ( FEMSigs );
+ if ( !UseGradientConstraints.set ) Execute< Real , MultiPointStreamData< Real , PointStreamValue< Real > > >( argc , argv , UIntPack< FEMSigs ... >() );
+ else if( NoValueConstraints.set ) Execute< Real , MultiPointStreamData< Real , PointStreamNormal< Real , Dim > > >( argc , argv , UIntPack< FEMSigs ... >() );
+ else Execute< Real , MultiPointStreamData< Real , PointStreamValue< Real > , PointStreamNormal< Real , Dim > > >( argc , argv , UIntPack< FEMSigs ... >() );
+}
+
+#ifndef FAST_COMPILE
+template< unsigned int Dim , class Real >
+void Execute( int argc , char* argv[] )
+{
+ switch( BType.value )
+ {
+ case BOUNDARY_FREE+1:
+ {
+ switch( Degree.value )
+ {
+// case 1: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 1 , BOUNDARY_FREE >::Signature >() );
+ case 2: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 2 , BOUNDARY_FREE >::Signature >() );
+ case 3: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 3 , BOUNDARY_FREE >::Signature >() );
+// case 4: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 4 , BOUNDARY_FREE >::Signature >() );
+ default: ERROR_OUT( "Only B-Splines of degree 1 - 3 are supported" );
+ }
+ }
+ case BOUNDARY_NEUMANN+1:
+ {
+ switch( Degree.value )
+ {
+// case 1: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 1 , BOUNDARY_NEUMANN >::Signature >() );
+ case 2: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 2 , BOUNDARY_NEUMANN >::Signature >() );
+ case 3: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 3 , BOUNDARY_NEUMANN >::Signature >() );
+// case 4: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 4 , BOUNDARY_NEUMANN >::Signature >() );
+ default: ERROR_OUT( "Only B-Splines of degree 1 - 3 are supported" );
+ }
+ }
+ case BOUNDARY_DIRICHLET+1:
+ {
+ switch( Degree.value )
+ {
+// case 1: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 1 , BOUNDARY_DIRICHLET >::Signature >() );
+ case 2: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 2 , BOUNDARY_DIRICHLET >::Signature >() );
+ case 3: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 3 , BOUNDARY_DIRICHLET >::Signature >() );
+// case 4: return Execute< Real >( argc , argv , IsotropicUIntPack< Dim , FEMDegreeAndBType< 4 , BOUNDARY_DIRICHLET >::Signature >() );
+ default: ERROR_OUT( "Only B-Splines of degree 1 - 3 are supported" );
+ }
+ }
+ default: ERROR_OUT( "Not a valid boundary type: " , BType.value );
+ }
+}
+#endif // !FAST_COMPILE
+
+int main( int argc , char* argv[] )
+{
+ Timer timer;
+#ifdef USE_SEG_FAULT_HANDLER
+ WARN( "using seg-fault handler" );
+ StackTracer::exec = argv[0];
+ signal( SIGSEGV , SignalHandler );
+#endif // USE_SEG_FAULT_HANDLER
+#ifdef ARRAY_DEBUG
+ WARN( "Array debugging enabled" );
+#endif // ARRAY_DEBUG
+
+ cmdLineParse( argc-1 , &argv[1] , params );
+#if 0
+ if( !In.set || !Out.set ) ERROR_OUT( "Need input and output" );
+ unsigned int width , height;
+ unsigned char *pixels = ImageReader::ReadColor( In.value , width , height );
+ FILE *fp = fopen( Out.value , "wb" );
+ if( !fp ) ERROR_OUT( "Failed to open file for reading: %s" , Out.value );
+ for( int i=0 ; i<10000 ; i++ )
+ {
+ int x = rand() % width , y = rand() % height;
+ double gray = (double)( pixels[ 3*(y*width+x) + 0 ] + pixels[ 3*(y*width+x) + 1 ] + pixels[ 3*(y*width+x) + 2 ] ) / ( 255. * 3 );
+ fprintf( fp , "%d %d %f\n" , x , y , gray );
+ }
+ fclose( fp );
+#else
+
+ if( MaxMemoryGB.value>0 ) SetPeakMemoryMB( MaxMemoryGB.value<<10 );
+ ThreadPool::DefaultChunkSize = ThreadChunkSize.value;
+ ThreadPool::DefaultSchedule = (ThreadPool::ScheduleType)ScheduleType.value;
+ messageWriter.echoSTDOUT = Verbose.set;
+
+ if( !In.set )
+ {
+ ShowUsage( argv[0] );
+ return 0;
+ }
+ if( NoValueConstraints.set ) ValueWeight.value = 0;
+ if( !UseGradientConstraints.set ) GradientWeight.value = 0;
+
+ if( ValueWeight.value<0 ) ERROR_OUT( "Value weight must be non-negative: " , ValueWeight.value , "> 0" );
+ if( GradientWeight.value<0 ) ERROR_OUT( "Gradient weight must be non-negative: " , GradientWeight.value , "> 0" );
+ if( !ValueWeight.value && !GradientWeight.value ) ERROR_OUT( "Either value or gradient weight must be positive" );
+
+ if( LapWeight.value<0 ) ERROR_OUT( "Laplacian weight must be non-negative: " , LapWeight.value , " > 0" );
+ if( BiLapWeight.value<0 ) ERROR_OUT( "Bi-Laplacian weight must be non-negative: " , BiLapWeight.value , " > 0" );
+ if( !LapWeight.value && !BiLapWeight.value ) ERROR_OUT( "Eiter Laplacian or bi-Laplacian weight must be positive" );
+
+ if( BaseDepth.value>FullDepth.value )
+ {
+ if( BaseDepth.set ) WARN( "Base depth must be smaller than full depth: " , BaseDepth.value , " <= " , FullDepth.value );
+ BaseDepth.value = FullDepth.value;
+ }
+
+#ifdef USE_DOUBLE
+ typedef double Real;
+#else // !USE_DOUBLE
+ typedef float Real;
+#endif // USE_DOUBLE
+
+#ifdef FAST_COMPILE
+ static const int Dimension = DIMENSION;
+ static const int Degree = DEFAULT_FEM_DEGREE;
+ static const BoundaryType BType = DEFAULT_FEM_BOUNDARY;
+ typedef IsotropicUIntPack< Dimension , FEMDegreeAndBType< Degree , BType >::Signature > FEMSigs;
+ WARN( "Compiled for degree-" , Degree , ", boundary-" , BoundaryNames[ BType ] , ", " , sizeof(Real)==4 ? "single" : "double" , "-precision _only_" );
+ Execute< Real >( argc , argv , FEMSigs() );
+#else // !FAST_COMPILE
+ if( Dimension.value==2 ) Execute< 2 , Real >( argc , argv );
+ else if( Dimension.value==3 ) Execute< 3 , Real >( argc , argv );
+ else ERROR_OUT( "Only Degrees 2 and 3 are supported" );
+#endif // FAST_COMPILE
+ if( Performance.set )
+ {
+ printf( "Time (Wall/CPU): %.2f / %.2f\n" , timer.wallTime() , timer.cpuTime() );
+ printf( "Peak Memory (MB): %d\n" , MemoryInfo::PeakMemoryUsageMB() );
+ }
+#endif
+ return EXIT_SUCCESS;
+}
diff --git a/Src/PointStream.inl b/Src/PointStream.inl
index a14e130..5cafe61 100644
--- a/Src/PointStream.inl
+++ b/Src/PointStream.inl
@@ -150,8 +150,8 @@ void PLYInputPointStream< Real , Dim >::reset( void )
{
foundVertices = true;
_pCount = num_elems , _pIdx = 0;
- for( int i=0 ; i::ReadComponents ; i++ )
- if( !_ply->get_property( elem_name , &(PlyVertex< Real , Dim >::Properties()[i]) ) ) ERROR_OUT( "Failed to find property in ply file: %s" , PlyVertex< Real , Dim >::Properties()[i].name );
+ for( int i=0 ; i::PlyReadNum ; i++ )
+ if( !_ply->get_property( elem_name , &(PlyVertex< Real , Dim >::PlyReadProperties()[i]) ) ) ERROR_OUT( "Failed to find property in ply file: %s" , PlyVertex< Real , Dim >::PlyReadProperties()[i].name );
}
for( int j=0 ; j
protected:
XForm< Real , Dim+1 > _xForm;
};
- static void readASCII( FILE* fp , PointStreamPosition& p )
+ static void ReadASCII( FILE* fp , PointStreamPosition& p )
{
float f;
for( int i=0 ; i >
};
static void ReadASCII( FILE* fp , PointStreamColor& p )
{
- unsigned char c[3];
- if( fscanf( fp , " %c %c %c " , &c[0] , &c[1] , &c[2] )!=3 ) ERROR_OUT( "Failed to read color" );
+ unsigned int c[3];
+ if( fscanf( fp , " %d %d %d " , &c[0] , &c[1] , &c[2] )!=3 ) ERROR_OUT( "Failed to read color" );
p.psData[0] = (Real)c[0] , p.psData[1] = (Real)c[1] , p.psData[2] = (Real)c[2];
};
static void ReadBinary( FILE* fp , PointStreamColor& p )
diff --git a/Src/PoissonRecon.cpp b/Src/PoissonRecon.cpp
index cb2e8e9..a5659ed 100644
--- a/Src/PoissonRecon.cpp
+++ b/Src/PoissonRecon.cpp
@@ -338,7 +338,7 @@ void ExtractMesh( UIntPack< FEMSigs ... > , std::tuple< SampleData ... > , FEMTr
typename IsoSurfaceExtractor< Dim , Real , Vertex >::IsoStats isoStats;
if( sampleData )
{
- SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL );
+ SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setMultiDepthDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL );
for( const RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) )
{
ProjectiveData< TotalPointSampleData , Real >* clr = _sampleData( n );
@@ -582,9 +582,27 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
{
profiler.start();
normalInfo = new SparseNodeData< Point< Real , Dim > , NormalSigs >();
- if( ConfidenceBias.value>0 ) *normalInfo = tree.setNormalField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , [&]( Real conf ){ return (Real)( log( conf ) * ConfidenceBias.value / log( 1<<(Dim-1) ) ); } );
- else *normalInfo = tree.setNormalField( NormalSigs() , *samples , *sampleData , density , pointWeightSum );
- ThreadPool::Parallel_for( 0 , normalInfo->size() , [&]( unsigned int , size_t i ){ (*normalInfo)[i] *= (Real)-1.; } );
+ std::function< bool ( TotalPointSampleData , Point< Real , Dim >& ) > ConversionFunction = []( TotalPointSampleData in , Point< Real , Dim > &out )
+ {
+ Point< Real , Dim > n = in.template data<0>();
+ Real l = (Real)Length( n );
+ // It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
+ if( !l ) return false;
+ out = n / l;
+ return true;
+ };
+ std::function< bool ( TotalPointSampleData , Point< Real , Dim >& , Real & ) > ConversionAndBiasFunction = []( TotalPointSampleData in , Point< Real , Dim > &out , Real &bias )
+ {
+ Point< Real , Dim > n = in.template data<0>();
+ Real l = (Real)Length( n );
+ // It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
+ if( !l ) return false;
+ out = n / l;
+ bias = (Real)( log( l ) * ConfidenceBias.value / log( 1<<(Dim-1) ) );
+ return true;
+ };
+ if( ConfidenceBias.value>0 ) *normalInfo = tree.setDataField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , ConversionAndBiasFunction );
+ else *normalInfo = tree.setDataField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , ConversionFunction ); ThreadPool::Parallel_for( 0 , normalInfo->size() , [&]( unsigned int , size_t i ){ (*normalInfo)[i] *= (Real)-1.; } );
profiler.dumpOutput2( comments , "# Got normal field:" );
messageWriter( "Point weight / Estimated Area: %g / %g\n" , pointWeightSum , pointCount*pointWeightSum );
}
@@ -672,7 +690,7 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
if( !fp ) ERROR_OUT( "Failed to open file for writing: " , Tree.value );
FEMTree< Dim , Real >::WriteParameter( fp );
DenseNodeData< Real , Sigs >::WriteSignatures( fp );
- tree.write( fp );
+ tree.write( fp , xForm );
solution.write( fp );
fclose( fp );
}
diff --git a/Src/PreProcessor.h b/Src/PreProcessor.h
index b7d76b9..c07380d 100644
--- a/Src/PreProcessor.h
+++ b/Src/PreProcessor.h
@@ -44,7 +44,7 @@ DAMAGE.
#define USE_DEEP_TREE_NODES // Chances are that if you are using big data, you want to support a tree with depth>15.
#endif // BIG_DATA
-#define VERSION "11.02" // The version of the code
+#define VERSION "12.00" // The version of the code
#define MEMORY_ALLOCATOR_BLOCK_SIZE 1<<12 // The chunk size for memory allocation
#endif // PRE_PROCESSOR_INCLUDED
\ No newline at end of file
diff --git a/Src/SSDRecon.cpp b/Src/SSDRecon.cpp
index 3068f40..3a0b55c 100644
--- a/Src/SSDRecon.cpp
+++ b/Src/SSDRecon.cpp
@@ -358,7 +358,7 @@ void ExtractMesh( UIntPack< FEMSigs ... > , std::tuple< SampleData ... > , FEMTr
typename IsoSurfaceExtractor< Dim , Real , Vertex >::IsoStats isoStats;
if( sampleData )
{
- SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL );
+ SparseNodeData< ProjectiveData< TotalPointSampleData , Real > , IsotropicUIntPack< Dim , DataSig > > _sampleData = tree.template setMultiDepthDataField< DataSig , false >( *samples , *sampleData , (DensityEstimator*)NULL );
for( const RegularTreeNode< Dim , FEMTreeNodeData , depth_and_offset_type >* n = tree.tree().nextNode() ; n ; n=tree.tree().nextNode( n ) )
{
ProjectiveData< TotalPointSampleData , Real >* clr = _sampleData( n );
@@ -603,8 +603,27 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
{
profiler.start();
normalInfo = new SparseNodeData< Point< Real , Dim > , NormalSigs >();
- if( ConfidenceBias.value>0 ) *normalInfo = tree.setNormalField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , [&]( Real conf ){ return (Real)( log( conf ) * ConfidenceBias.value / log( 1<<(Dim-1) ) ); } );
- else *normalInfo = tree.setNormalField( NormalSigs() , *samples , *sampleData , density , pointWeightSum );
+ std::function< bool ( TotalPointSampleData , Point< Real , Dim >& ) > ConversionFunction = []( TotalPointSampleData in , Point< Real , Dim > &out )
+ {
+ Point< Real , Dim > n = in.template data<0>();
+ Real l = (Real)Length( n );
+ // It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
+ if( !l ) return false;
+ out = n / l;
+ return true;
+ };
+ std::function< bool ( TotalPointSampleData , Point< Real , Dim >& , Real & ) > ConversionAndBiasFunction = []( TotalPointSampleData in , Point< Real , Dim > &out , Real &bias )
+ {
+ Point< Real , Dim > n = in.template data<0>();
+ Real l = (Real)Length( n );
+ // It is possible that the samples have non-zero normals but there are two co-located samples with negative normals...
+ if( !l ) return false;
+ out = n / l;
+ bias = (Real)( log( l ) * ConfidenceBias.value / log( 1<<(Dim-1) ) );
+ return true;
+ };
+ if( ConfidenceBias.value>0 ) *normalInfo = tree.setDataField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , ConversionAndBiasFunction );
+ else *normalInfo = tree.setDataField( NormalSigs() , *samples , *sampleData , density , pointWeightSum , ConversionFunction );
profiler.dumpOutput2( comments , "# Got normal field:" );
messageWriter( "Point weight / Estimated Area: %g / %g\n" , pointWeightSum , pointCount*pointWeightSum );
}
@@ -674,7 +693,7 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
if( !fp ) ERROR_OUT( "Failed to open file for writing: " , Tree.value );
FEMTree< Dim , Real >::WriteParameter( fp );
DenseNodeData< Real , Sigs >::WriteSignatures( fp );
- tree.write( fp );
+ tree.write( fp , xForm );
solution.write( fp );
fclose( fp );
}
diff --git a/Src/SparseMatrix.h b/Src/SparseMatrix.h
index 44a9c2c..ee2487c 100644
--- a/Src/SparseMatrix.h
+++ b/Src/SparseMatrix.h
@@ -25,8 +25,8 @@ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
DAMAGE.
*/
-#ifndef __SPARSEMATRIX_HPP
-#define __SPARSEMATRIX_HPP
+#ifndef SPARSE_MATRIX_INCLUDED
+#define SPARSE_MATRIX_INCLUDED
#include "SparseMatrixInterface.h"
#include "Array.h"
@@ -151,4 +151,4 @@ public:
};
#include "SparseMatrix.inl"
-#endif /* __SPARSEMATRIX_HPP */
+#endif // SPARSE_MATRIX_INCLUDED