gvl_calc2.c

Go to the documentation of this file.

/*!
   \file lib/ogsf/gvl_calc2.c
 
   \brief OGSF library - loading and manipulating volumes, MarchingCubes 33
   Algorithm (lower level functions)
 
   GRASS OpenGL gsurf OGSF Library
 
   Based on implementation of MarchingCubes 33 Algorithm by
   Thomas Lewiner, thomas.lewiner@polytechnique.org, Math Dept, PUC-Rio
 
   (C) 1999-2008 by the GRASS Development Team
 
   This program is free software under the
   GNU General Public License (>=v2).
   Read the file COPYING that comes with GRASS
   for details.
 
   \author Tomas Paudits, 2004
   \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
 */
 
#include <float.h>
 
#include <grass/gis.h>
#include <grass/ogsf.h>
 
#include "mc33_table.h"
 
unsigned char m_case, m_config, m_subconfig;
 
/*!
   \brief ADD
 
   \param face
   \param v
 
   \return
 */
int mc33_test_face(char face, float *v)
{
    float A, B, C, D;
 
    switch (face) {
    case -1:
    case 1:
        A = v[0];
        B = v[4];
        C = v[5];
        D = v[1];
        break;
 
    case -2:
    case 2:
        A = v[1];
        B = v[5];
        C = v[6];
        D = v[2];
        break;
 
    case -3:
    case 3:
        A = v[2];
        B = v[6];
        C = v[7];
        D = v[3];
        break;
 
    case -4:
    case 4:
        A = v[3];
        B = v[7];
        C = v[4];
        D = v[0];
        break;
 
    case -5:
    case 5:
        A = v[0];
        B = v[3];
        C = v[2];
        D = v[1];
        break;
 
    case -6:
    case 6:
        A = v[4];
        B = v[7];
        C = v[6];
        D = v[5];
        break;
 
    default:
        fprintf(stderr, "Invalid face code %d\n", face);
        A = B = C = D = 0;
    };
 
    return face * A * (A * C - B * D) >= 0;
}
 
/*!
   \brief ADD
 
   \param s
   \param v
 
   \return
 */
