source: GTP/trunk/Lib/Geom/shared/GTGeometry/src/libs/vmi/src/simplify.cpp @ 1007

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

#include "../include/metrics.h"
#include "../include/global.h"
#include "../include/simplify.h"
#include "../include/saliency.h"

//#define RECOMPUTE_THE_WHOLE_HEAP
//#define MAKE_FULL_COPY

using namespace VMI;

///////////////////////////////////////////////////////////////////////////////

void VMI::bh_mydump(Mesh *mesh, bheap_t *h) {
    int i, d;

    printf("Heap status.\n");
    for(i = 1; i <= h->n; i++) {
        d = h->a[i].item;
        printf ("Heap [%d] = pri %f, e%d(%d, %d)\n", i,
                  h->a[i].key,d,
                  mesh->edges[d].u,
                  mesh->edges[d].v);
    }

    printf("\n");
    
    for(i = 2; i <= h->n; i++) {
        if(h->a[i].key < h->a[i/2].key) {
            printf("key error at entry %d, value %f\n", i, h->a[i].key);
            exit(1);
        }
    }
    for(i = 1; i <= h->n; i++) {
        if(h->p[h->a[i].item] != i) {
            printf("indexing error at entry %d", i); exit(1);
        }
    }    
}

bheap_t *VMI::initHeap(Mesh *mesh) {
    GLuint i;
    double cost;
    bheap_t *h = NULL;
#ifdef MAKE_FULL_COPY
    Triangle *triangleCopy = NULL;
    GLuint lastNumTriangles = mesh->numTriangles;

    // Allocate memory
    triangleCopy = (Triangle *)malloc(lastNumTriangles * sizeof(Triangle));  
    if (triangleCopy == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }
    // Make a copy of current triangle list
    memcpy(triangleCopy, mesh->triangles, lastNumTriangles * sizeof(Triangle));
#endif
    
    printf("Creating a heap for simplification...");
    
    h = bh_alloc(mesh->numEdges);

                //      Init percent update count.
                float   percent =       60.0 / mesh->numEdges;

    for (i = 0; i <mesh->numEdges; i++) {

                                //      Update progress bar.
                                mUPB(percent);

        if (mesh->edges[i].enable == TRUE) {

            cost = computeEdgeCost(mesh, i);
        
            // Add only valid edges
            if (cost != FLT_MAX)
                bh_insert(h, i, cost);

#ifdef MAKE_FULL_COPY
            // Restore initial triangle list
            memcpy(mesh->triangles, triangleCopy, lastNumTriangles * sizeof(Triangle));
            mesh->numTriangles = lastNumTriangles;
#endif
        }
    }
    
#ifdef MAKE_FULL_COPY
    // Free memory
    if (triangleCopy != NULL) free(triangleCopy);
#endif
    printf("Ok\n");

#ifdef DEBUG 
    bh_mydump(mesh, h);
    //getchar();
#endif

    return h;
}

