Version 18.40

mkazhdan · mkazhdan · commit 913f13b90517 · 2024-10-30T12:14:32.000-04:00
diff --git a/README.md b/README.md
@@ -1,4 +1,4 @@
-<center><h2>Adaptive Multigrid Solvers (Version 18.35)</h2></center>
+<center><h2>Adaptive Multigrid Solvers (Version 18.40)</h2></center>
 <center>
 <a href="#LINKS">links</a>
 <a href="#COMPILATION">compilation</a>
@@ -29,10 +29,11 @@ This code-base was born from the Poisson Surface Reconstruction code. It has evo
 <a href="https://www.cs.jhu.edu/~misha/MyPapers/CGF23.pdf">[Kazhdan and Hoppe, 2023]</a>
 <br>
 <b>Executables: </b>
-<a href="https://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version18.35/AdaptiveSolvers.x64.zip">Win64</a><br>
+<a href="https://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version18.40/AdaptiveSolvers.x64.zip">Win64</a><br>
 <b>Source Code:</b>
-<a href="https://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version18.35/AdaptiveSolvers.zip">ZIP</a> <a href="https://github.com/mkazhdan/PoissonRecon">GitHub</a><br>
+<a href="https://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version18.40/AdaptiveSolvers.zip">ZIP</a> <a href="https://github.com/mkazhdan/PoissonRecon">GitHub</a><br>
 <b>Older Versions:</b>
+<a href="https://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version18.35/">V18.35</a>,
 <a href="https://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version18.31/">V18.31</a>,
 <a href="https://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version18.30/">V18.30</a>,
 <a href="https://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version18.20/">V18.20</a>,
@@ -1054,11 +1055,11 @@ If both the <b>--keep</b> flag and the <b>--fraction</b> flag are set, the <b>--
 In addition to executables, the reconstruction code can be interfaced into through the functionality implemented in <CODE>Reconstructors.h</CODE> and <CODE>Extrapolator.h</CODE>
 Using the functionality requires requires choosing a finite element type and defining input oriented-point stream and output vertex and face streams. In the descriptions below, the template parameter <CODE>Real</CODE> is the floating point type used to represent data (typically <code>float</code>) and <CODE>Dim</CODE> is the integer dimension of the space (fixed at <CODE>Dim</CODE>=3).
 <UL>
-<LI>The template parameter <CODE>FEMSig</CODE> describes the 1D finite element type, which is a composite of the degree of the finite element and the boundary conditions it satisfies. Given an integer valued <CODE>Degree</CODE> and boundary type <CODE>BType</CODE> (one of <CODE>BOUNDARY_FREE</CODE>, <CODE>BOUNDARY_DIRICHLET</CODE>, and <CODE>BOUNDARY_NEUMANN</CODE> defined in <CODE>BSplineData.h</CODE>), the signature is defined as (<CODE>Reconstruction.example.cppy</CODE> line 308):
+<LI>The template parameter <CODE>FEMSig</CODE> describes the 1D finite element type, which is a composite of the degree of the finite element and the boundary conditions it satisfies. Given an integer valued <CODE>Degree</CODE> and boundary type <CODE>BType</CODE> (one of <CODE>BOUNDARY_FREE</CODE>, <CODE>BOUNDARY_DIRICHLET</CODE>, and <CODE>BOUNDARY_NEUMANN</CODE> defined in <CODE>BSplineData.h</CODE>), the signature is defined as:
 <PRE>
 <CODE>static const unsigned int FEMSig = FEMDegreeAndBType&lt; Degree , BoundaryType &gt;::Signature;</CODE>
 </PRE>
-<LI>The template parameter <CODE>FEMSIgs</CODE> describes the tensor-product finite element type, typically defined as an "isotropic" element with the same 1D finite element type across all <CODE>Dim</CODE> dimensions. It is defined as (<CODE>Reconstruction.example.cppy</CODE> line 314):
+<LI>The template parameter <CODE>FEMSIgs</CODE> describes the tensor-product finite element type, typically defined as an "isotropic" element with the same 1D finite element type across all <CODE>Dim</CODE> dimensions. It is defined as:
 <PRE>
 <CODE>using FEMSigs = IsotropicUIntPack&lt; Dim , FEMSig &gt;;</CODE>
 </PRE>
