<div dir="ltr">Hi Oliver,<div><br></div><div>That's a nice method - the mean is gradually accumulated, preserving the accuracy of the floating point value. I'm sure there are other places in VTK that would benefit from this code change...</div>
<div><br></div><div>If you're familiar with Gerrit, could you submit this change? I'll have a poke around the code base to see if there are other similar mean calculations.</div><div><br></div><div>Goodwin</div></div>
<div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, Aug 12, 2014 at 1:54 PM, Oliver Weinheimer <span dir="ltr"><<a href="mailto:mail@oliwe.com" target="_blank">mail@oliwe.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi all,<br>
<br>
the problem with<br>
<br>
void vtkOBBTree::ComputeOBB(vtkPoints *pts, double corner[3], double max[3],<br>
double mid[3], double min[3], double size[3])<br>
<br>
is a problem of the accuracy of floating point arithmetic.<br>
The objects I am dealing with contain several million points.<br>
This leads to problems with the mean value calculations in ComputeOBB.<br>
At two locations mean values are calculated in ComputeOBB, the second<br>
location introduced the greatest part of the errors in my case.<br>
I replaced the mean calculations with the method you can find in Knuth's Vol<br>
2 of The Art of Computer Programming, 1998 edition, page 232.<br>
<br>
This solved at least my problems with the method!<br>
<br>
(1)<br>
I replaced<br>
<br>
--- source code begin<br>
        for (pointId=0; pointId < numPts; pointId++ )<br>
        {<br>
                pts->GetPoint(pointId, x);<br>
                for (i=0; i < 3; i++)<br>
                {<br>
                        mean[i] += x[i];<br>
                }<br>
        }<br>
        for (i=0; i < 3; i++)<br>
        {<br>
                mean[i] /= numPts;<br>
        }<br>
--- source code end<br>
<br>
by these lines:<br>
<br>
--- source code begin<br>
for (pointId = 0; pointId < numPts; pointId++)<br>
{<br>
        pts->GetPoint(pointId, x);<br>
        for (i = 0; i < 3; i++){<br>
                mean[i] += (x[i] - mean[i])/(pointId+1);<br>
        }<br>
}<br>
<br>
--- source code end<br>
<br>
(2)<br>
<br>
and I replaced<br>
<br>
--- source code begin<br>
for (pointId=0; pointId < numPts; pointId++ )<br>
{<br>
        pts->GetPoint(pointId, x);<br>
        xp[0] = x[0] - mean[0]; xp[1] = x[1] - mean[1]; xp[2] = x[2] - mean[2];<br>
        for (i=0; i < 3; i++)<br>
        {<br>
                a0[i] += xp[0] * xp[i];<br>
                a1[i] += xp[1] * xp[i];<br>
                a2[i] += xp[2] * xp[i];<br>
        }<br>
}//for all points<br>
<br>
for (i=0; i < 3; i++)<br>
{<br>
        a0[i] /= numPts;<br>
        a1[i] /= numPts;<br>
        a2[i] /= numPts;<br>
}<br>
--- source code end<br>
<br>
by<br>
<br>
--- source code begin<br>
for (pointId = 0; pointId < numPts; pointId++)<br>
{<br>
        pts->GetPoint(pointId, x);<br>
        xp[0] = x[0] - mean[0];<br>
        xp[1] = x[1] - mean[1];<br>
        xp[2] = x[2] - mean[2];<br>
        for (i = 0; i < 3; i++)<br>
        {<br>
                a0[i] += (xp[0] * xp[i] - a0[i])/(pointId+1);<br>
                a1[i] += (xp[1] * xp[i] - a1[i])/(pointId+1);<br>
                a2[i] += (xp[2] * xp[i] - a2[i])/(pointId+1);<br>
        }<br>
}<br>
--- source code end<br>
<br>
Best Wishes,<br>
Oliver<br>
<br>
<br>
<br>
--<br>
View this message in context: <a href="http://vtk.1045678.n5.nabble.com/Problems-with-ComputeOBB-wrong-OBB-computed-tp5728174p5728187.html" target="_blank">http://vtk.1045678.n5.nabble.com/Problems-with-ComputeOBB-wrong-OBB-computed-tp5728174p5728187.html</a><br>

<div class="HOEnZb"><div class="h5">Sent from the VTK - Users mailing list archive at Nabble.com.<br>
_______________________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Please keep messages on-topic and check the VTK FAQ at: <a href="http://www.vtk.org/Wiki/VTK_FAQ" target="_blank">http://www.vtk.org/Wiki/VTK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://public.kitware.com/mailman/listinfo/vtkusers" target="_blank">http://public.kitware.com/mailman/listinfo/vtkusers</a><br>
</div></div></blockquote></div><br></div>