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

#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <string.h>
#include <math.h>

#include "../include/mesh.h"
#include "../include/area.h"
#include "../include/global.h"

using namespace VMI;

#define MAX_NUM_TRI 450

void VMI::printItemList(int *list, int n) {
    int i;
    
    for (i=0; i<n; i++)
        printf("%d ",list[i]);

    printf("\n");
}

int VMI::findItem(int *list, int n, int item) {
    int i;

    for (i=0; i<n; i++)
        if (list[i] == item) return TRUE;
        return FALSE;
}

void VMI::addItem(int *list, int *n, int item) {

    if (!findItem(list, *n, item))
        list[(*n)++] = item;
}

int VMI::findEdge(Edge *e, GLuint num, GLuint _u, GLuint _v) {
    GLuint i, u, v;
    int found = -1; 

    for (i=0; i<num; i++) {

        if (e[i].enable == TRUE) {

            u = e[i].u;
            v = e[i].v;

            if (((u == _u) && (v == _v)) ||
                ((u == _v) && (v == _u))) {
                found = i;
                break;
            }
        }
    }

    return found;
}

GLdouble VMI::computeTriangleArea(Vertex *vertices, Triangle *t) {
    GLuint i;
    GLdouble v0[3], v1[3], v2[3], a;
    
    // Compute triangle area
    i = t->indices[0];
    v0[0] = vertices[i].x;
    v0[1] = vertices[i].y;
    v0[2] = vertices[i].z;

    i = t->indices[1];
    v1[0] = vertices[i].x;
    v1[1] = vertices[i].y;
    v1[2] = vertices[i].z;

    i = t->indices[2];
    v2[0] = vertices[i].x;
    v2[1] = vertices[i].y;
    v2[2] = vertices[i].z;

    a = triangleArea(v0, v1, v2);
    
    //printf("Triangle area: %f\n", a)
    return a;  
}

void VMI::computeTriangleNormal(Vertex *vertices, Triangle *t) {
    GLuint i;
    GLfloat x0, y0, z0,
            x1, y1, z1,
            x2, y2, z2;
    GLfloat x10, y10, z10,
            x12, y12, z12,
            cpx, cpy, cpz, r;

    i = t->indices[0];
    x0 = vertices[i].x;
    y0 = vertices[i].y;
    z0 = vertices[i].z;

    i = t->indices[1];
    x1 = vertices[i].x;
    y1 = vertices[i].y;
    z1 = vertices[i].z;

    i = t->indices[2];
    x2 = vertices[i].x;
    y2 = vertices[i].y;
    z2 = vertices[i].z;

    /* Compute edge vectors */
    x10 = x1 - x0;
    y10 = y1 - y0;
    z10 = z1 - z0;
    x12 = x1 - x2;
    y12 = y1 - y2;
    z12 = z1 - z2;

    /* Compute the cross product */
    cpx = (z10 * y12) - (y10 * z12);
    cpy = (x10 * z12) - (z10 * x12);
    cpz = (y10 * x12) - (x10 * y12);

    /* Normalize the result to get the unit-length facet normal */
    r = sqrt(cpx * cpx + cpy * cpy + cpz * cpz);

    t->normal[0] = cpx / r;
    t->normal[1] = cpy / r;
    t->normal[2] = cpz / r;
}

