gcoptimization.cpp

markov random field in matlab code

#include "energy.h"
#include "graph.h"
#include "GCoptimization.h"
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include "string.h"
#define MAX_INTT 1000000000


/**************************************************************************************/

void GCoptimization::initialize_memory()
{

    m_lookupPixVar = (PixelType *) new PixelType[m_nPixels];
    m_labelTable   = (LabelType *) new LabelType[m_nLabels];

    terminateOnError( !m_lookupPixVar || !m_labelTable ,"Not enough memory");

    for ( int i = 0; i < m_nLabels; i++ )
        m_labelTable[i] = i;
}

/**************************************************************************************/

void GCoptimization::commonGridInitialization(PixelType width, PixelType height, int nLabels)
{
    terminateOnError( (width < 0) || (height <0) || (nLabels <0 ),"Illegal negative parameters");

    m_width       = width;
    m_height      = height;
    m_nPixels     = width*height;
    m_nLabels     = nLabels;
    m_grid_graph  = 1;
        

}

/**************************************************************************************/

void GCoptimization::setParameters(int numParam, void *param)
{
    if (numParam != 1 ) 
        printf("\nInvalid number of parameters, can only set one parameter\\that is boolean label order\n");
    else
    {
        m_random_label_order = *((bool *) param);
    }
}

/**************************************************************************************/

void GCoptimization::commonNonGridInitialization(PixelType nupixels, int num_labels)
{
    terminateOnError( (nupixels <0) || (num_labels <0 ),"Illegal negative parameters");

    m_nLabels        = num_labels;
    m_nPixels        = nupixels;
    m_grid_graph         = 0;

    m_neighbors = (LinkedBlockList *) new LinkedBlockList[nupixels];

    terminateOnError(!m_neighbors,"Not enough memory");

}

/**************************************************************************************/

void GCoptimization::commonInitialization()
{
    m_needToFreeV        = 0;
    m_random_label_order = 1;
    initialize_memory();
}

/**************************************************************************************/
/* Use this constructor only for grid graphs                                          */
GCoptimization::GCoptimization(PixelType width,PixelType height,LabelType nLabels,EnergyFunction *eng):MRF(width,height,nLabels,eng)
{
    commonGridInitialization(width,height,nLabels);
    
    m_labeling           = (LabelType *) new LabelType[m_nPixels];
    terminateOnError(!m_labeling,"out of memory");
    for ( int i = 0; i < m_nPixels; i++ ) m_labeling[i] = (LabelType) 0;

    commonInitialization(); 
}

/**************************************************************************************/
/* Use this constructor for general graphs                                            */
GCoptimization::GCoptimization(PixelType nPixels,int nLabels,EnergyFunction *eng):MRF(nPixels,nLabels,eng) 
{

    commonNonGridInitialization(nPixels, nLabels);
    
    m_labeling           = (LabelType *) new LabelType[m_nPixels];
    terminateOnError(!m_labeling,"out of memory");
    for ( int i = 0; i < nPixels; i++ ) m_labeling[i] = (LabelType) 0;

    commonInitialization(); 
}

/**************************************************************************************/

void GCoptimization::setData(EnergyTermType* dataArray)
{
    m_datacost  = dataArray;
}


/**************************************************************************************/

void GCoptimization::setData(DataCostFn dataFn)
{
    m_dataFnPix = dataFn;

}

/**************************************************************************************/

void GCoptimization::setSmoothness(EnergyTermType* V)
{
    m_smoothcost = V;
}

/**************************************************************************************/

void GCoptimization::setSmoothness(int smoothExp,CostVal smoothMax, CostVal lambda)
{
    int i, j;
    CostVal cost;

    m_needToFreeV = 1;

    m_smoothcost = (CostVal *) new CostVal[m_nLabels*m_nLabels*sizeof(CostVal)];
    if (!m_smoothcost) { fprintf(stderr, "Not enough memory!\n"); exit(1); }


    for (i=0; i<m_nLabels; i++)
        for (j=i; j<m_nLabels; j++)
        {
            cost = (CostVal) ((smoothExp == 1) ? j - i : (j - i)*(j - i));
            if (cost > smoothMax) cost = smoothMax;
            m_smoothcost[i*m_nLabels + j] = m_smoothcost[j*m_nLabels + i] = cost*lambda;
        }

}

/**************************************************************************************/

void GCoptimization::setCues(EnergyTermType* hCue, EnergyTermType* vCue)
{

    m_horizWeights    = hCue;
    m_vertWeights     = vCue;
}


/**************************************************************************************/

void GCoptimization::setSmoothness(SmoothCostGeneralFn cost)
{
    m_smoothFnPix   = cost;
}

