JointHistogram.cpp

#include "JointHistogram.hh"
#include <Base/Math/Random.hh>
#include <Filter/Gauss.hh>

using namespace BIAS;
using namespace std;

template<class StorageType>
JointHistogram<StorageType>::JointHistogram(){


}
template<class StorageType>
JointHistogram<StorageType>::~JointHistogram(){


}

template<class StorageType> int JointHistogram<StorageType>::
ComputeParzenWindow(BIAS::Image<StorageType>& image1,
        BIAS::Image<StorageType>& image2,
        unsigned int nrOfSamples)
{

  dNumberOfSamples_ = nrOfSamples;

  BIAS::Gauss<StorageType,float> gaussFilter_;
  BIAS::Random randomizer;

  dWidth_ = image1.GetWidth();
  dHeight_ = image1.GetHeight();
  dChannelCount_= image1.GetChannelCount();
  
  if(dWidth_ != image2.GetWidth() || dHeight_!= image2.GetHeight()
     || dChannelCount_ !=image2.GetChannelCount() ){
    BIASERR("Give Images of equal dimension!");
    return -1;
  }

  //gaussfilter the input image
  gaussFilter_.Filter7x7Grey(image1, gaussImage1_);
  gaussFilter_.Filter7x7Grey(image2, gaussImage2_);

  // init the histogram image
  if(!histImage_.IsEmpty() || !nrImage_.IsEmpty() ){
    histImage_.Release();  
    nrImage_.Release();
  }
  
  float min1,max1,min2,max2;  
  gaussImage1_.GetMinMaxPixelValue(min1,max1);
  gaussImage2_.GetMinMaxPixelValue(min2,max2);

  if(min2>min1) min1=min2;
  if(max2>max1) max1=max2;

  histImage_.Init((unsigned int)(rint(max1+1)),(unsigned int)(rint(max1+1)),1);
  nrImage_.Init((unsigned int)(rint(max1+1)),(unsigned int)(rint(max1+1)),1);

  float ** ida1 = NULL;
  float ** ida2 = NULL; 
  float ** idaH = NULL;
  float ** idaNr = NULL;
 
  ida1 = gaussImage1_.GetImageDataArray(); //take gaussfiltered image 
  ida2 = gaussImage2_.GetImageDataArray(); //take gaussfiltered image
  idaNr = nrImage_.GetImageDataArray();
  idaH = histImage_.GetImageDataArray();

  Vector2<unsigned int> coordsTmp;

  for(unsigned int i=0;i<dNumberOfSamples_;i++){
    coordsTmp[0] = randomizer.GetUniformDistributedInt(0,(int)dWidth_-1);
    coordsTmp[1] = randomizer.GetUniformDistributedInt(0,(int)dHeight_-1);  
 
    float pixel1=0.0,pixel2=0.0;
    pixel1 = ida1[coordsTmp[1]][coordsTmp[0]]; // get value at y,x  
    pixel2 = ida2[coordsTmp[1]][coordsTmp[0]]; // get value at y,x  

    idaH[(unsigned int)(pixel1)][(unsigned int)(pixel2)]+=1.0;
    idaNr[(unsigned int)(pixel1)][(unsigned int)(pixel2)] +=1.0;
  }
  
  NormalizeParzen_();

  return 0;

}



template<class StorageType> int JointHistogram<StorageType>::
Compute(BIAS::Image<StorageType>& image1,
        BIAS::Image<StorageType>& image2)
{
  dWidth_ = image1.GetWidth();
  dHeight_ = image1.GetHeight();
  dChannelCount_= image1.GetChannelCount();

  if(dWidth_ != image2.GetWidth() || dHeight_!= image2.GetHeight()
     || dChannelCount_ !=image2.GetChannelCount() ){
    BIASERR("Give Images of equal dimension!");
    return -1;
  }
  
  if(!histImage_.IsEmpty())
    histImage_.Release();
  
  StorageType min1,max1,min2,max2;  
  image1.GetMinMaxPixelValue(min1,max1);
  image2.GetMinMaxPixelValue(min2,max2);

  if(min2>min1) min1=min2;
  if(max2>max1) max1=max2;
  
  if(!histImage_.IsEmpty())
    histImage_.Release();
  histImage_.Init((unsigned int)(max1+1),(unsigned int)(max1+1),1);
  
  float ** idaHist = histImage_.GetImageDataArray();
  StorageType ** ida1 = image1.GetImageDataArray();
  StorageType ** ida2 = image2.GetImageDataArray();

  //count total number of pixels
  pixelCount_=0;

  for(unsigned int x=0;x< dWidth_;x++){
    for(unsigned int y=0;y< dHeight_;y++){
   
      idaHist[(unsigned int)(ida1[y][x])][(unsigned int)(ida2[y][x])]+=1.0;
      pixelCount_++;
    }
  }
  drawImage_ = histImage_;
  //normalize image to enable Shannon entropy
  Normalize_();
  
  return 0;
}

template<class StorageType> void JointHistogram<StorageType>::
Draw(BIAS::Image<float>& image){
  image = drawImage_;
}

template<class StorageType> float JointHistogram<StorageType>::
CalcShannonEntropy(){
  float entropy=0.0;
  unsigned int width = histImage_.GetWidth();
  unsigned int height = histImage_.GetHeight();
 
  float ** idaHist = histImage_.GetImageDataArray();

  for(unsigned int x=0;x< width;x++){
    for(unsigned int y=0;y< height;y++){

      if(idaHist[y][x] != 0.0){
        entropy += (idaHist[y][x]*log(idaHist[y][x]));
       
      }
    }
  }
  return (-1)*entropy;
}

template<class StorageType> float JointHistogram<StorageType>::
CalcRenyiEntropy(double alpha){
  float entropy=0.0;
  unsigned int width = histImage_.GetWidth();
  unsigned int height = histImage_.GetHeight();
 
  float ** idaHist = histImage_.GetImageDataArray();

  for(unsigned int x=0;x< width;x++){
    for(unsigned int y=0;y< height;y++){

      if(idaHist[y][x] != 0.0){
        entropy += pow((float)idaHist[y][x],(float)alpha);
       
      }
    }
  }
  return (float)( (1/(1-alpha))*log(entropy) ); 
}


template<class StorageType> void JointHistogram<StorageType>::
NormalizeParzen_(){
  unsigned int width = histImage_.GetWidth();
  unsigned int height = histImage_.GetHeight();

  float ** idaHist = histImage_.GetImageDataArray();
  float ** idaNr = nrImage_.GetImageDataArray();
  for(unsigned int x=0;x< width;x++){
    for(unsigned int y=0;y< height;y++){
      if(idaNr[(int)y][(int)x] != 0.0)
        idaHist[y][x] /=  idaNr[(int)y][(int)x]*dNumberOfSamples_;
      else if(idaHist[y][x] != 0.0)
        idaHist[y][x] /=dNumberOfSamples_;
        
    }
  }
}

template<class StorageType> void JointHistogram<StorageType>::
Normalize_(){
  unsigned int width = histImage_.GetWidth();
  unsigned int height = histImage_.GetHeight();

  float ** idaHist = histImage_.GetImageDataArray();

  for(unsigned int x=0;x< width;x++){
    for(unsigned int y=0;y< height;y++){
      idaHist[y][x] /= (float)pixelCount_;
    }
  }
}

#define INST(type) template class JointHistogram<type>;
namespace BIAS{
INST(float)
INST(unsigned char)
}