Mesh *VMI::initMesh(GLMmodel* pmodel) {

    GLuint i, j, v1, v2, v3, n, e, t;
    Mesh *mesh;

    mesh = (Mesh *)malloc (sizeof(Mesh));
    
    if (mesh == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }

    mesh->vertices = (Vertex *)malloc (sizeof(Vertex) * pmodel->numvertices);
    
    if (mesh->vertices == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }

    mesh->triangles = (Triangle *)malloc (sizeof(Triangle) * pmodel->numtriangles);
    
    if (mesh->triangles == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }
    
    printf("Adding vertices...");
    for (i=1; i<=pmodel->numvertices; i++) { // Vertices start at 1

        mesh->vertices[i - 1].x = pmodel->vertices[3 * i + 0];
        mesh->vertices[i - 1].y = pmodel->vertices[3 * i + 1];
        mesh->vertices[i - 1].z = pmodel->vertices[3 * i + 2];
        mesh->vertices[i - 1].numTriangles = 0;
        mesh->vertices[i - 1].triangles = NULL;
        mesh->vertices[i - 1].enable = GL_TRUE;
    }
    printf("Ok\n");
    mesh->numVertices = pmodel->numvertices;
    mesh->currentNumVertices = pmodel->numvertices;

    printf("Adding triangles...");
    for (i=0; i<pmodel->numtriangles; i++) {
        mesh->triangles[i].id = i;
        mesh->triangles[i].indices[0] = pmodel->triangles[i].vindices[0] - 1;
        mesh->triangles[i].indices[1] = pmodel->triangles[i].vindices[1] - 1;
        mesh->triangles[i].indices[2] = pmodel->triangles[i].vindices[2] - 1;
        mesh->triangles[i].area = computeTriangleArea(mesh->vertices, &mesh->triangles[i]);
        //printf("\n%d a: %f",i , mesh->triangles[i].area);
        computeTriangleNormal(mesh->vertices, &mesh->triangles[i]);

        mesh->triangles[i].saliency = 0.0;

        mesh->triangles[i].enable = GL_TRUE;
        
        for (j=0; j<3; j++) {
            // Adding triangle i adjacent to 3 vertices
            v1 = mesh->triangles[i].indices[j];
   
            // Reallocate memory for the new adjacent triangle
            mesh->vertices[v1].triangles = (GLuint *)realloc(mesh->vertices[v1].triangles, (mesh->vertices[v1].numTriangles + 1) * sizeof(GLuint));
            addItem((int *)mesh->vertices[v1].triangles, (int *)&mesh->vertices[v1].numTriangles, i);
        }
    }
    //printf("\n");
    printf("Ok\n");

    mesh->numTriangles = pmodel->numtriangles;
    mesh->currentNumTriangles = pmodel->numtriangles;

    mesh->edges = (Edge *)malloc (sizeof(Edge) * mesh->numTriangles * 3);  // E = 3 T / 2
    
    if (mesh->edges == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }
    
    // Init edges
    for (i=0; i<mesh->numTriangles * 3; i++) {
        
        mesh->edges[i].triangles    = NULL;
        mesh->edges[i].numTriangles = 0;
    }

    printf("Adding edges...");
    n = 0;

    for (i=0; i<mesh->numTriangles; i++) {
        t = 0;
        v1 = mesh->triangles[i].indices[0];
        v2 = mesh->triangles[i].indices[1];
        v3 = mesh->triangles[i].indices[2];
        
        if ((e = findEdge(mesh->edges, n, v1, v2)) == -1) {
            mesh->edges[n].u = v1;
            mesh->edges[n].v = v2;
            mesh->edges[n].enable = GL_TRUE;
            
            // Reallocate memory for the new adjacent triangle
            mesh->edges[n].triangles = (GLuint *)realloc(mesh->edges[n].triangles, (mesh->edges[n].numTriangles + 1) * sizeof(GLuint));
            // Adding triangle i adjacent to edge n
            addItem((int *)mesh->edges[n].triangles, (int *)&mesh->edges[n].numTriangles, i);
            //printf("n:%d i:%d\n", n, i);
            
            // Adding edge n adjacent to triangle i
            addItem((int *)mesh->triangles[i].edges, (int *)&t, n);
            n++;
        } else {
            // Reallocate memory for the new adjacent triangle
            mesh->edges[e].triangles = (GLuint *)realloc(mesh->edges[e].triangles, (mesh->edges[e].numTriangles + 1) * sizeof(GLuint));
            // Adding triangle i adjacent to edge e
            addItem((int *)mesh->edges[e].triangles, (int *)&mesh->edges[e].numTriangles, i);
            //printf("n:%d i:%d\n", e, i);
            
            // Adding edge e adjacent to triangle i
            addItem((int *)mesh->triangles[i].edges, (int *)&t, e);
        }
        if (( e = findEdge(mesh->edges, n, v2, v3)) == -1) {
            mesh->edges[n].u = v2;
            mesh->edges[n].v = v3;
            mesh->edges[n].enable = GL_TRUE;
            
            // Reallocate memory for the new adjacent triangle
            mesh->edges[n].triangles = (GLuint *)realloc(mesh->edges[n].triangles, (mesh->edges[n].numTriangles + 1) * sizeof(GLuint));
            // Adding triangle i adjacent to edge n
            addItem((int *)mesh->edges[n].triangles, (int *)&mesh->edges[n].numTriangles, i);
            //printf("n:%d i:%d\n", n, i);
            
            // Adding edge n adjacent to triangle i
            addItem((int *)mesh->triangles[i].edges, (int *)&t, n);
            n++;
        } else {
            // Reallocate memory for the new adjacent triangle
            mesh->edges[e].triangles = (GLuint *)realloc(mesh->edges[e].triangles, (mesh->edges[e].numTriangles + 1) * sizeof(GLuint));
            // Adding triangle i adjacent to edge e
            addItem((int *)mesh->edges[e].triangles, (int *)&mesh->edges[e].numTriangles, i);
            //printf("n:%d i:%d\n", e, i);
            
            // Adding edge e adjacent to triangle i
            addItem((int *)mesh->triangles[i].edges, (int *)&t, e);
        }
        if ((e = findEdge(mesh->edges, n, v3, v1)) == -1) {
            mesh->edges[n].u = v3;
            mesh->edges[n].v = v1;
            mesh->edges[n].enable = GL_TRUE;
            
            // Reallocate memory for the new adjacent triangle
            mesh->edges[n].triangles = (GLuint *)realloc(mesh->edges[n].triangles, (mesh->edges[n].numTriangles + 1) * sizeof(GLuint));
            // Adding triangle i adjacent to edge n
            addItem((int *)mesh->edges[n].triangles, (int *)&mesh->edges[n].numTriangles, i);
            //printf("n:%d i:%d\n", n, i);
            
            // Adding edge n adjacent to triangle i
            addItem((int *)mesh->triangles[i].edges, (int *)&t, n);
            n++;
        } else  {
            // Reallocate memory for the new adjacent triangle
            mesh->edges[e].triangles = (GLuint *)realloc(mesh->edges[e].triangles, (mesh->edges[e].numTriangles + 1) * sizeof(GLuint));
            // Adding triangle i adjacent to edge e
            addItem((int *)mesh->edges[e].triangles, (int *)&mesh->edges[e].numTriangles, i);
            //printf("n:%d i:%d\n", e, i);
            
            // Adding edge e adjacent to triangle i
            addItem((int *)mesh->triangles[i].edges, (int *)&t, e);
        }
    }
    printf("Ok\n");
    mesh->numEdges = n;

    return mesh;
}

