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

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

#define LOG2 (double)0.69314718055994530941723212145818
#define log2(x) (log((x))/LOG2)

#include "../include/metrics.h"
#include "../include/histogram.h"
#include "../include/global.h"

using namespace VMI;

GLdouble VMI::computeMeanProjAreaNoBG(GLuint **histogram, GLuint numCameras, int t) {
    GLuint i;
    GLdouble mean_proj_area = 0.0;

    for (i=0; i<numCameras; i++) {
        mean_proj_area += histogram[i][t];
    }

    return (mean_proj_area / (GLdouble)numCameras);
}

GLdouble VMI::computeMeanProjArea(GLuint **histogram, GLuint numCameras, int t) {
        GLuint i;
        GLdouble mean_proj_area = 0.0, total_proj_area, resolution = width * height;

        for (i=0; i<numCameras; i++) {
                total_proj_area = resolution - histogram[i][0];

                mean_proj_area += histogram[i][t] / total_proj_area;
        }

        return (mean_proj_area / (GLdouble)numCameras);
}
///////////////////////////////////////////////////////////////////////////////
// Mutual Information
GLdouble VMI::computeMI(Mesh *mesh, GLuint **histogram, GLuint numCameras, GLuint cam) {
    GLdouble I = 0.0, pov, po, total_proj_area = width * height;
    GLuint i;
    
    for (i=0; i<mesh->numTriangles; i++) { 
        
        if (mesh->triangles[i].enable == TRUE) {
            
            pov = (GLdouble)histogram[cam][i + 1] / total_proj_area;
            
            po = computeMeanProjAreaNoBG(histogram, numCameras, i + 1) / total_proj_area;
            
            if ((pov > 0.0) && (po > 0.0))
                I += (pov * log2(pov / po));
        }
    }

    // take into account the background
    pov = (GLdouble)histogram[cam][0] / total_proj_area;
    
    po = computeMeanProjAreaNoBG(histogram, numCameras, 0) / total_proj_area;
    
    if ((pov > 0.0) && (po > 0.0))
        I += (pov * log2(pov / po));
    
    return (I / numCameras);
}

GLdouble VMI::decMI(GLdouble I, GLuint **histogram, GLuint numCameras, GLuint cam, Change *c) {
    GLdouble newI = I * numCameras, 
             pov = 0.0, po;
    GLsizei total_proj_area = width * height;
    GLint i, t;

    // decrease entropy of deleted triangles
    for (i=0; i<c->numDel; i++) { 
        
        t = c->deleted[i].id;
        // projected pixels of triangle t is at t + 1 position
        pov = (GLdouble)histogram[cam][t + 1] / total_proj_area;
        
        po = computeMeanProjAreaNoBG(histogram, numCameras, t + 1) / total_proj_area;
        
        if ((pov > 0.0) && (po > 0.0))
            newI -= (pov * log2(pov / po));
    }

    // decrease entropy of modified triangles
    for (i=0; i<c->numMod; i++) { 
        
        t = c->modified[i].id;
        // projected pixels of triangle t is at t + 1 position
        pov = (GLdouble)histogram[cam][t + 1] / total_proj_area;
        
        po = computeMeanProjAreaNoBG(histogram, numCameras, t + 1) / total_proj_area;
        
        if ((pov > 0.0) && (po > 0.0))
            newI -= (pov * log2(pov / po));
    }

    // take into account the background
    pov = (GLdouble)histogram[cam][0] / total_proj_area;
    
    po = computeMeanProjAreaNoBG(histogram, numCameras, 0) / total_proj_area;
    
    if ((pov > 0.0) && (po > 0.0))
        newI -= (pov * log2(pov / po));

    return (newI / numCameras);
}

