ComputeOBB floating point arithmetic accumulates excessive error
Since I ran across the problem myself, and it appears to have not yet been submitted as an issue after 4 years nor corrected, I'm posting Oliver Weinheimer's discovery as noted in the VTK-Users mailing list on Aug 12, 2014: http://vtk.1045678.n5.nabble.com/Problems-with-ComputeOBB-wrong-OBB-computed-td5728174.html
From Oliver's second post:
The problem with
void vtkOBBTree::ComputeOBB(vtkPoints *pts, double corner[3], double max[3], double mid[3], double min[3], double size[3])
is a problem of the accuracy of floating point arithmetic. The objects I am dealing with contain several million points. This leads to problems with the mean value calculations in ComputeOBB. At two locations mean values are calculated in ComputeOBB, the second location introduced the greatest part of the errors in my case. I replaced the mean calculations with the method you can find in Knuth's Vol 2 of The Art of Computer Programming, 1998 edition, page 232.This solved at least my problems with the method!
I replaced
for (pointId=0; pointId < numPts; pointId++ )
{
pts->GetPoint(pointId, x);
for (i=0; i < 3; i++)
{
mean[i] += x[i];
}
}
for (i=0; i < 3; i++)
{
mean[i] /= numPts;
}
by these lines:
for (pointId = 0; pointId < numPts; pointId++)
{
pts->GetPoint(pointId, x);
for (i = 0; i < 3; i++){
mean[i] += (x[i] - mean[i])/(pointId+1);
}
}
and I replaced
for (pointId=0; pointId < numPts; pointId++ )
{
pts->GetPoint(pointId, x);
xp[0] = x[0] - mean[0]; xp[1] = x[1] - mean[1]; xp[2] = x[2] - mean[2];
for (i=0; i < 3; i++)
{
a0[i] += xp[0] * xp[i];
a1[i] += xp[1] * xp[i];
a2[i] += xp[2] * xp[i];
}
}//for all points
for (i=0; i < 3; i++)
{
a0[i] /= numPts;
a1[i] /= numPts;
a2[i] /= numPts;
}
by
for (pointId = 0; pointId < numPts; pointId++)
{
pts->GetPoint(pointId, x);
xp[0] = x[0] - mean[0];
xp[1] = x[1] - mean[1];
xp[2] = x[2] - mean[2];
for (i = 0; i < 3; i++)
{
a0[i] += (xp[0] * xp[i] - a0[i])/(pointId+1);
a1[i] += (xp[1] * xp[i] - a1[i])/(pointId+1);
a2[i] += (xp[2] * xp[i] - a2[i])/(pointId+1);
}
}
Later in Oliver's conversation with Goodwin Lawlor there was a suggestion to combine the operations into a single central function: static int vtkMath::GetMean(vtkDataArray *array, int dim, double* mean, double* variance)