[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