We use the poisson reconstruction on some point clouds with normals. Sometimes we get a spike in the output that is annoying to have to deal with. It is interesting to note that this bug usually occurs on point clouds with large holes / empty areas. Complete point clouds don't generally have this problem. They generally appear in a region of the mesh that lies on the boundary of the point cloud.
The attached screenshot shows an example of the spike. The left, denser portion of the mesh corresponds to data in the point cloud. The right, less-dense portion of the mesh did not have any corresponding data in the point cloud. The spike has a single vertex at its tip and three faces that intersect three other triangles on the mesh. It is like a very long tetrahedron but the base does not have a triangle defined for it. If the vertex were relocated to someplace on the surface of the mesh in the middle of its base then the mesh would be smooth, like it should be.
A coworker who knows more about these things than I do suspected the following code as a culprit:
File: MultiGridOctreeData.inl
Function: template<int degree="">
int Octree<degree>::GetRoot(const RootInfo& ri,const Real& isoValue,Point3D<real> & position,hash_map<long long,std::pair\<real,point3d\<real=""> > >& normalHash,const int& nonLinearFit)</long></real></degree></int>
Code:
if(rCount && nonLinearFit) {averageRoot/=rCount;}
else {averageRoot=Real((x0-isoValue)/(x0-x1));}
The concern is the "else" block: if x0 is close to x1 then it is almost division by zero and averageRoot becomes very large. To test, we added some debug code to print the value of "averageRoot" in the "else" block. Normally, averageRoot would be somewhere in the range (-1, 1) but sometimes it would be very large. The spike in the screenshot had a value of about 25, for example. We have seen values as high as 200 (which creates very long spikes). There seems to be a correlation between the number of high-magnitude values of averageRoot and the size & number of spikes observed in the mesh. Any calculated value of averageRoot that has a magnitude greater than 2 appears to cause noticeable spikes/bumps in the mesh. Probably values less than 2 cause similar issues but they are probably too small for us to find when glancing through a large mesh.
To test this hypothesis, I added the following code to clamp averageRoot:
if(rCount && nonLinearFit) {averageRoot/=rCount;}
else {averageRoot=Real((x0-isoValue)/(x0-x1));
if (averageRoot > 1) averageRoot=1;
if (averageRoot < -1) averageRoot=-1;
}
This significantly reduced the size of the spikes to more manageable sizes, such that they would not even be noticeable if I didn't know where to look for them on the mesh. However, there is still a small, very artificial-looking "bump" on the mesh with extremely sharp angles. It is not smooth.
Unfortunately, neither my coworker or I know enough about these algorithms to correct this bug once and for all to get completely smooth meshes with no such sharp angles in the result. This calculation is apparently problematic...
This bug exists in both MeshLab and the original research code at http://www.cs.jhu.edu/~misha/Code/PoissonRecon/. We tested and found spikes in the results of both programs. Our testing was primarily done with the original research code (just because it was easier to recompile & test with a single-purpose console program).
Spike screenshot
Close up view of the start of the spike
Another close-up view of the start of the spike
Result when averageRoot clamped to (-1, 1)
Bug fix from Michael Kazhdan (one of the original authors)
I e-mailed one of the original authors of this algorithm (Michael Kazhdan) regarding this issue, and he suggested some changes which led to a small patch. The patch is attached as "MultiGridOctreeData Patch.txt" and should be applied to SVN revision 1608 of MultiGridOctreeData.inl. I tested it on several data sets that were previously affected by these spikes and there are no longer any spikes with this patch applied.