[vtkusers] TriangleArea() return the not correct value.
Cory Quammen
cory.quammen at kitware.com
Thu Aug 28 23:10:23 EDT 2014
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
More information about the vtkusers
mailing list