GLdouble VMI::incMI(GLdouble I, GLuint **histogram, GLuint numCameras, GLuint cam, Change *c) {
    GLdouble newI = I * numCameras, 
             pov = 0.0, po;
    GLsizei total_proj_area = width * height;
    GLuint i, t;

    // increase entropy of modified triangles
    for (i=0; i<(GLuint)c->numMod; i++) { 
        
        t = c->modified[i].id;
        // projected pixels of triangle t is at t + 1 position
        pov = (GLdouble)histogram[cam][t + 1] / total_proj_area;
        
        po = computeMeanProjAreaNoBG(histogram, numCameras, t + 1) / total_proj_area;
        
        if ((pov > 0.0) && (po > 0.0))
            newI += (pov * log2(pov / po));
    }

    // take into account the background
    pov = (GLdouble)histogram[cam][0] / total_proj_area;
    
    po = computeMeanProjAreaNoBG(histogram, numCameras, 0) / total_proj_area;
    
    if ((pov > 0.0) && (po > 0.0))
        newI += (pov * log2(pov / po));

    return (newI / numCameras);
}
///////////////////////////////////////////////////////////////////////////////
// Kullback-Leibler divergence
GLdouble VMI::computeKL(Mesh *mesh, GLuint *histogram) {
    GLdouble I = 0.0, 
             total_real_area = 0.0, pi, tri_area;
    GLsizei  total_proj_area = (width * height) - histogram[0];
    GLuint   i;

    // Compute total real area of all triangles
    for (i=0; i<mesh->numTriangles; i++) {
        if (mesh->triangles[i].enable == TRUE) {
            total_real_area += mesh->triangles[i].area;
        }
    }

    for (i=0; i<mesh->numTriangles; i++) { // don't take into account the background
        
        if (mesh->triangles[i].enable == TRUE) {
            pi = (GLdouble)histogram[i + 1] / total_proj_area;
            
            tri_area = mesh->triangles[i].area;
            
            if ((pi > 0.0) && (tri_area != 0.0))
                I += (pi * log2((pi * total_real_area) / tri_area));
        }
    }

    return I;
}

// Hellinger divergence (square root)
GLdouble VMI::computeHE(Mesh *mesh, GLuint *histogram) {
    GLdouble I = 0.0, temp, 
             total_real_area = 0.0, pi, qi;
    GLsizei  total_proj_area = (width * height) - histogram[0];
    GLuint i;
    
    // Compute total real area of all triangles
    for (i=0; i<mesh->numTriangles; i++) {
        if (mesh->triangles[i].enable == TRUE) {
            total_real_area += mesh->triangles[i].area;
        }
    }

    // sqrt( 1/2 * \Sum_i (sqrt(p_i) - sqrt(q_i))^2)
    for (i=0; i<mesh->numTriangles; i++) { // + 1 because the first is the background
        if (mesh->triangles[i].enable == TRUE) {
            pi = (GLdouble)histogram[i + 1] / total_proj_area;
            
            qi = mesh->triangles[i].area / total_real_area;
            
            if ((pi > 0.0) && (qi > 0.0)) {
                
                temp = sqrt(pi) - sqrt(qi);
                
                I += (temp * temp);
            }
        }
    }

    return sqrt(I / 2.0);

}

// Chi-Square divergence (square root)
GLdouble VMI::computeCS(Mesh *mesh, GLuint *histogram) {
    GLdouble I = 0.0, 
             total_real_area = 0.0, pi, qi;
    GLsizei  total_proj_area = (width * height) - histogram[0];
    GLuint i;
    
    /// Compute total real area of all triangles
    for (i=0; i<mesh->numTriangles; i++) {
        if (mesh->triangles[i].enable == TRUE) {
            total_real_area += mesh->triangles[i].area;
        }
    }

    // sqrt( \Sum_i ((p_i/A) - (p'_i/A'))^2 / (p'_i/A')))
    for (i=0; i<mesh->numTriangles; i++) { // + 1 because the first is the background
        if (mesh->triangles[i].enable == TRUE) {
            pi = (GLdouble)histogram[i + 1] / total_proj_area;
            
            qi = mesh->triangles[i].area / total_real_area;
            
            if ((qi > 0.0) && (pi > 0.0))
                I += (((pi - qi) * (pi - qi)) / qi);
        }
    }

    return sqrt(I);
}
///////////////////////////////////////////////////////////////////////////////
GLdouble VMI::computeEntropy(GLuint **histogram, GLuint numCameras, GLuint k) {
        GLdouble H = 0.0, POk, pi,
                                         total_proj_area, resolution = width * height, pv = 1.0 / (GLdouble)numCameras;
        GLuint i;

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

                total_proj_area = resolution - histogram[i][0];

                POk = computeMeanProjArea(histogram, numCameras, k + 1);

                if (POk > 0.0)
                        pi = (pv * ((double)histogram[i][k + 1]) / total_proj_area ) / POk;
                else pi = 0.0;

                if (pi > 0.0)
                        H += (pi * log2(pi));
        }

        return -H;
}

