## 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)`