bheap_t *VMI::updateHeap(bheap_t *h, Mesh *mesh, Change *c) {
    int e, i, n0, n1, n2, u, v;
#ifdef RECOMPUTE_THE_WHOLE_HEAP
    bheap_t *nh;
#else
    int *lv, *le, numV, numE;
#endif
    double cost;
#ifdef MAKE_FULL_COPY
    GLuint lastNumTriangles = mesh->numTriangles;
    Triangle *triangleCopy = NULL;
    
    // Allocate memory
    triangleCopy = (Triangle *)malloc(lastNumTriangles * sizeof(Triangle));  
    if (triangleCopy == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }
    // Make a copy of current triangle list
    memcpy(triangleCopy, mesh->triangles, lastNumTriangles * sizeof(Triangle));
#endif

    for (i=0; i<c->numDel; i++) {

        n0 = c->deleted[i].indices[0];
        n1 = c->deleted[i].indices[1];
        n2 = c->deleted[i].indices[2];

        if (((n0 == c->u) || (n1 == c->u)) &&
            ((e = findEdge(mesh->edges, mesh->numEdges, n0, n1)) != -1)) {

            mesh->edges[e].enable = FALSE;
            bh_delete(h, e);
        }

        if (((n1 == c->u) || (n2 == c->u)) &&
            ((e = findEdge(mesh->edges, mesh->numEdges, n1, n2)) != -1)) {
             
            mesh->edges[e].enable = FALSE;
            bh_delete(h, e);
        }

        if (((n0 == c->u) || (n2 == c->u)) &&
            ((e = findEdge(mesh->edges, mesh->numEdges, n0, n2)) != -1)) {
              
            mesh->edges[e].enable = FALSE;
            bh_delete(h, e);
        }
    }

    // Update mesh and heap
    for(i = 1; i <= h->n; i++) {
        e = h->a[i].item;

        if ((int)mesh->edges[e].u == c->u)
            mesh->edges[e].u = c->v;

        if ((int)mesh->edges[e].v == c->u)
            mesh->edges[e].v = c->v;
        
        // Check edge
        u = mesh->edges[e].u;
        v = mesh->edges[e].v;
    
        // if it is not a valid edge we simply delete it
        if ((u == v) || 
            (u == c->u) || (v == c->u)) {

            mesh->edges[e].enable = FALSE;
            bh_delete(h, e);
        }
    }

#ifdef RECOMPUTE_THE_WHOLE_HEAP

    nh = bh_alloc(mesh->numEdges);

    // Recompute heap
    for(i = 1; i <= h->n; i++) {
        e = h->a[i].item;
        
        cost = computeEdgeCost(mesh, e);
        
        // Add only valid edges
        if (cost != FLT_MAX)
            bh_insert(nh, e, cost);
        
#ifdef MAKE_FULL_COPY
        // Restore initial triangle list
        memcpy(mesh->triangles, triangleCopy, lastNumTriangles * sizeof(Triangle));
        mesh->numTriangles = lastNumTriangles;
#endif
    }

#ifdef MAKE_FULL_COPY
    // Free memory
    if (triangleCopy != NULL) free(triangleCopy);
#endif
    bh_free(h);

    h = nh;

#else

    // Compute vertices adjacent to a vertex
    lv = verticesAdjToVertex(mesh, c->v, &numV);

    // Compute edges adjacent to vertices
    le = edgesAdjToVertices(mesh, lv, numV, &numE);

    free(lv);

    // Delete adjacent edges from the heap
    for(i = 0; i <numE; i++) {
        e = le[i];

        mesh->edges[e].enable = FALSE;
        bh_delete(h, e);
    }

    // Recompute cost only for the adjacent edges
    for(i = 0; i <numE; i++) {
        e = le[i];
        
        cost = computeEdgeCost(mesh, e);

        // Add only valid edges
        if (cost != FLT_MAX) {
            bh_insert(h, e, cost);
            mesh->edges[e].enable = TRUE;
        }
        
#ifdef MAKE_FULL_COPY
        // Restore initial triangle list
        memcpy(mesh->triangles, triangleCopy, lastNumTriangles * sizeof(Triangle));
        mesh->numTriangles = lastNumTriangles;
#endif

    }

#ifdef MAKE_FULL_COPY
    // Free memory
    if (triangleCopy != NULL) free(triangleCopy);
#endif

    free(le);
#endif

#ifdef DEBUG 
    bh_mydump(mesh, h);
    //getchar();
#endif

    return h;
}

GLdouble VMI::computeEdgeLength(Vertex *vertices, int u, int v) {
    GLfloat ux, uy, uz, vx, vy ,vz, rx, ry, rz;
    
    ux = vertices[u].x;
    uy = vertices[u].y;
    uz = vertices[u].z;

    vx = vertices[v].x;
    vy = vertices[v].y;
    vz = vertices[v].z;

    rx = ux - vx;
    ry = uy - vy;
    rz = uz - vz;

    return sqrt((rx * rx) + (ry * ry) + (rz * rz));
}

///////////////////////////////////////////////////////////////////////////////

void VMI::swap(unsigned int *i, unsigned int *j) {
    unsigned int t;

    t = *i;
    *i = *j;
    *j = t;
}