GLdouble VMI::computeMixedEntropy(GLdouble *mixed, GLuint numCameras) {
    GLdouble H = 0.0, pi;
    GLuint i;
    
    for (i=0; i<numCameras; i++) {
        
        pi = mixed[i];
        
        if (pi > 0.0)
            H += (pi * log2(pi));
    }
    
    return -H;
}
///////////////////////////////////////////////////////////////////////////////
// Jensen-Shannon divergence
GLdouble VMI::computeJS(GLuint **histogram, GLuint numCameras, GLuint j, GLuint k) {
        GLdouble js, Wj = 0.0, Wk = 0.0, POj, POk, pj, pk, total_proj_area,
                                         resolution = width * height, pv = 1.0 / (GLdouble)numCameras,
                                         *mixing = (GLdouble *)malloc(numCameras * sizeof(GLdouble));
        GLuint i;

        POj = computeMeanProjArea(histogram, numCameras, j + 1);
        POk = computeMeanProjArea(histogram, numCameras, k + 1);

        //printf("%f %f\n", POj, POk);

        // weights
        if ((POj + POk) > 0.0) {
                Wj = POj / (POj + POk);
                Wk = POk / (POj + POk);
        }

        //printf("%f %f\n", Wj, Wk);

        // Compute mixing distribution
        for (i=0; i<numCameras; i++) {

                total_proj_area = resolution - histogram[i][0];

                if (POj > 0.0) pj = (pv * ((double)histogram[i][j + 1] / total_proj_area)) / POj;
                else pj = 0.0;

                if (POk > 0.0) pk = (pv * ((double)histogram[i][k + 1] / total_proj_area)) / POk;
                else pk = 0.0;

                mixing[i] = (Wj * pj) + (Wk * pk);
        }

        js = /*(POj + POk) * */(computeMixedEntropy(mixing, numCameras) 
                        - (Wj * computeEntropy(histogram, numCameras, j)) 
                        - (Wk * computeEntropy(histogram, numCameras, k)));

        free(mixing);

        return js;
}
///////////////////////////////////////////////////////////////////////////////
void VMI::getProjectedAreas(GLuint **histogram, GLuint numCameras) {

    renderScene(histogram, numCameras);
}

void VMI::getProjectedAreasWin(GLuint **histogram, GLuint numCameras, Change *c) {

    renderSceneWin(histogram, numCameras, c);
}

void VMI::resetProjectedAreas(GLuint **histogram, GLuint numCameras) {
    GLuint i = 0;

    // Reset the projected areas for all cameras
    for (i=0; i<numCameras; i++)
        resetSWHistogram(histogram[i], mesh->numTriangles);
}

GLdouble *VMI::initIs(GLuint numCameras) {
    GLdouble *initialIs;

     // Allocating memory for the metric
    initialIs = (GLdouble *)malloc(sizeof(GLdouble) * numCameras);
    
    if (initialIs == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }

    // Set metric to 0.0
    memset(initialIs, 0, sizeof(GLdouble) * numCameras);

    return initialIs;
}

void VMI::computeCameraIs(GLuint **histogram, GLuint numCameras, GLdouble *mis) {
    GLuint i = 0;
    GLdouble meanI = 0.0;

    for (i=0;i <numCameras; i++) {
        
#ifdef KL // Kullback-Leibler
        mis[i] = computeKL(mesh, histogram[i]);
#endif
#ifdef MI // Mutual Information
        mis[i] = computeMI(mesh, histogram, numCameras, i);
#endif
#ifdef HE // Hellinger
        mis[i] = computeHE(mesh, histogram[i]);
#endif
#ifdef CS // Chi-Square
        mis[i] = computeCS(mesh, histogram[i]);
#endif
        
        meanI += mis[i];
    }
#ifndef MI 
    meanI /= numCameras;
#endif
    //printIs(mis, numCameras);
    printf("I0= %f\n", meanI); 
}

void VMI::printIs(GLdouble *mis, GLuint numCameras) {
    GLuint i;

    printf("\n");
    for (i=0; i<numCameras; i++) {
    
        printf("Camera %d  I = %f\n", i, mis[i]);
    }
}