@@ -1117,22 +1118,22 @@ This method returns the extrapolated value at the prescribed position.
 <UL>
 These steps can be found in the <code>Reconstruction.example.cpp</code> code.
 <UL>
-<LI>The finite-elements signature is created in line 308.
-<LI>An input sample stream generating a specified number of random points on the surface of the sphere is defined in lines 85-122 and constructed in line 374.
-<LI>An output polygon stream that pushes the polygon to an <code>std::vector</code> of <code>std::vector&lt; int &gt;</code>s is defined in lines 213-230 and constructed in line 384.
-<LI>An output vertex stream that pushes just the position information to an <code>std::vector</code> of <code>Real</code>s is desfined in lines 223-244 and constructed in line 385.
-<LI>The reconstructor is constructed in line 377.
-<LI>The level-set extraction is performed on line 388.
-<LI>The evaluator is created on line 334.
-<LI>The evaluator is used to query the values and gradients of the implicit function in line 327.
-<LI>The extrapolator is constructed on line 406.
-<LI>The extrapolator is evaluated at the vertex positions at line 416.
+<LI>The finite-elements signatures are created in lines 308 and 311.
+<LI>An input sample stream generating a specified number of random points on the surface of the sphere is defined in lines 85-122 and constructed in line 379.
+<LI>An output polygon stream that pushes the polygon to an <code>std::vector</code> of <code>std::vector&lt; int &gt;</code>s is defined in lines 213-230 and constructed in line 389.
+<LI>An output vertex stream that pushes just the position information to an <code>std::vector</code> of <code>Real</code>s is desfined in lines 223-244 and constructed in line 390.
+<LI>The implicit representation is computed in line 382.
+<LI>The level-set extraction is performed on line 393.
+<LI>The evaluator is created on line 340.
+<LI>The evaluator is used to query the values and gradients of the implicit function in line 333.
+<LI>The extrapolator is constructed on line 411.
+<LI>The extrapolator is evaluated at the vertex positions at line 421.
 </UL>
-Note that a similar approach can be used to perform the <A HREF="http://mesh.brown.edu/ssd/">Smoothed Signed Distance</A> reconstruction (lines 452 and 453).<BR>
+Note that a similar approach can be used to perform the <A HREF="http://mesh.brown.edu/ssd/">Smoothed Signed Distance</A> reconstruction (lines 459 and 460).<BR>
 The approach also supports reconstruction of meshes with auxiliary information like color, with the only constraint that the auxiliary data type supports the computation affine combinations (e.g. the <CODE>RGBColor</CODE> type defined in lines 67-82). The auxiliary inform is derived in one of two ways:
 <OL>
-<LI> As part of the reconstruction process, so that the level-set extraction phase also sets the auxiliary information. (Lines 343-470)
-<LI> Independently by constructing the extrapolation field from the samples and then evaluating at the positions of the level-set vertices. (Lines 397-418)
+<LI> As part of the reconstruction process, so that the level-set extraction phase also sets the auxiliary information. (Lines 349-375)
+<LI> Independently by constructing the extrapolation field from the samples and then evaluating at the positions of the level-set vertices. (Lines 400-427)
 </OL>
 </UL>
 </DL>
@@ -1712,6 +1713,11 @@ Similarly, to reduce compilation times, support for specific degrees can be remo
 <LI> Removed the confidence bias option.
 </OL>
 
+<a href="https://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version18.31/">Version 18.40</a>:
+<OL>
+<LI> Added support for extrapolator accumulation.
+</OL>
+
 </DETAILS>
 
 <hr>
diff --git a/Src/Extrapolator.h b/Src/Extrapolator.h
@@ -98,8 +98,10 @@ namespace PoissonRecon
 			XForm< Real , Dim+1 > worldToUnitCubeTransform( void ) const { return _worldToUnitCube; }
 
 			AuxData operator()( unsigned int t , Point< Real , Dim > p ){ return _value(t,p); }
+			AuxData operator()(                  Point< Real , Dim > p ){ return _value(0,p); }
 
