[vtkusers] TriangleArea() return the not correct value.
John Platt
jcplatt at dsl.pipex.com
Mon Sep 1 04:53:16 EDT 2014
Hi,
It could be worth looking at "Miscalculating Area and Angles of a
Needle-like Triangle" W. Kahan (24 March 2000)
http://www.cs.berkeley.edu/~wkahan/Triangle.pdf for implementation of
Heron's formula.
HTH
John.
----- Original Message -----
From: "Masato Asou" <asou at soum.co.jp>
To: <david.gobbi at gmail.com>
Cc: <vtkusers at vtk.org>
Sent: Friday, August 29, 2014 7:18 AM
Subject: Re: [vtkusers] TriangleArea() return the not correct value.
>I didn't aware about efficiency.
>
> I hope that I can choose whichever efficiency or quality.
> Thanks.
> --
> ASOU Masato
>
>> Note that there is a difference in efficiency, because Heron's formula
>> requires applying sqrt() four times, while VTK's formula requires
>> applying
>> sqrt() only once.
>
>> It might be worthwhile to use both Heron's formula and the current VTK
>> formula, with a discriminator to choose which based on the quality of the
>> triangle.
>> - David
>>
>>
>> On Thu, Aug 28, 2014 at 9:10 PM, Cory Quammen <cory.quammen at kitware.com>
>> wrote:
>>> Masato,
>>>
>>> It looks like you are running into numerical precision problems
>>> because two of your points are extremely close, almost close enough to
>>> be indistinguishable by double precision.
>>>
>>> The inconsistencies that depend on the order the points are given are
>>> caused by which sides are considered in the cross product and area
>>> calculation. Heron's formula, which I wasn't aware of, seems to
>>> eliminate this inconsistency. We should consider using this in VTK, in
>>> particular the numerically more stable implementation:
>>> http://en.wikipedia.org/wiki/Heron's_formula#Numerical_stability
>>>
>>> I don't know if my answer helps, but thank you for pointing out Heron's
>>> formula.
>>>
>>> Thanks,
>>> Cory
>>>
>>> On Thu, Aug 28, 2014 at 10:41 PM, Masato Asou <asou at soum.co.jp> wrote:
>>>> Hi,
>>>>
>>>> I have a probrem about vtkTriangle::TriangleArea(). I use VTK-5.10.1.
>>>>
>>>> The triangle to be calculated is what one side is very short.
>>>>
>>>> I calculate the area of the triangle in vtkTriangle::TriangleArea(),
>>>> the result are different to that using the Heron's formula and cross
>>>> products. I also calculate tha area is regarded as a rectanglar
>>>> triangle,
>>>> this value is also different vtkTriangle::TriangleArea()
>>>>
>>>> In addition, vtkTriangle::TriangleArea() return a different value when
>>>> change the order of the vertices.
>>>>
>>>> This problem be resolved?
>>>>
>>>> My system as:
>>>>
>>>> $ cat /etc/redhat-release
>>>> CentOS release 5.6 (Final)
>>>> $ uname -m -r -s -v
>>>> Linux 2.6.18-238.19.1.el5 #1 SMP Fri Jul 15 07:31:24 EDT 2011 x86_64
>>>>
>>>> My program as below:
>>>>
>>>> $ cat triangle.cpp
>>>> #include <math.h>
>>>> #include <ostream>
>>>> using namespace std;
>>>> #include <vtkTriangle.h>
>>>>
>>>> double
>>>> lineLength(const double p0[3], const double p1[3])
>>>> {
>>>> double len2 = 0.0;
>>>>
>>>> for(int i = 0; i < 3; i++) {
>>>> len2 += (p0[i] - p1[i]) * (p0[i] - p1[i]);
>>>> }
>>>> return sqrt(len2);
>>>> }
>>>>
>>>> double
>>>> CrossProduct(double p1[3], double p2[3], double p3[3])
>>>> {
>>>> double a[3], b[3], o[3], v[3];
>>>>
>>>> // Create the vector a and b.
>>>> a[0] = p2[0] - p1[0];
>>>> a[1] = p2[1] - p1[1];
>>>> a[2] = p2[2] - p1[2];
>>>> b[0] = p3[0] - p1[0];
>>>> b[1] = p3[1] - p1[1];
>>>> b[2] = p3[2] - p1[2];
>>>>
>>>> // Calculate the cross product of vector a and b.
>>>> o[0] = 0.0;
>>>> o[1] = 0.0;
>>>> o[2] = 0.0;
>>>> v[0] = (a[1] * b[2]) - (a[2] * b[1]);
>>>> v[1] = (a[2] * b[0]) - (a[0] * b[2]);
>>>> v[2] = (a[0] * b[1]) - (a[1] * b[0]);
>>>> return fabs(lineLength(o, v)) / 2;
>>>> }
>>>>
>>>> double
>>>> HeronsFormula(double p1[3], double p2[3], double p3[3])
>>>> {
>>>> double a = lineLength(p1, p2);
>>>> double b = lineLength(p2, p3);
>>>> double c = lineLength(p3, p1);
>>>>
>>>> return sqrt((a+b+c) * (-a+b+c) * (a-b+c) * (a+b-c)) * 0.25;
>>>> }
>>>>
>>>> int
>>>> main(void)
>>>> {
>>>> double A[3] = {-3.58117425311239667707, -3.58117425311239667707,
>>>> 1.0};
>>>> double B[3] = {-3.04058712655619833853, -3.04058712655619833853,
>>>> 0.95};
>>>> double C[3] = {-3.04058712655619833853, -3.04058712655619789444,
>>>> 0.95};
>>>> double a = lineLength(A, B);
>>>> double b = lineLength(B, C);
>>>> double c = lineLength(C, A);
>>>>
>>>> cout << fixed << setprecision(20);
>>>> cout << "TriangleArea(A, B, C): " << vtkTriangle::TriangleArea(A,
>>>> B, C) << endl;
>>>> cout << "TriangleArea(B, C, A): " << vtkTriangle::TriangleArea(B,
>>>> C, A) << endl;
>>>> cout << "CrossProduct(A, B, C): " << CrossProduct(A, B, C) <<
>>>> endl;
>>>> cout << "CrossProduct(B, C, A): " << CrossProduct(B, C, A) <<
>>>> endl;
>>>> cout << "HeronsFormula(A, B, C): " << HeronsFormula(A, B, C) <<
>>>> endl;
>>>> cout << "HeronsFormula(B, C, A): " << HeronsFormula(B, C, A) <<
>>>> endl;
>>>> cout << "RectangularTriangle: " << (a * b / 2) << endl;
>>>> cout << "line length: a: " << a << ", b: " << b << ", c: " << c <<
>>>> endl;
>>>>
>>>> return 0;
>>>> }
>>>> $ g++ -g -I../../src -I/home/asou/freesw/vtk-5.10.1/include/vtk-5.10
>>>> triangle.cpp -Xlinker -rpath -Xlinker
>>>> ../../src -L/home/asou/freesw/vtk-5.10.1/lib/vtk-5.10 -lvtkCommon -lvtkWidgets
>>>> -Xlinker -rpath -Xlinker /home/asou/freesw/vtk-5.10.1/lib/vtk-5.10
>>>> In file included from
>>>> /usr/lib/gcc/x86_64-redhat-linux/4.1.2/../../../../include/c++/4.1.2/backward/strstream:51,
>>>> from
>>>> /home/asou/freesw/vtk-5.10.1/include/vtk-5.10/vtkIOStream.h:108,
>>>> from
>>>> /home/asou/freesw/vtk-5.10.1/include/vtk-5.10/vtkSystemIncludes.h:40,
>>>> from
>>>> /home/asou/freesw/vtk-5.10.1/include/vtk-5.10/vtkIndent.h:24,
>>>> from
>>>> /home/asou/freesw/vtk-5.10.1/include/vtk-5.10/vtkObjectBase.h:43,
>>>> from
>>>> /home/asou/freesw/vtk-5.10.1/include/vtk-5.10/vtkObject.h:41,
>>>> from
>>>> /home/asou/freesw/vtk-5.10.1/include/vtk-5.10/vtkCell.h:40,
>>>> from
>>>> /home/asou/freesw/vtk-5.10.1/include/vtk-5.10/vtkTriangle.h:23,
>>>> from triangle.cpp:4:
>>>> /usr/lib/gcc/x86_64-redhat-linux/4.1.2/../../../../include/c++/4.1.2/backward/backward_warning.h:32:2:
>>>> warning: #warning This file includes at least one deprecated or
>>>> antiquated header. Please consider using one of the 32 headers found in
>>>> section 17.4.1.2 of the C++ standard. Examples include substituting the
>>>> <X> header for the <X.h> header for C++ includes, or <iostream> instead
>>>> of the deprecated header <iostream.h>. To disable this warning
>>>> use -Wno-deprecated.
>>>> $ ./a.out
>>>> TriangleArea(A, B, C): 0.00000000526835606386
>>>> TriangleArea(B, C, A): 0.00000000000000009839
>>>> CrossProduct(A, B, C): 0.00000000000000011168
>>>> CrossProduct(B, C, A): 0.00000000000000012055
>>>> HeronsFormula(A, B, C): 0.00000000000000011252
>>>> HeronsFormula(B, C, A): 0.00000000000000011252
>>>> RectangularTriangle: 0.00000000000000017012
>>>> line length: a: 0.76613894483740641039, b: 0.00000000000000044409, c:
>>>> 0.76613894483740674346
>>>> $
>>>>
>>>> Thanks in advance any help or suggestions.
>>>> --
>>>> ASOU Masato
>>>> _______________________________________________
>>>> Powered by www.kitware.com
>>>>
>>>> Visit other Kitware open-source projects at
>>>> http://www.kitware.com/opensource/opensource.html
>>>>
>>>> Please keep messages on-topic and check the VTK FAQ at:
>>>> http://www.vtk.org/Wiki/VTK_FAQ
>>>>
>>>> Follow this link to subscribe/unsubscribe:
>>>> http://public.kitware.com/mailman/listinfo/vtkusers
>>> _______________________________________________
>>> Powered by www.kitware.com
>>>
>>> Visit other Kitware open-source projects at
>>> http://www.kitware.com/opensource/opensource.html
>>>
>>> Please keep messages on-topic and check the VTK FAQ at:
>>> http://www.vtk.org/Wiki/VTK_FAQ
>>>
>>> Follow this link to subscribe/unsubscribe:
>>> http://public.kitware.com/mailman/listinfo/vtkusers
>>
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at:
> http://www.vtk.org/Wiki/VTK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/vtkusers
>
More information about the vtkusers
mailing list