int mc33_test_interior(char s, float *v)
{
    float t, At = 0, Bt = 0, Ct = 0, Dt = 0, a, b;
    char test = 0;
    char edge = -1;
 
    switch (m_case) {
    case 4:
    case 10:
        a = (v[4] - v[0]) * (v[6] - v[2]) - (v[7] - v[3]) * (v[5] - v[1]);
        b = v[2] * (v[4] - v[0]) + v[0] * (v[6] - v[2]) - v[1] * (v[7] - v[3]) -
            v[3] * (v[5] - v[1]);
        t = -b / (2 * a);
 
        if (t < 0 || t > 1)
            return s > 0;
 
        At = v[0] + (v[4] - v[0]) * t;
        Bt = v[3] + (v[7] - v[3]) * t;
        Ct = v[2] + (v[6] - v[2]) * t;
        Dt = v[1] + (v[5] - v[1]) * t;
        break;
 
    case 6:
    case 7:
    case 12:
    case 13:
        switch (m_case) {
        case 6:
            edge = cell_table[OFFSET_T6_1_1 + m_config].polys[0];
            break;
        case 7:
            edge = cell_table[OFFSET_T7_4_1 + m_config].polys[13];
            break;
        case 12:
            edge = cell_table[OFFSET_T12_2_S1 + m_config].polys[14];
            break;
        case 13:
            edge = cell_table[OFFSET_T13_5_1 + m_subconfig + m_config * 4]
                       .polys[2];
            break;
        }
 
        switch (edge) {
        case 0:
            t = v[0] / (v[0] - v[1]);
            At = 0;
            Bt = v[3] + (v[2] - v[3]) * t;
            Ct = v[7] + (v[6] - v[7]) * t;
            Dt = v[4] + (v[5] - v[4]) * t;
            break;
        case 1:
            t = v[1] / (v[1] - v[2]);
            At = 0;
            Bt = v[0] + (v[3] - v[0]) * t;
            Ct = v[4] + (v[7] - v[4]) * t;
            Dt = v[5] + (v[6] - v[5]) * t;
            break;
        case 2:
            t = v[2] / (v[2] - v[3]);
            At = 0;
            Bt = v[1] + (v[0] - v[1]) * t;
            Ct = v[5] + (v[4] - v[5]) * t;
            Dt = v[6] + (v[7] - v[6]) * t;
            break;
        case 3:
            t = v[3] / (v[3] - v[0]);
            At = 0;
            Bt = v[2] + (v[1] - v[2]) * t;
            Ct = v[6] + (v[5] - v[6]) * t;
            Dt = v[7] + (v[4] - v[7]) * t;
            break;
        case 4:
            t = v[4] / (v[4] - v[5]);
            At = 0;
            Bt = v[7] + (v[6] - v[7]) * t;
            Ct = v[3] + (v[2] - v[3]) * t;
            Dt = v[0] + (v[1] - v[0]) * t;
            break;
        case 5:
            t = v[5] / (v[5] - v[6]);
            At = 0;
            Bt = v[4] + (v[7] - v[4]) * t;
            Ct = v[0] + (v[3] - v[0]) * t;
            Dt = v[1] + (v[2] - v[1]) * t;
            break;
        case 6:
            t = v[6] / (v[6] - v[7]);
            At = 0;
            Bt = v[5] + (v[4] - v[5]) * t;
            Ct = v[1] + (v[0] - v[1]) * t;
            Dt = v[2] + (v[3] - v[2]) * t;
            break;
        case 7:
            t = v[7] / (v[7] - v[4]);
            At = 0;
            Bt = v[6] + (v[5] - v[6]) * t;
            Ct = v[2] + (v[1] - v[2]) * t;
            Dt = v[3] + (v[0] - v[3]) * t;
            break;
        case 8:
            t = v[0] / (v[0] - v[4]);
            At = 0;
            Bt = v[3] + (v[7] - v[3]) * t;
            Ct = v[2] + (v[6] - v[2]) * t;
            Dt = v[1] + (v[5] - v[1]) * t;
            break;
        case 9:
            t = v[1] / (v[1] - v[5]);
            At = 0;
            Bt = v[0] + (v[4] - v[0]) * t;
            Ct = v[3] + (v[7] - v[3]) * t;
            Dt = v[2] + (v[6] - v[2]) * t;
            break;
        case 10:
            t = v[2] / (v[2] - v[6]);
            At = 0;
            Bt = v[1] + (v[5] - v[1]) * t;
            Ct = v[0] + (v[4] - v[0]) * t;
            Dt = v[3] + (v[7] - v[3]) * t;
            break;
        case 11:
            t = v[3] / (v[3] - v[7]);
            At = 0;
            Bt = v[2] + (v[6] - v[2]) * t;
            Ct = v[1] + (v[5] - v[1]) * t;
            Dt = v[0] + (v[4] - v[0]) * t;
            break;
        default:
            fprintf(stderr, "Invalid edge %d\n", edge);
            break;
        }
        break;
 
    default:
        fprintf(stderr, "Invalid ambiguous case %d\n", m_case);
        break;
    }
 
    if (At >= 0)
        test++;
    if (Bt >= 0)
        test += 2;
    if (Ct >= 0)
        test += 4;
    if (Dt >= 0)
        test += 8;
 
    switch (test) {
    case 0:
        return s > 0;
    case 1:
        return s > 0;
    case 2:
        return s > 0;
    case 3:
        return s > 0;
    case 4:
        return s > 0;
    case 5:
        if (At * Ct < Bt * Dt)
            return s > 0;
        break;
    case 6:
        return s > 0;
    case 7:
        return s < 0;
    case 8:
        return s > 0;
    case 9:
        return s > 0;
    case 10:
        if (At * Ct >= Bt * Dt)
            return s > 0;
        break;
    case 11:
        return s < 0;
    case 12:
        return s > 0;
    case 13:
        return s < 0;
    case 14:
        return s < 0;
    case 15:
        return s < 0;
    }
 
    return s < 0;
}
 
/*!
   \brief ADD
 
   \param c_ndx
   \param v
 
   \return
 */
