[surface] Poisson surface reconstruction bug in Octree< Real >::UpdateWeightContribution

[menelausyjl]

In the changelist of the original authors' website, they stated that there was a bug in sample contribution scaling, but it is still kept in PCL.

The fixed version(v7.0):

template< class Real >
int Octree< Real >::UpdateWeightContribution( std::vector< Real >& kernelDensityWeights , TreeOctNode* node , const Point3D<Real>& position , typename TreeOctNode::NeighborKey3& neighborKey , Real weight )
{
	typename TreeOctNode::Neighbors3& neighbors = neighborKey.setNeighbors( node );
	if( (int)kernelDensityWeights.size()<TreeNodeData::NodeCount ) kernelDensityWeights.resize( TreeNodeData::NodeCount , 0 );
	double x , dxdy , dx[DIMENSION][3] , width;
	Point3D< Real > center;
	Real w;
	node->centerAndWidth( center , w );
	width = w;

	for( int i=0 ; i<DIMENSION ; i++ )
	{
		x = ( center[i] - position[i] - width ) / width;
		dx[i][0] = 1.125 + 1.500*x + 0.500*x*x;
		dx[i][1] = -0.25 - 2.*x - x*x;
		dx[i][2] = 1. - dx[i][1] - dx[i][0];
	}

	// Suppose that the samples are uniformly placed along the middle of the three slices.
	// Then splatting the points we get coefficients:
	//		0.125 / 0.75 / 0.125 across the three slices.
	// Sampling at the center slice we get:
	//		0.125^2 + 0.75^2 + 0.125^2 = 19/32
	const double SAMPLE_SCALE = 1. / ( 0.125 * 0.125 + 0.75 * 0.75 + 0.125 * 0.125 );
	weight *= (Real)SAMPLE_SCALE;

	for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ )
	{
		dxdy = dx[0][i] * dx[1][j] * weight;
		TreeOctNode** _neighbors = neighbors.neighbors[i][j];
		for( int k=0 ; k<3 ; k++ ) if( _neighbors[k] ) kernelDensityWeights[ _neighbors[k]->nodeData.nodeIndex ] += Real( dxdy * dx[2][k] );
	}
	return 0;
}

The PCL version:

    template<int Degree>
    int Octree<Degree>::NonLinearUpdateWeightContribution( TreeOctNode* node , const Point3D<Real>& position , Real weight )
    {
      TreeOctNode::Neighbors3& neighbors = neighborKey.setNeighbors( node );
      double x,dxdy,dx[DIMENSION][3];
      double width;
      Point3D<Real> center;
      Real w;
      node->centerAndWidth( center , w );
      width=w;
      const double SAMPLE_SCALE = 1. / ( 0.125 * 0.125 + 0.75 * 0.75 + 0.125 * 0.125 );

      for( int i=0 ; i<DIMENSION ; i++ )
      {
        x = ( center[i] - position[i] - width ) / width;
        dx[i][0] = 1.125 + 1.500*x + 0.500*x*x;
        x = ( center[i] - position[i] ) / width;
        dx[i][1] = 0.750           -       x*x;
        dx[i][2] = 1. - dx[i][1] - dx[i][0];
        // Note that we are splatting along a co-dimension one manifold, so uniform point samples
        // do not generate a unit sample weight.
        dx[i][0] *= SAMPLE_SCALE;
      }

      for( int i=0 ; i<3 ; i++ ) for( int j=0 ; j<3 ; j++ )
      {
        dxdy = dx[0][i] * dx[1][j] * weight;
        for( int k=0 ; k<3 ; k++ ) if( neighbors.neighbors[i][j][k] )
          neighbors.neighbors[i][j][k]->nodeData.centerWeightContribution += Real( dxdy * dx[2][k] );
      }
      return 0;
    }

[mvieth]

Hi, yes the Poisson code that PCL uses is quite outdated, but unfortunately nobody found time to update it yet