[vtkusers] TriangleArea() return the not correct value.
Masato Asou
asou at soum.co.jp
Fri Aug 29 02:16:59 EDT 2014
> It looks like you are running into numerical precision problems
I understand.
> We should consider using this in VTK,
I hope so.
Thanks.
--
ASOU Masato
From: Cory Quammen <cory.quammen at kitware.com>
Date: Thu, 28 Aug 2014 23:10:23 -0400
> 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