void VMI::chooseBestEndPoints(Mesh *mesh, int e) {
    int v = mesh->edges[e].v;
    int u = mesh->edges[e].u;
    int numTu;
    int *Tu  = trianglesAdjToVertex(mesh, u, &numTu);
    int numTuv = 0;
    int *Tuv = (int *)malloc(numTu * sizeof(int));
    int numTv;
    int *Tv  = trianglesAdjToVertex(mesh, v, &numTv);
    int i, j, f, n, n0, n1, n2;
    float cost_u_v, cost_v_u, mincurve, curve, dot;
    
    //printItemList(Tv, numTv);
    //printItemList(Tu, numTu);
    
    // Compute Tuv
    for (i=0; i<(int)numTu; i++) {
        
        n = Tu[i];
        
        n0 = mesh->triangles[n].indices[0];
        n1 = mesh->triangles[n].indices[1];
        n2 = mesh->triangles[n].indices[2];
        
        if ((n0 == v) || (n1 == v) || (n2 == v))
            
            Tuv[numTuv++] = n;
    }
    //printItemList(Tuv, numTuv);
    
    curve = 0.0f;
    for (i=0; i<numTu; i++) {
        n = Tu[i];
        mincurve = 1.0f;
        for (j=0; j<numTuv; j++) {
            f = Tuv[j];
            dot = glmDot(mesh->triangles[f].normal, mesh->triangles[n].normal);
            //printf("dot: %f \n", dot);
            mincurve = MIN(mincurve, (1.0f - dot) / 2.0f);
        }
        //printf("mincurve: %f \n", mincurve);
        curve = MAX(curve, mincurve);
    }
    cost_u_v = curve;
    
    curve = 0.0f;
    for (i=0; i<numTv; i++) {
        n = Tv[i];
        mincurve = 1.0f;
        for (j=0; j<numTuv; j++) {
            f = Tuv[j];
            dot = glmDot(mesh->triangles[f].normal, mesh->triangles[n].normal);
            //printf("%f \n", dot);
            mincurve = MIN(mincurve, (1.0f - dot) / 2.0f);
        }
        curve = MAX(curve, mincurve);
    }
    cost_v_u = curve;
    
    //printf("e%d(%d,%d) cost: %f \n", e, u, v, cost_u_v);
    //printf("e%d(%d,%d) cost: %f \n", e, v, u, cost_v_u);
    
    if ( (cost_v_u + 0.000001) < cost_u_v) {
        //printf("Swap edge\n");
        //printf("e%d(%d,%d) cost: %f \n", e, u, v, cost_u_v);
        //printf("e%d(%d,%d) cost: %f \n", e, v, u, cost_v_u);
        swap(&mesh->edges[e].u, &mesh->edges[e].v);
    }

    free(Tuv);
    free(Tu);
    free(Tv);
}
///////////////////////////////////////////////////////////////////////////////

GLdouble VMI::computeEdgeCost(Mesh *mesh, int e) {
    GLdouble cost = 0.0, newI, sal = 1.0, /*length,*/ error;
    GLuint i;
    Change *c;

    chooseBestEndPoints(mesh, e);
    //getchar();

    c = createChange(mesh, e);
    
    // Apply the edge collapse
    computeChanges(mesh, c);
    
    // Compute cost only for boundary or manifold edges
    if ((c->numDel <3) && (c->numDel > 0)) {

        //length = computeEdgeLength(mesh->vertices, c->u, c->v);

        doChange(mesh, c);
        
        // Compute the cost of an edge collapse
        getProjectedAreas(histogram, numCameras);
        //printFullHistogram(histogram, numCameras, mesh->numTriangles);
        //getchar();
#ifdef SALIENCY
        sal = computeEdgeSaliency(mesh, c, percentile);
        //printf("  e:%d s: %f\n", c->e, sal);
#endif
        for (i=0;i <numCameras; i++) {

#ifdef KL // Kullback-Leibler
            newI = computeKL(mesh, histogram[i]);
#endif
#ifdef MI // Mutual Information
            newI = computeMI(mesh, histogram, numCameras, i);
#endif
#ifdef HE // Hellinger
            newI = computeHE(mesh, histogram[i]);
#endif
#ifdef CS // Chi-Square
            newI = computeCS(mesh, histogram[i]);
#endif

            error = ABS(initialIs[i] - newI);
            
            if (error > 0.0) 
                cost += ((error * cameras[i].weight * sal) /*+ (length * 0.05)*/);
        }

        undoChange(mesh, c);

    } else cost = FLT_MAX;


    deleteChange(c);

    return cost;
}

