source: GTP/trunk/Lib/Geom/shared/GTGeometry/src/libs/vmi/src/saliency.cpp @ 2090

#include <stdio.h>
#include <stdlib.h>

#include <math.h>

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

#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))

#define MAX_NUM_TRI 450

namespace VMI
{
double maxSal = 0.0, minSal = 0.0;
int alpha = 30; // 30th percentile
int lambda = /*100*/1;
double percentile = 0.0; // percentile
}

using namespace VMI;

void VMI::computeSaliency(Mesh *mesh, GLuint **histogram, GLuint numCameras) {
        GLuint i = 0;
        double sal;

        maxSal = minSal = computeTriangleSaliency(mesh, histogram, numCameras, 0);

        mesh->triangles[0].saliency = maxSal;

        for (i=1; i<mesh->numTriangles; i++) {

                sal = computeTriangleSaliency(mesh, histogram, numCameras, i);

                if (sal > maxSal) maxSal = sal;
                if (sal < minSal) minSal = sal;

                mesh->triangles[i].saliency = sal;
        }

        printf("\nMax Sal: %f  Min Sal: %f\n", maxSal, minSal);

        percentile = computePercentile(alpha, mesh->triangles, mesh->numTriangles);
        printf("\n%dth percentile: %f\n", alpha, percentile);
}

void updateTriangleSaliency(Mesh *mesh, GLuint **histogram, GLuint numCameras, int v) {
        int i, t;

        for(i=0; i<(int)mesh->vertices[v].numTriangles; i++) {

                t = mesh->vertices[v].triangles[i];

                mesh->triangles[t].saliency = computeTriangleSaliency(mesh, histogram, numCameras, t);
        }
}

double VMI::computeTriangleSaliency(Mesh *mesh, GLuint **histogram, GLuint numCameras, GLuint k) {
     GLuint i, l, v0, v1, v2;
    int *triangles = NULL, n = 0;
    double sal = 0.0;
    
    v0 = mesh->triangles[k].indices[0];
    v1 = mesh->triangles[k].indices[1];
    v2 = mesh->triangles[k].indices[2];
    
    /* Allocating memory */
    triangles = (int *)malloc((mesh->vertices[v0].numTriangles +
                               mesh->vertices[v1].numTriangles +
                               mesh->vertices[v2].numTriangles) * sizeof(int));

    for(i=0; i<mesh->vertices[v0].numTriangles; i++) {
        l = mesh->vertices[v0].triangles[i];
        
        if (l != k) 
            addItem(triangles, &n, l);
    }
    
    for(i=0; i<mesh->vertices[v1].numTriangles; i++) {
        l = mesh->vertices[v1].triangles[i];
        
        if (l != k) 
            addItem(triangles, &n, l);
    }
    
    for(i=0; i<mesh->vertices[v2].numTriangles; i++) {
        l = mesh->vertices[v2].triangles[i];
        
        if (l != k) 
            addItem(triangles, &n, l);
    }
    
    //printItemList(triangles, n);
    
    for(i=0; i<(GLuint)n; i++) {
        //printf("\n%d %d", k, triangles[i]);
        //sal += computeJS(histogram, numCameras, k, triangles[i]);
                sal = MAX ( sal, computeJS(histogram, numCameras, k, triangles[i]));
    }
    //printf("\nT%d Sal: %f\n", k, (sal / n));
    //getchar();
    
    free(triangles);

    return (sal /*/ n*/);
}

double VMI::computeEdgeSaliency(Mesh *mesh, Change *c, double p) {
    int i, t;
    double sal, edgeSal = 0.0;
    
    // edge saliency
    for(i=0; i<c->numDel; i++) {

        t = c->deleted[i].id;

        sal = mesh->triangles[t].saliency;

        //printf("\nt%d: osal: %f\n",t, mesh->triangles[t].saliency);

        // weight map derived from saliency map
        if (sal>= p) sal *= lambda;

        //printf("\nt%d: nsal: %f\n",t, sal);

        edgeSal += sal;
    }
    
    //printf("\ne:%d sal: %f\n",c->e, edgeSal);

    //getchar();

    return edgeSal;
}

void VMI::saveSaliencyMap(Mesh *mesh, char* filename) {
    FILE *fp;
    GLuint i;


    if((fp= fopen(filename, "wt"))== NULL) {
        printf("Can't open file %s\n", filename);
        exit(1);
    }

    printf("Writing file %s\n", filename);

    fprintf (fp, "%d\n", mesh->numTriangles);
    //printf("%d\n", num);

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

        fprintf(fp, "t%d %f\n",  i, mesh->triangles[i].saliency);

    }

    fclose(fp);
    printf("File writen.\n");
}