/**************************************************************************************/

GCoptimization::EnergyType GCoptimization::dataEnergy()
{
    
    if ( m_dataType == ARRAY) 
        return(giveDataEnergyArray());
    else return(giveDataEnergyFnPix());

    return(0);
}

/**************************************************************************************/
    
GCoptimization::EnergyType GCoptimization::giveDataEnergyArray()
{
    EnergyType eng = (EnergyType) 0;


    for ( int i = 0; i < m_nPixels; i++ )
        eng = eng + m_datacost(i,m_labeling[i]);

    return(eng);
}

/**************************************************************************************/
    
GCoptimization::EnergyType GCoptimization::giveDataEnergyFnPix()
{

    EnergyType eng = (EnergyType) 0;

    for ( int i = 0; i < m_nPixels; i++ )
        eng = eng + m_dataFnPix(i,m_labeling[i]);

    return(eng);
}

/**************************************************************************************/

GCoptimization::EnergyType GCoptimization::smoothnessEnergy()
{
    
    if ( m_grid_graph )
    {
        if ( m_smoothType != FUNCTION )
        {
            if (m_varWeights) return(giveSmoothEnergy_G_ARRAY_VW());
            else return(giveSmoothEnergy_G_ARRAY());
        }
        else return(giveSmoothEnergy_G_FnPix());
    }
    else
    {
        if ( m_smoothType != FUNCTION ) return(giveSmoothEnergy_NG_ARRAY());
        else  return(giveSmoothEnergy_NG_FnPix());
    }
    return(0);

}

/**************************************************************************************/

GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_NG_FnPix()
{

    EnergyType eng = (EnergyType) 0;
    int i;
    Neighbor *temp; 

    for ( i = 0; i < m_nPixels; i++ )
        if ( !m_neighbors[i].isEmpty() )
        {
            m_neighbors[i].setCursorFront();
            while ( m_neighbors[i].hasNext() )
            {
                temp = (Neighbor *) m_neighbors[i].next();
                if ( i < temp->to_node )
                    eng = eng + m_smoothFnPix(i,temp->to_node, m_labeling[i],m_labeling[temp->to_node]);
            }
        }
        
    return(eng);
    
}

/**************************************************************************************/

GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_NG_ARRAY()
{
    EnergyType eng = (EnergyType) 0;
    int i;
    Neighbor *temp; 

    for ( i = 0; i < m_nPixels; i++ )
        if ( !m_neighbors[i].isEmpty() )
        {
            m_neighbors[i].setCursorFront();
            while ( m_neighbors[i].hasNext() )
            {
                temp = (Neighbor *) m_neighbors[i].next();

                if ( i < temp->to_node )
                    eng = eng + m_smoothcost(m_labeling[i],m_labeling[temp->to_node])*(temp->weight);

            }
        }

    return(eng);
}

/**************************************************************************************/

GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_G_ARRAY_VW()
{

    EnergyType eng = (EnergyType) 0;
    int x,y,pix;

    
    for ( y = 0; y < m_height; y++ )
        for ( x = 1; x < m_width; x++ )
        {
            pix = x+y*m_width;
            eng = eng + m_smoothcost(m_labeling[pix],m_labeling[pix-1])*m_horizWeights[pix-1];
        }

    for ( y = 1; y < m_height; y++ )
        for ( x = 0; x < m_width; x++ )
        {
            pix = x+y*m_width;
            eng = eng + m_smoothcost(m_labeling[pix],m_labeling[pix-m_width])*m_vertWeights[pix-m_width];
        }

    
    return(eng);
    
}
/**************************************************************************************/

GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_G_ARRAY()
{

    EnergyType eng = (EnergyType) 0;
    int x,y,pix;


    for ( y = 0; y < m_height; y++ )
        for ( x = 1; x < m_width; x++ )
        {
            pix = x+y*m_width;
            eng = eng + m_smoothcost(m_labeling[pix],m_labeling[pix-1]);
        }

    for ( y = 1; y < m_height; y++ )
        for ( x = 0; x < m_width; x++ )
        {
            pix = x+y*m_width;
            eng = eng + m_smoothcost(m_labeling[pix],m_labeling[pix-m_width]);
        }

    return(eng);
    
}
/**************************************************************************************/

GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_G_FnPix()
{

    EnergyType eng = (EnergyType) 0;
    int x,y,pix;


    for ( y = 0; y < m_height; y++ )
        for ( x = 1; x < m_width; x++ )
        {
            pix = x+y*m_width;
            eng = eng + m_smoothFnPix(pix,pix-1,m_labeling[pix],m_labeling[pix-1]);