///////////////////////////////////////////////////////////////////////////////

void VMI::simplifyModel(Mesh *mesh, GLuint numDemandedTri)
{
        int                     e;
        Change  *c;
        FILE            *file                                   = NULL;
        FILE            *file_simp_seq  =       NULL;
        char            s[MAX_CHAR];
        double  cost                                            = 0.0;
        bheap_t *h                                                      = NULL;
        float           percent =       (40.0) / (float)(mesh->numTriangles - numDemandedTri);

        if (numDemandedTri > mesh->currentNumTriangles)
        {
                return;
        }

        // Save changes into a file
        if (bSaveLog == TRUE)
        {
#ifdef SALIENCY       
#ifdef KL // Kullback-Leibler
                sprintf(s,"%s_VKL_S.log", filename);
#endif
#ifdef MI // Mutual Information
                sprintf(s,"%s_VMI_S.log", filename);
#endif
#ifdef HE // Hellinger
                sprintf(s,"%s_VHE_S.log", filename);
#endif
#ifdef CS // Chi-Square
                sprintf(s,"%s_VCS_S.log", filename);
#endif
#else
#ifdef KL // Kullback-Leibler
                sprintf(s,"%s_VKL.log", filename);
#endif
#ifdef MI // Mutual Information
                sprintf(s,"%s_VMI.log", filename);
#endif
#ifdef HE // Hellinger
                sprintf(s,"%s_VHE.log", filename);
#endif
#ifdef CS // Chi-Square
                sprintf(s,"%s_VCS.log", filename);
#endif
#endif

                glmWriteOBJ(pmodel, s, GLM_NONE);

                /* open the file */
                file = fopen(s, "a+");
                
                if (!file)
                {
                        fprintf(stderr, "simplifyModel() failed: can't open file \"%s\" to write.\n", s);
                        exit(1);
                }
        }

        //      Open file of simplification sequence.
        file_simp_seq   =       fopen("SimplifSequence.txt", "w");
        
        if (!file_simp_seq)
        {
                fprintf(stderr,
                                                "simplifyModel() failed: ",
                                                "can't open file \"SimplifSequence\" to write.\n");
        }
        
        h = initHeap(mesh);

        printf("Starting simplification...\n");

        while (mesh->currentNumTriangles > numDemandedTri /*&& cost == 0.0*/)
        {
                // Get the edge that has the minimum cost
                e = bh_min(h);

                c = createChange(mesh, e);

                // Apply the edge collapse
                computeChanges(mesh, c);
                doChange(mesh, c);

                if (bSaveLog == TRUE) writeChange(file, c);

                //      Write Simplification sequence.
                saveSimplificationSequence(c);
                
                cost = h->a[h->p[e]].key;
                
                printf( "Edge collapse e%d(%d,%d) %f MIN d %d m %d\n",
                                                c->e,
                                                c->u,
                                                c->v,
                                                cost,
                                                c->numDel,
                                                c->numMod);

                mesh->edges[e].enable = FALSE;
                
                bh_delete(h, e);

                // Get projected areas after the edge collapse
                getProjectedAreas(histogram, numCameras);

                computeCameraIs(histogram, numCameras, initialIs);

                // Update the heap according to the edge collapse
                h = updateHeap(h, mesh, c);

                mUPB((float)c->numDel * percent);

                deleteChange(c);

                printf("t %d\n", mesh->currentNumTriangles);
        }

        if (bSaveLog == TRUE)
        {
                fclose(file);
                printf("Log file written...Ok\n");
        }

        //      Close simplification sequence file.
        fclose(file_simp_seq);

        bh_free(h);
}