source: GTP/trunk/Lib/Geom/shared/GTGeometry/src/libs/vmi/src/area.cpp @ 983

#include <stdio.h>
#include <math.h>

#include "../include/area.h"

using namespace VMI;

/*
   m00 m01 m02
   m10 m11 m12
   m20 m21 m22
   
   |a1 a2 a3|
   |b1 b2 b3| = a1b2c3 - a1b3c2 - a2b1c3 + a2b3c1 + a3b1c2 - a3b2c1
   |c1 c2 c3|
*/
// Compute the determinant of a 3x3 matrix
double VMI::computeMatDet(double m[3][3]) {
    return ((m[0][0] * m[1][1] * m[2][2]) - (m[0][0] * m[1][2] * m[2][1]) - (m[0][1] * m[1][0] * m[2][2]) +
            (m[0][1] * m[1][2] * m[2][0]) + (m[0][2] * m[1][0] * m[2][1]) - (m[0][2] * m[1][1] * m[2][0]));
}
/*
   v00 x1
   v01 y1
   v02 z1

   v10 x2
   v11 y2
   v12 z2

   v20 x3
   v21 y3
   v22 z3
                  |y1 z1 1|       |z1 x1 1|       |x1 y1 1|
   A = 1/2 sqrt( (|y2 z2 1|)^2 + (|z2 x2 1|)^2 + (|x2 y2 1|)^2 )
                  |y3 z3 1|       |z3 x3 1|       |x3 y3 1|
*/
// Compute the area of a 3D triangle
double VMI::triangleArea(double v0[3], double v1[3], double v2[3]) {
    double m0[3][3], m1[3][3], m2[3][3],
           area, part0, part1, part2,
           d0, d1, d2;
    int i;
    
    for (i=0; i<3; i++) {
        
        m0[i][2] = 1.0;
        m1[i][2] = 1.0;
        m2[i][2] = 1.0;
    }
    // m0
    m0[0][0] = v0[1];
    m0[1][0] = v1[1];
    m0[2][0] = v2[1];

    m0[0][1] = v0[2];
    m0[1][1] = v1[2];
    m0[2][1] = v2[2];

    // m1
    m1[0][0] = v0[2];
    m1[1][0] = v1[2];
    m1[2][0] = v2[2];

    m1[0][1] = v0[0];
    m1[1][1] = v1[0];
    m1[2][1] = v2[0];

    // m2
    m2[0][0] = v0[0];
    m2[1][0] = v1[0];
    m2[2][0] = v2[0];

    m2[0][1] = v0[1];
    m2[1][1] = v1[1];
    m2[2][1] = v2[1];

    d0 = computeMatDet(m0);
    d1 = computeMatDet(m1);
    d2 = computeMatDet(m2);
    
    part0 = d0 * d0;
    part1 = d1 * d1;
    part2 = d2 * d2;
    
    area = sqrt(part0 + part1 + part2) / 2.0;
    
    return area;
}
/*
int main (int argc, char **argv) {
    double v0[] = {0.0, 0.0, 0.0},
           v1[] = {4.0, 4.0, 0.0},
           v2[] = {4.0, 0.0, 0.0},
           m[3][3] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
    
    printf("Det: %f\n", computeMatDet(m));
    
    printf("Area: %f\n", triangleArea(v0, v1, v2));
    
    getchar();
    return 0;
}
*/