int mc33_process_cube(int c_ndx, float *v)
{
    m_case = cases[c_ndx][0];
    m_config = cases[c_ndx][1];
    m_subconfig = 0;
 
    switch (m_case) {
    case 0:
        return -1;
 
    case 1:
        return OFFSET_T1 + m_config;
 
    case 2:
        return OFFSET_T2 + m_config;
 
    case 3:
        if (mc33_test_face(test[OFFSET_TEST3 + m_config][0], v))
            return OFFSET_T3_2 + m_config; /* 3.2 */
        else
            return OFFSET_T3_1 + m_config; /* 3.1 */
 
    case 4:
        if (mc33_test_interior(test[OFFSET_TEST4 + m_config][0], v))
            return OFFSET_T4_1 + m_config; /* 4.1 */
        else
            return OFFSET_T4_2 + m_config; /* 4.2 */
 
    case 5:
        return OFFSET_T5 + m_config;
 
    case 6:
        if (mc33_test_face(test[OFFSET_TEST6 + m_config][0], v))
            return OFFSET_T6_2 + m_config; /* 6.2 */
        else {
            if (mc33_test_interior(test[OFFSET_TEST6 + m_config][1], v))
                return OFFSET_T6_1_1 + m_config; /* 6.1.1 */
            else
                return OFFSET_T6_1_2 + m_config; /* 6.1.2 */
        }
 
    case 7:
        if (mc33_test_face(test[OFFSET_TEST7 + m_config][0], v))
            m_subconfig += 1;
        if (mc33_test_face(test[OFFSET_TEST7 + m_config][1], v))
            m_subconfig += 2;
        if (mc33_test_face(test[OFFSET_TEST7 + m_config][2], v))
            m_subconfig += 4;
 
        switch (subconfig7[m_subconfig]) {
        case 0:
            if (mc33_test_interior(test[OFFSET_TEST7 + m_config][3], v))
                return OFFSET_T7_4_2 + m_config; /* 7.4.2 */
            else
                return OFFSET_T7_4_1 + m_config; /* 7.4.1 */
        case 1:
            return OFFSET_T7_3_S1 + m_config; /* 7.3 */
        case 2:
            return OFFSET_T7_3_S2 + m_config; /* 7.3 */
        case 3:
            return OFFSET_T7_3_S3 + m_config; /* 7.3 */
        case 4:
            return OFFSET_T7_2_S1 + m_config; /* 7.2 */
        case 5:
            return OFFSET_T7_2_S2 + m_config; /* 7.2 */
        case 6:
            return OFFSET_T7_2_S3 + m_config; /* 7.2 */
        case 7:
            return OFFSET_T7_1 + m_config; /* 7.1 */
        };
        break; /* will not reach this as previous switch is exhaustive */
 
    case 8:
        return OFFSET_T8 + m_config;
 
    case 9:
        return OFFSET_T9 + m_config;
 
    case 10:
        if (mc33_test_face(test[OFFSET_TEST10 + m_config][0], v)) {
            if (mc33_test_face(test[OFFSET_TEST10 + m_config][1], v))
                return OFFSET_T10_1_1_S2 + m_config; /* 10.1.1 */
            else {
                return OFFSET_T10_2_S1 + m_config; /* 10.2 */
            }
        }
        else {
            if (mc33_test_face(test[OFFSET_TEST10 + m_config][1], v)) {
                return OFFSET_T10_2_S2 + m_config; /* 10.2 */
            }
            else {
                if (mc33_test_interior(test[OFFSET_TEST10 + m_config][2], v))
                    return OFFSET_T10_1_1_S1 + m_config; /* 10.1.1 */
                else
                    return OFFSET_T10_1_2 + m_config; /* 10.1.2 */
            }
        }
 
    case 11:
        return OFFSET_T11 + m_config;
 
    case 12:
        if (mc33_test_face(test[OFFSET_TEST12 + m_config][0], v)) {
            if (mc33_test_face(test[OFFSET_TEST12 + m_config][1], v))
                return OFFSET_T12_1_1_S2 + m_config; /* 12.1.1 */
            else {
                return OFFSET_T12_2_S1 + m_config; /* 12.2 */
            }
        }
        else {
            if (mc33_test_face(test[OFFSET_TEST12 + m_config][1], v)) {
                return OFFSET_T12_2_S2 + m_config; /* 12.2 */
            }
            else {
                if (mc33_test_interior(test[OFFSET_TEST12 + m_config][2], v))
                    return OFFSET_T12_1_1_S1 + m_config; /* 12.1.1 */
                else
                    return OFFSET_T12_1_2 + m_config; /* 12.1.2 */
            }
        }
 
    case 13:
        if (mc33_test_face(test[OFFSET_TEST13 + m_config][0], v))
            m_subconfig += 1;
        if (mc33_test_face(test[OFFSET_TEST13 + m_config][1], v))
            m_subconfig += 2;
        if (mc33_test_face(test[OFFSET_TEST13 + m_config][2], v))
            m_subconfig += 4;
        if (mc33_test_face(test[OFFSET_TEST13 + m_config][3], v))
            m_subconfig += 8;
        if (mc33_test_face(test[OFFSET_TEST13 + m_config][4], v))
            m_subconfig += 16;
        if (mc33_test_face(test[OFFSET_TEST13 + m_config][5], v))
            m_subconfig += 32;
 
        switch (subconfig13[m_subconfig]) {
        case 0: /* 13.1 */
            return OFFSET_T13_1_S1 + m_config;
 
        case 1: /* 13.2 */
            return OFFSET_T13_2_S1 + 0 + m_config * 6;
        case 2: /* 13.2 */
            return OFFSET_T13_2_S1 + 1 + m_config * 6;
        case 3: /* 13.2 */
            return OFFSET_T13_2_S1 + 2 + m_config * 6;
        case 4: /* 13.2 */
            return OFFSET_T13_2_S1 + 3 + m_config * 6;
        case 5: /* 13.2 */
            return OFFSET_T13_2_S1 + 4 + m_config * 6;
        case 6: /* 13.2 */
            return OFFSET_T13_2_S1 + 5 + m_config * 6;
 
        case 7: /* 13.3 */
            return OFFSET_T13_3_S1 + 0 + m_config * 12;
        case 8: /* 13.3 */
            return OFFSET_T13_3_S1 + 1 + m_config * 12;
        case 9: /* 13.3 */
            return OFFSET_T13_3_S1 + 2 + m_config * 12;
        case 10: /* 13.3 */
            return OFFSET_T13_3_S1 + 3 + m_config * 12;
        case 11: /* 13.3 */
            return OFFSET_T13_3_S1 + 4 + m_config * 12;
        case 12: /* 13.3 */
            return OFFSET_T13_3_S1 + 5 + m_config * 12;
        case 13: /* 13.3 */
            return OFFSET_T13_3_S1 + 6 + m_config * 12;
        case 14: /* 13.3 */
            return OFFSET_T13_3_S1 + 7 + m_config * 12;
        case 15: /* 13.3 */
            return OFFSET_T13_3_S1 + 8 + m_config * 12;
        case 16: /* 13.3 */
            return OFFSET_T13_3_S1 + 9 + m_config * 12;
        case 17: /* 13.3 */
            return OFFSET_T13_3_S1 + 10 + m_config * 12;
        case 18: /* 13.3 */
            return OFFSET_T13_3_S1 + 11 + m_config * 12;
 
        case 19: /* 13.4 */
            return OFFSET_T13_4 + 0 + m_config * 4;
        case 20: /* 13.4 */
            return OFFSET_T13_4 + 1 + m_config * 4;
        case 21: /* 13.4 */
            return OFFSET_T13_4 + 2 + m_config * 4;
        case 22: /* 13.4 */
            return OFFSET_T13_4 + 3 + m_config * 4;
 
        case 23: /* 13.5 */
            m_subconfig = 0;
            if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
                return OFFSET_T13_5_1 + 0 + m_config * 4;
            else
                return OFFSET_T13_5_2 + 0 + m_config * 4;
        case 24: /* 13.5 */
            m_subconfig = 1;
            if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
                return OFFSET_T13_5_1 + 1 + m_config * 4;
            else
                return OFFSET_T13_5_2 + 1 + m_config * 4;
        case 25: /* 13.5 */
            m_subconfig = 2;
            if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
                return OFFSET_T13_5_1 + 2 + m_config * 4;
            else
                return OFFSET_T13_5_2 + 2 + m_config * 4;
        case 26: /* 13.5 */
            m_subconfig = 3;
            if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
                return OFFSET_T13_5_1 + 3 + m_config * 4;
            else
                return OFFSET_T13_5_2 + 3 + m_config * 4;
 
        case 27: /* 13.3 */
            return OFFSET_T13_3_S2 + 0 + m_config * 12;
        case 28: /* 13.3 */
            return OFFSET_T13_3_S2 + 1 + m_config * 12;
        case 29: /* 13.3 */
            return OFFSET_T13_3_S2 + 2 + m_config * 12;
        case 30: /* 13.3 */
            return OFFSET_T13_3_S2 + 3 + m_config * 12;
        case 31: /* 13.3 */
            return OFFSET_T13_3_S2 + 4 + m_config * 12;
        case 32: /* 13.3 */
            return OFFSET_T13_3_S2 + 5 + m_config * 12;
        case 33: /* 13.3 */
            return OFFSET_T13_3_S2 + 6 + m_config * 12;
        case 34: /* 13.3 */
            return OFFSET_T13_3_S2 + 7 + m_config * 12;
        case 35: /* 13.3 */
            return OFFSET_T13_3_S2 + 8 + m_config * 12;
        case 36: /* 13.3 */
            return OFFSET_T13_3_S2 + 9 + m_config * 12;
        case 37: /* 13.3 */
            return OFFSET_T13_3_S2 + 10 + m_config * 12;
        case 38: /* 13.3 */
            return OFFSET_T13_3_S2 + 11 + m_config * 12;
 
        case 39: /* 13.2 */
            return OFFSET_T13_2_S2 + 0 + m_config * 6;
        case 40: /* 13.2 */
            return OFFSET_T13_2_S2 + 1 + m_config * 6;
        case 41: /* 13.2 */
            return OFFSET_T13_2_S2 + 2 + m_config * 6;
        case 42: /* 13.2 */
            return OFFSET_T13_2_S2 + 3 + m_config * 6;
        case 43: /* 13.2 */
            return OFFSET_T13_2_S2 + 4 + m_config * 6;
        case 44: /* 13.2 */
            return OFFSET_T13_2_S2 + 5 + m_config * 6;
 
        case 45: /* 13.1 */
            return OFFSET_T13_1_S2 + m_config;
 
        default:
            fprintf(stderr, "Marching Cubes: Impossible case 13?\n");
        }
        break; /* will not reach this as previous switch is exhaustive */
 
    case 14:
        return OFFSET_T14 + m_config;
    }
 
    return -1;
}