-			AuxData operator()( Point< Real , Dim > p ){ return operator()( 0 , p ); }
+			void evaluate( unsigned int t , Point< Real , Dim > p , AuxData &data ){ return _evaluate( t , p , data ); }
+			void evaluate(                  Point< Real , Dim > p , AuxData &data ){ return _evaluate( 0 , p , data ); }
 
 		protected:
 			typename FEMTree< Dim , Real >::template MultiThreadedSparseEvaluator< IsotropicUIntPack< Dim , DataSig > , ProjectiveData< AuxData , Real > > *_auxEvaluator = nullptr;
@@ -113,6 +115,23 @@ namespace PoissonRecon
 				_auxEvaluator->addValue( q , pData );
 				return pData.value();
 			}
+
+			void _evaluate( unsigned int t , Point< Real , Dim > p , AuxData &data )
+			{
+				Point< Real , Dim > q = _worldToUnitCube * p;
+				for( unsigned int d=0 ; d<Dim ; d++ ) if( q[d]<0 || q[d]>1 ) throw OutOfUnitCubeException(p,q);
+
+				Real weight = (Real)0.;
+				data *= 0;
+
+				auto Accumulation = [&weight,&data]( const ProjectiveData< AuxData , Real > &pData , Real scale )
+					{
+						weight += pData.weight * scale;
+						data += pData.data * scale;
+					};
+				_auxEvaluator->accumulate( q , Accumulation );
+				if( weight ) data /= weight;
+			}
 		};
 
 		template< typename Real , unsigned int Dim , typename AuxData , unsigned int DataDegree >
diff --git a/Src/FEMTree.Evaluation.inl b/Src/FEMTree.Evaluation.inl
@@ -550,15 +550,43 @@ void FEMTree< Dim , Real >::MultiThreadedSparseEvaluator< UIntPack< FEMSigs ...
 	_tree->template _addEvaluation< T , SparseNodeData< T , FEMSignatures > , 0 >( _coefficients , p , _tree->_globalToLocal( node->depth() ) , *_pointEvaluator , nKey , t );
 }
 
+template< unsigned int Dim , class Real >
+template< unsigned int ... FEMSigs , typename T >
+template< typename AccumulationFunctor/*=std::function< void ( const T & , Real s ) > */ >
+void FEMTree< Dim , Real >::MultiThreadedSparseEvaluator< UIntPack< FEMSigs ... > , T >::accumulate( Point< Real , Dim > p , AccumulationFunctor &Accumulate , int thread , const FEMTreeNode *node )
+{
+	if( !node ) node = _tree->leaf( p );
+	ConstPointSupportKey< FEMDegrees >& nKey = _pointNeighborKeys[thread];
+	nKey.getNeighbors( node );
+	_tree->template _accumulate< T , SparseNodeData< T , FEMSignatures > , 0 , AccumulationFunctor >( _coefficients , p , _tree->_globalToLocal( node->depth() ) , *_pointEvaluator , nKey , Accumulate );
+}
+
 template< unsigned int Dim , class Real >
 template< class V , class Coefficients , unsigned int D , unsigned int ... DataSigs >
 void FEMTree< Dim , Real >::_addEvaluation( const Coefficients& coefficients , Point< Real , Dim > p , const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , V &value ) const
 {
-	_addEvaluation< V , Coefficients , D , DataSigs ... >( coefficients , p , _globalToLocal( dataKey.depth() ) , pointEvaluator , dataKey , value );
+	auto AF = [&value]( const V &w , Real s ){ value += w * s; };
+	_accumulate< V , Coefficients , D , decltype(AF) , DataSigs ... >( coefficients , p , pointEvaluator , dataKey , AF );
 }
+
 template< unsigned int Dim , class Real >
 template< class V , class Coefficients , unsigned int D , unsigned int ... DataSigs >
 void FEMTree< Dim , Real >::_addEvaluation( const Coefficients& coefficients , Point< Real , Dim > p , LocalDepth pointDepth , const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , V &value ) const
