CODEWITHSUNDEEP
c++double
Cecilia
Cecilia

Reputation: 4721

Oh where has my precision gone with OpenMesh vector arithmetic?

Using doubles I would expect to have about 15 decimal points of precision. I know that many decimal numbers are not exactly representable in floating point notation, so I would get an approximation for 1/3 for example. However, using a double I would expect an approximation that was correct to about 15 decimal points. I would also expect to retain that level of accuracy when doing arithmetic.

However, in the following example, I try to calculate the area of a triangle using Heron's formula and OpenMesh::Vec3d which are backed by OpenMesh::VectorDataT<double,3> and end up with a result that is only accurate to 5 decimal points.

The correct result is area = 8.19922e-8, but I'm getting area=8.1992238711962083e-8. Any ideas where this is coming from?

The suggestion that this might result from the instability in Heron's Formula is a good one, but unfortunately is not the case in this example. I have added code which calculates the stable variation on Heron for those who might be interested. In this example, u.norm()>v.norm()>w.norm().

#include <OpenMesh/Core/Mesh/PolyMesh_ArrayKernelT.hh>

int main()
{
    //triangle vertices
    OpenMesh::Vec3d x(0.051051, 0.057411, 0.001355);
    OpenMesh::Vec3d y(0.050981, 0.057337, -0.000678);
    OpenMesh::Vec3d z(0.050949, 0.057303, 0.0);

    //edge vectors
    OpenMesh::Vec3d u = x-y;
    OpenMesh::Vec3d v = x-z;
    OpenMesh::Vec3d w = y-z;

    //Heron's Formula
    double semiP = (u.norm() + v.norm() + w.norm())/2.0; 
    double area = sqrt(semiP * (semiP - u.norm()) * (semiP - v.norm()) * (semiP - w.norm()) ); 

    //Heron's Formula for small angles
    double areaSmall = sqrt((u.norm() + (v.norm()+w.norm()))*(w.norm()-(u.norm()-v.norm()))*(w.norm()+(u.norm()-v.norm()))*(u.norm()+(v.norm()-w.norm())))/4.0;
}

Upvotes: 1

Views: 300

Answers (3)

tmyklebu
tmyklebu

Reputation: 14205

To 75 decimal places, the correct area of your triangle is

0.000000081992238711963087279421583920293974467992093148008322378721298327364.

If I replace the nine double constants you have with their decimal equivalents, I get

0.000000081992238711965902754749500279615357792172906541206211853522524016959

It would appear that you are not getting what you're expecting because you're expecting something unreasonable.

Upvotes: 3

Mark Ransom
Mark Ransom

Reputation: 308091

Any calculation involving subtraction will result in a loss of precision, if the values are at all close to each other. How many significant digits do you expect from this subtraction?

  1.23456789012345
- 1.23456789000000
  ----------------
  0.00000000012345

Both operands have 15 digits of precision, but the result only has 5.

Upvotes: 2

alain
alain

Reputation: 12037

Heron's formula is numerically unstable. If you have a very "flat" triangle with small angles, the sum of the two small sides is almost the long side, so one of the terms gets very small. If, for example, a and b are the small sides,

(s - c)

will be very small, because

s = (a + b + c)/2

is nearly equal to c.

The wikipedia article about herons formula mentions a stable alternative:

Arrange the sides such that a > b > c and use

A = 1/4*sqrt((a + (b + c))*(c - (a - b))*(c + (a - b))*(a + (b - c)))

Upvotes: 3

Related Questions