void VMI::printMesh(Mesh *mesh) {
    GLuint i, j;

    printf("Vertices (%d)\n", mesh->numVertices);
    printf("------------------------\n");
    
    for (i=0; i<mesh->numVertices; i++) {
        printf("v%d (%f,%f,%f)\n", i, mesh->vertices[i].x, mesh->vertices[i].y,
                                     mesh->vertices[i].z);
        printf("  t(%d)[", mesh->vertices[i].numTriangles);
        if (mesh->vertices[i].triangles != NULL) {
            printf("%d", mesh->vertices[i].triangles[0]);
            for (j=1; j<mesh->vertices[i].numTriangles; j++) {
                printf(",%d", mesh->vertices[i].triangles[j]);
            }
        }
        printf("]\n");
    }

    printf("Edges (%d)\n", mesh->numEdges);
    printf("------------------------\n");
    for (i=0; i<mesh->numEdges; i++) {
        printf("e%d (%d,%d)\n", i, mesh->edges[i].u, mesh->edges[i].v);
        printf("  t(%d)[", mesh->edges[i].numTriangles);
        if (mesh->edges[i].triangles != NULL) {
            printf("%d", mesh->edges[i].triangles[0]);
            for (j=1; j<mesh->edges[i].numTriangles; j++) {
                printf(",%d", mesh->edges[i].triangles[j]);
            }
        }
        printf("]\n");
    }

    printf("Triangles (%d)\n", mesh->numTriangles);
    printf("------------------------\n");
    for (i=0; i<mesh->numTriangles; i++) {
        printf("t%d (%d,%d,%d)\n", i, mesh->triangles[i].indices[0], mesh->triangles[i].indices[1],
                                     mesh->triangles[i].indices[2]);
        printf("  e(3)[");
        printf("%d", mesh->triangles[i].edges[0]);
        for (j=1; j<3; j++) {
            printf(",%d", mesh->triangles[i].edges[j]);
        }
        printf("]\n");

        printf("  n(3)[");
        printf("%f", mesh->triangles[i].normal[0]);
        for (j=1; j<3; j++) {
            printf(",%f", mesh->triangles[i].normal[j]);
        }
        printf("]\n");
       
        printf("  area: %f\n", mesh->triangles[i].area);
    }
}

void VMI::deleteMesh(Mesh *mesh) {
        GLuint i;

        if (NULL != mesh) {
                for (i=0; i<mesh->numVertices; i++) {
                        if (mesh->vertices[i].triangles != NULL)
                                free(mesh->vertices[i].triangles);

                }
                if (mesh->vertices != NULL) free(mesh->vertices);
                mesh->numVertices = 0;

                for (i=0; i<mesh->numEdges; i++) {
                        if (mesh->edges[i].triangles != NULL)
                                free(mesh->edges[i].triangles);
                }
                if (mesh->edges != NULL) free(mesh->edges);
                mesh->numEdges = 0;

                if (mesh->triangles != NULL) free(mesh->triangles);
                mesh->numTriangles = 0;

                free(mesh);

                mesh = NULL;
        }
}