+{
+	auto AF = [&value]( const V &w , Real s ){ value += w * s; };
+	_accumulate< V , Coefficients , D , decltype(AF) , DataSigs... >( coefficients , p , pointDepth , pointEvaluator , dataKey , AF );
+}
+
+template< unsigned int Dim , class Real >
+template< typename V , typename Coefficients , unsigned int D , typename AccumulationFunctor , unsigned int ... DataSigs >
+void FEMTree< Dim , Real >::_accumulate( const Coefficients& coefficients , Point< Real , Dim > p , const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , AccumulationFunctor &Accumulate ) const
+{
+	_accumulate< V , Coefficients , D , AccumulationFunctor , DataSigs ... >( coefficients , p , _globalToLocal( dataKey.depth() ) , pointEvaluator , dataKey , Accumulate );
+}
+
+template< unsigned int Dim , class Real >
+template< typename V , typename Coefficients , unsigned int D , typename AccumulationFunctor , unsigned int ... DataSigs >
+void FEMTree< Dim , Real >::_accumulate( const Coefficients& coefficients , Point< Real , Dim > p , LocalDepth pointDepth , const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , AccumulationFunctor &Accumulate ) const
 {
 	typedef UIntPack< BSplineSupportSizes< FEMSignature< DataSigs >::Degree >::SupportSize ... > SupportSizes;
 	PointEvaluatorState< UIntPack< DataSigs ... > , ZeroUIntPack< Dim > > state;
@@ -582,7 +610,7 @@ void FEMTree< Dim , Real >::_addEvaluation( const Coefficients& coefficients , P
 			if( v )
 			{
 				LocalDepth d ; LocalOffset off ; _localDepthAndOffset( nodes[i] , d , off );
-				value += (*v) * (Real)state.value( off , derivatives );
+				Accumulate( *v , (Real)state.value( off , derivatives ) );
 			}
 		}
 	}
diff --git a/Src/FEMTree.h b/Src/FEMTree.h
@@ -2302,8 +2302,14 @@ namespace PoissonRecon
 		template< bool CreateNodes , bool ThreadSafe , unsigned int WeightDegree , class V , unsigned int ... DataSigs > Point< Real , 2 >      _splatPointData( Allocator< FEMTreeNode > *nodeAllocator , const DensityEstimator< WeightDegree >& densityWeights , Real minDepthCutoff ,                     Point< Real , Dim > point , V v , SparseNodeData< V , UIntPack< DataSigs ... > >& data , PointSupportKey< IsotropicUIntPack< Dim , WeightDegree > >& weightKey , PointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , LocalDepth minDepth , std::function< int ( Point< Real , Dim > ) > &pointDepthFunctor , int dim , Real depthBias );
 		template< bool CreateNodes , bool ThreadSafe , unsigned int WeightDegree , class V , unsigned int ... DataSigs > Point< Real , 2 > _multiSplatPointData( Allocator< FEMTreeNode > *nodeAllocator , const DensityEstimator< WeightDegree >* densityWeights , FEMTreeNode* node , Point< Real , Dim > point , V v , SparseNodeData< V , UIntPack< DataSigs ... > >& data , PointSupportKey< IsotropicUIntPack< Dim , WeightDegree > >& weightKey , PointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , int dim );
 		template< unsigned int WeightDegree , class V , unsigned int ... DataSigs > Real _nearestMultiSplatPointData( const DensityEstimator< WeightDegree >* densityWeights , FEMTreeNode* node , Point< Real , Dim > point , V v , SparseNodeData< V , UIntPack< DataSigs ... > >& data , PointSupportKey< IsotropicUIntPack< Dim , WeightDegree > >& weightKey , int dim=Dim );
-		template< class V , class Coefficients , unsigned int D , unsigned int ... DataSigs > void _addEvaluation( const Coefficients& coefficients , Point< Real , Dim > p , const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , V &value ) const;
-		template< class V , class Coefficients , unsigned int D , unsigned int ... DataSigs > void _addEvaluation( const Coefficients& coefficients , Point< Real , Dim > p , LocalDepth pointDepth , const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , V &value ) const;
+		template< class V , class Coefficients , unsigned int D , unsigned int ... DataSigs >
+		void _addEvaluation( const Coefficients& coefficients , Point< Real , Dim > p ,                         const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , V &value ) const;
+		template< class V , class Coefficients , unsigned int D , unsigned int ... DataSigs >
+		void _addEvaluation( const Coefficients& coefficients , Point< Real , Dim > p , LocalDepth pointDepth , const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , V &value ) const;
+		template< typename V , typename Coefficients , unsigned int D , typename AccumulationFunctor /*=std::function< void ( const V & , Real s ) > */ , unsigned int ... DataSigs >
+		void _accumulate( const Coefficients& coefficients , Point< Real , Dim > p ,                         const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , AccumulationFunctor &AF ) const;
+		template< typename V , typename Coefficients , unsigned int D , typename AccumulationFunctor /*=std::function< void ( const V & , Real s ) > */ , unsigned int ... DataSigs >
+		void _accumulate( const Coefficients& coefficients , Point< Real , Dim > p , LocalDepth pointDepth , const PointEvaluator< UIntPack< DataSigs ... > , IsotropicUIntPack< Dim , D > >& pointEvaluator , const ConstPointSupportKey< UIntPack< FEMSignature< DataSigs >::Degree ... > >& dataKey , AccumulationFunctor &AF ) const;
 
 	public:
 