void VMI::loadSaliencyMap(Mesh *mesh, char* filename) {
    FILE *fp;
    unsigned int i, num;
    float s;


    if((fp= fopen(filename, "rt"))== NULL) {
        printf("Can't open file %s\n", filename);
                                getchar();
                                exit(1);
    }

    printf("Reading file %s\n", filename);

    maxSal = minSal = 0.0;

    // Reading number of triangles
    fscanf (fp, "%d\n", &num);
    //printf("%d\n", num);

    while (!feof(fp)) {

        fscanf(fp, "t%d %f\n",  &i, &s);

        if (s > maxSal) maxSal = s;
        if (s < minSal) minSal = s;

        // Save the saliency into the triangle i
        mesh->triangles[i].saliency = s;

    }

    fclose(fp);
    printf("File read.\n");

    printf("\nMax Sal: %f  Min Sal: %f\n", maxSal, minSal);

    percentile = computePercentile(alpha, mesh->triangles, mesh->numTriangles);
    printf("\n%dth percentile: %f\n", alpha, percentile);
}

void VMI::viewSaliency(Mesh *mesh, Camera *cameras, int cam) {
    GLuint i, v1, v2, v3;
    float r, g, b;
    
    // Clear color and depth buffers
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    
    // Camera  i
    gluLookAt(cameras[cam].eyeX, cameras[cam].eyeY, cameras[cam].eyeZ,
              cameras[cam].centerX, cameras[cam].centerY, cameras[cam].centerZ,
              cameras[cam].upX, cameras[cam].upY, cameras[cam].upZ);
    
    // Immediate mode
    glBegin (GL_TRIANGLES);
    for (i=0; i<mesh->numTriangles; i++) {
        if (mesh->triangles[i].enable == TRUE) {
            
            v1= mesh->triangles[i].indices[0];
            v2= mesh->triangles[i].indices[1];
            v3= mesh->triangles[i].indices[2];
            
            computeRGB(minSal, maxSal, mesh->triangles[i].saliency, &r, &g, &b);
            
            glColor3f(r, g, b);
            glVertex3f(mesh->vertices[v1].x, mesh->vertices[v1].y, mesh->vertices[v1].z);
            glVertex3f(mesh->vertices[v2].x, mesh->vertices[v2].y, mesh->vertices[v2].z);
            glVertex3f(mesh->vertices[v3].x, mesh->vertices[v3].y, mesh->vertices[v3].z);
        }
    }
    glEnd();
}

void VMI::computeRGB(double min, double max,double value,float *r,float *g,float *b)
{
        /* Given the maximum, the minimum value and the demanded color value returns
         * the respective RGB color
         */
        if(value>max) value=max;
        if(value<min) value=min;

        if(max==min) value=1;
        else value=(float)(value-min)/(float)(max-min);

        if (value<=0.25)
        {
                (*r)=0;
                (*g)=4*value;
                (*b)=1;
        }
        else if (value<=0.50)
        {
                (*r)=0;
                (*g)=1;
                (*b)=2-4*value;
        }
        else if (value<=0.75)
        {
                (*r)=4*value-2;
                (*g)=1;
                (*b)=0;
        }
        else
        {
                (*r)=1;
                (*g)=4-4*value;
                (*b)=0;
        }
}

int VMI::compare(const void *arg1, const void *arg2)
{
    if(*(double *)arg1 < *(double *)arg2) return -1;
    else if(*(double *)arg1 > *(double *)arg2) return 1;
    else return 0;
}

int VMI::isDecimal(double value) {
    long r = round(value);

    //printf ("%d\n",r);

    if ((double)r != value) return TRUE;
    else return FALSE;
}

double VMI::computePercentile(int q, Triangle *triangles, int numTriangles) {
    double p = (numTriangles * q) / 100.0;
    double *list = (double *)malloc(numTriangles * sizeof(double)), l, per = 0.0;
    int i;
    
    //printf("%f\n",p);

    // filling the list
    for(i = 0; i < numTriangles; i++) 
        list[i] = triangles[i].saliency;
    
    /* Quick-Sort */
    qsort(list, numTriangles, sizeof(list[0]), compare);
    
    /* showing the sorted list */
    /*for(i = 0; i < numTriangles; i++) 
        printf("%f\n", list[i]);
    printf("\n");*/

    if (isDecimal(p) == TRUE) per = list[(int)p + 1];
    else {
        l = (list[(int)p + 1] - list[(int)p]) * (q /100.0);

        per = list[(int)p] + l;
    }

    free(list);
    
    return per;
}