int *VMI::trianglesAdjToEdge(Mesh *mesh, int e, int *n) {
    int triangles[MAX_NUM_TRI];
    int i, n0, n1, n2, 
        u = mesh->edges[e].u, 
        v = mesh->edges[e].v;
    int *list, num = 0;

    for (i=0; i<(int)mesh->numTriangles; i++) {
        
        if (mesh->triangles[i].enable == TRUE) {
            
            n0 = mesh->triangles[i].indices[0];
            n1 = mesh->triangles[i].indices[1];
            n2 = mesh->triangles[i].indices[2];
            
            if ((n0 == v) || (n1 == v) || (n2 == v) ||
                (n0 == u) || (n1 == u) || (n2 == u))
                
                triangles[num++] = i;
        }
    }

    list = (int*) malloc (sizeof(int) * num);
    if (NULL == list) {
        fprintf (stderr, "no more memory for list");
        exit(1);
    }

    memcpy(list, triangles, sizeof(int) * num);

    *n = num;

    return list;
}

int *VMI::trianglesAdjToVertex(Mesh *mesh, int v, int *n) {
    int triangles[MAX_NUM_TRI];
    int i, n0, n1, n2;
    int *list, num = 0;

    for (i=0; i<(int)mesh->numTriangles; i++) {
        
        if (mesh->triangles[i].enable == TRUE) {
            
            n0 = mesh->triangles[i].indices[0];
            n1 = mesh->triangles[i].indices[1];
            n2 = mesh->triangles[i].indices[2];
            
            if ((n0 == v) || (n1 == v) || (n2 == v))
                
                triangles[num++] = i;
        }
    }

    list = (int*) malloc (sizeof(int) * num);
    if (NULL == list) {
        fprintf (stderr, "no more memory for list");
        exit(1);
    }

    memcpy(list, triangles, sizeof(int) * num);

    *n = num;

    return list;
}

int *VMI::verticesAdjToVertex(Mesh *mesh, int v, int *n) {
    int numTr;
    int *triangles = trianglesAdjToVertex(mesh, v, &numTr);
    int *vertices = (int*) malloc (sizeof(int) * numTr * 2);
    int i, t, n0, n1, n2;
    int num = 0;

    for (i=0; i<numTr; i++) {
        t = triangles[i];

        n0 = mesh->triangles[t].indices[0];
        n1 = mesh->triangles[t].indices[1];
        n2 = mesh->triangles[t].indices[2];

        if (n0 != v) addItem(vertices, &num, n0);
        if (n1 != v) addItem(vertices, &num, n1);
        if (n2 != v) addItem(vertices, &num, n2);
    }

    *n = num;

    free(triangles);

    return vertices;
}

int *VMI::edgesAdjToVertex(Mesh *mesh, int v, int *n) {
    int numTr;
    int *triangles = trianglesAdjToVertex(mesh, v, &numTr);
    int *edges = (int*) malloc (sizeof(int) * numTr * 2);
    int i, t, n0, n1, n2, e;
    int num = 0;

    for (i=0; i<numTr; i++) {
        t = triangles[i];

        n0 = mesh->triangles[t].indices[0];
        n1 = mesh->triangles[t].indices[1];
        n2 = mesh->triangles[t].indices[2];
        
        if (((n0 == v) || (n1 == v)) &&
            ((e = findEdge(mesh->edges, mesh->numEdges, n0, n1)) != -1)) addItem(edges, &num, e);
 
        if (((n1 == v) || (n2 == v)) &&
            ((e = findEdge(mesh->edges, mesh->numEdges, n1, n2)) != -1)) addItem(edges, &num, e);
        
        if (((n0 == v) || (n2 == v)) &&
            ((e = findEdge(mesh->edges, mesh->numEdges, n0, n2)) != -1)) addItem(edges, &num, e);  
    }

    *n = num;

    free(triangles);

    return edges;
}

int *VMI::edgesAdjToVertices(Mesh *mesh, int *vertices, int numVertices, int *n) {
    int edges[MAX_NUM_TRI*2];
    int i, j, v, num = 0, e;
    int *list, numEdges;

    for (i=0; i<numVertices; i++) {
        v = vertices[i];

        list = edgesAdjToVertex(mesh, v, &numEdges);

        for (j=0; j<numEdges; j++) {
            e = list[j];
            
            addItem(edges, &num, e);
        }

        free(list);
    }

    list = (int*) malloc (sizeof(int) * num);
    if (NULL == list) {
        fprintf (stderr, "no more memory for list");
        exit(1);
    }

    memcpy(list, edges, sizeof(int) * num);

    *n = num;

    return list;
}