@@ -2513,6 +2519,8 @@ namespace PoissonRecon
 			MultiThreadedSparseEvaluator( const FEMTree* tree , const SparseNodeData< T , FEMSignatures >& coefficients , int threads=std::thread::hardware_concurrency() );
 			~MultiThreadedSparseEvaluator( void ){ if( _pointEvaluator ) delete _pointEvaluator; }
 			void addValue( Point< Real , Dim > p , T &t , int thread=0 , const FEMTreeNode* node=NULL );
+			template< typename AccumulationFunctor/*=std::function< void ( const T & , Real s ) > */ >
+			void accumulate( Point< Real , Dim > p , AccumulationFunctor &Accumulate , int thread=0 , const FEMTreeNode* node=NULL );
 		};
 
 	protected:
diff --git a/Src/PoissonRecon.cpp b/Src/PoissonRecon.cpp
@@ -102,7 +102,6 @@ CmdLineParameter< float >
 	Confidence( "confidence" , 0.f ) ,
 	CGSolverAccuracy( "cgAccuracy" , 1e-3f ) ,
 	LowDepthCutOff( "lowDepthCutOff" , 0.f ) ,
-	TargetValue( "targetValue" , 0.5f ) ,
 	PointWeight( "pointWeight" );
 
 CmdLineReadable* params[] =
@@ -140,7 +139,6 @@ CmdLineReadable* params[] =
 	&ThreadChunkSize ,
 	&LowDepthCutOff ,
 	&AlignmentDir ,
-	&TargetValue ,
 	&GridCoordinates ,
 	NULL
 };
@@ -168,7 +166,6 @@ void ShowUsage(char* ex)
 	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 <interpolation weight>=%.3e * <b-spline degree>]\n" , PointWeight.name , Reconstructor::Poisson::WeightMultiplier * Reconstructor::Poisson::DefaultFEMDegree );
-	printf( "\t[--%s <target value>=%f]\n" , TargetValue.name , TargetValue.value );
 	printf( "\t[--%s <iterations>=%d]\n" , Iters.name , Iters.value );
 	printf( "\t[--%s]\n" , ExactInterpolation.name );
 	printf( "\t[--%s <pull factor>=%f]\n" , DataX.name , DataX.value );
@@ -327,7 +324,6 @@ void Execute( const AuxDataFactory &auxDataFactory )
 	sParams.pointWeight = (Real)PointWeight.value;
 	sParams.samplesPerNode = (Real)SamplesPerNode.value;
 	sParams.cgSolverAccuracy = (Real)CGSolverAccuracy.value;
-	sParams.targetValue = (Real)TargetValue.value;
 	sParams.perLevelDataScaleFactor = (Real)DataX.value;
 	sParams.depth = (unsigned int)Depth.value;
 	sParams.baseDepth = (unsigned int)BaseDepth.value;
diff --git a/Src/PreProcessor.h b/Src/PreProcessor.h
@@ -46,7 +46,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 ADAPTIVE_SOLVERS_VERSION "18.35"		// The version of the code
+#define ADAPTIVE_SOLVERS_VERSION "18.40"		// The version of the code
 #define MEMORY_ALLOCATOR_BLOCK_SIZE 1<<12		// The chunk size for memory allocation
 
 #define ADAPTIVE_PADDING						// Only pushes padding points deep enough so that they are "close" to the slab in terms of units at that depth
diff --git a/Src/Reconstructors.h b/Src/Reconstructors.h
