[vtkusers] TriangleArea() return the not correct value.
David Gobbi
david.gobbi at gmail.com
Thu Aug 28 23:19:08 EDT 2014
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
More information about the vtkusers
mailing list