GLboolean VMI::findVertex(GLfloat *vertices, GLuint num, GLfloat x, GLfloat y, GLfloat z, int *pos) {
    GLuint i;
    GLfloat _x, _y, _z;
    GLboolean found = GL_FALSE; 
    
    *pos = -1;

    for (i=1; i<num; i++) { // i=0

        _x = vertices[3 * i + 0];
        _y = vertices[3 * i + 1];
        _z = vertices[3 * i + 2];

        if (((_x == x) && (_y == y) && (_z == z))) {
            found = GL_TRUE;
            *pos = (int)i;
            break;
        }
    }

    return found;
}

void VMI::updateModel(GLMmodel* pmodel, Mesh *mesh, GLuint numVertices, GLuint numTriangles) {
    GLuint i , v1, v2 ,v3, numV = 1, numT = 0;
    int pos;
    GLfloat x, y, z;
    
    if (pmodel->vertices != NULL) free(pmodel->vertices);

    pmodel->vertices = (GLfloat *)calloc ((numVertices + 1) * 3,  sizeof(GLfloat));
    

    if (pmodel->triangles != NULL) free(pmodel->triangles);

    pmodel->triangles = (GLMtriangle *)malloc (numTriangles * sizeof(GLMtriangle));

    
    for (i=0; i<mesh->numTriangles; i++) {
        if (mesh->triangles[i].enable) {
            v1 = mesh->triangles[i].indices[0];
            v2 = mesh->triangles[i].indices[1];
            v3 = mesh->triangles[i].indices[2];
            
            x = mesh->vertices[v1].x;
            y = mesh->vertices[v1].y;
            z = mesh->vertices[v1].z;
            
            if (!findVertex(pmodel->vertices, numV, x, y ,z, &pos)) {
                pmodel->triangles[numT].vindices[0] = numV;
                
                pmodel->vertices[3 * numV + 0] = x;
                pmodel->vertices[3 * numV + 1] = y;
                pmodel->vertices[3 * numV + 2] = z;
                numV++;
            } else pmodel->triangles[numT].vindices[0] = (GLuint)pos;
            
            x = mesh->vertices[v2].x;
            y = mesh->vertices[v2].y;
            z = mesh->vertices[v2].z;
            
            if (!findVertex(pmodel->vertices, numV, x, y ,z, &pos)) {
                pmodel->triangles[numT].vindices[1] = numV;
                
                pmodel->vertices[3 * numV + 0] = x;
                pmodel->vertices[3 * numV + 1] = y;
                pmodel->vertices[3 * numV + 2] = z;
                numV++;
            } else pmodel->triangles[numT].vindices[1] = (GLuint)pos;
            
            x = mesh->vertices[v3].x;
            y = mesh->vertices[v3].y;
            z = mesh->vertices[v3].z;
            
            if (!findVertex(pmodel->vertices, numV, x, y ,z, &pos)) {
                pmodel->triangles[numT].vindices[2] = numV;
                
                pmodel->vertices[3 * numV + 0] = x;
                pmodel->vertices[3 * numV + 1] = y;
                pmodel->vertices[3 * numV + 2] = z;
                numV++;
            } else pmodel->triangles[numT].vindices[2] = (GLuint)pos;
            numT++;
        }
    }
    pmodel->numvertices  = numV - 1;
    pmodel->numtriangles = numT;
}

void VMI::saveModel(char *filename, GLMmodel* pmodel, Mesh *mesh) {
    GLuint i;

    // The output model contains only one group
    pmodel->numgroups = 1;

    if (pmodel->groups->name != NULL) free(pmodel->groups->name);
    pmodel->groups->name = strdup("default");
    pmodel->groups->material = 0;
    pmodel->groups->next = NULL;
    pmodel->groups->numtriangles = mesh->currentNumTriangles;

    if (pmodel->groups->triangles != NULL) free(pmodel->groups->triangles);
    pmodel->groups->triangles = (GLuint*)malloc(sizeof(GLuint) * mesh->currentNumTriangles);

    for (i=0; i<mesh->currentNumTriangles; i++)
        pmodel->groups->triangles[i] = i;

    updateModel(pmodel, mesh, mesh->currentNumTriangles * 3, mesh->currentNumTriangles);

    printf("glmWriteOBJ...");
    
    glmWriteOBJ(pmodel, filename, GLM_NONE);
    printf("Ok\n");
}