IntegralHistogram.cpp

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


#include "IntegralHistogram.hh"
#include <Base/Image/ImageBase.hh>
#include <Base/Image/ImageConvert.hh>
#include <Base/Debug/TimeMeasure.hh>

#include <Base/Common/BIASpragma.hh>

using namespace BIAS;
using namespace std;

BIAS::TimeMeasure usertime, rowtime;

/**
   This class creates an IntegralHistogram in the HSL color space.
   @author ingo schiller
   @date feb 2004 
 */
IntegralHistogram::IntegralHistogram()
{
  bExpDist_=true;
  nMinSat_=10;
  nLambda_=20;
  nBinsize_=32;
  nBinsizeSat_=64;
  nOldBinsize_=0;

}

IntegralHistogram::~IntegralHistogram() 
{
  Clear();
}

void IntegralHistogram::SetLambda(int lambda){
    nLambda_= lambda;
}

void IntegralHistogram::SetMinSat(int minsat){
    nMinSat_=minsat;
}

void IntegralHistogram::SetBinSize(int binsize){
  nBinsize_ = binsize;  
}
void IntegralHistogram::SetExpDist(bool expdist){
    bExpDist_=expdist;
}

void IntegralHistogram::Clear()
{ 
  bExpDist_=true;
  nMinSat_=10;
  nLambda_=20;
  nBinsize_=32;
  nBinsizeSat_=64;
  nOldBinsize_=0;
}

Vector<int> IntegralHistogram::
GetHistogram(int ux,int uy,int lx, int ly)
{
  //D = D-C-B+A 
  Vector<int> histo = mHists_[ly-1][lx-1];
  if(ux>0 && uy>0)
  {
    histo -= mHists_[ly-1][ux-1]; //lowery - (upperx-1) 
    histo -= mHists_[uy-1][lx-1]; // (uppery-1) - lowerx 
    histo += mHists_[uy-1][ux-1]; // (uppery-1) - (uppery-1) 
  }
  else if(uy==0 && ux >0)
  {
      histo -= mHists_[ly-1][ux-1];  
  }
  else if(ux==0 && uy >0)
  {
      histo -= mHists_[uy-1][lx-1];
  }  
  return histo;
}


int IntegralHistogram::GenerateIntegralHist(Image<unsigned char>& img,
                                            int nBinsize)
{
  
  
  //at first test if image is in HSL
  if (img.GetColorModel()!=ImageBase::CM_HSL){
    BIASERR("GenerateIntegralHist() needs an image with HSL Colormodel!");
    return -1;
  }
  //only resize if changed;
  if(nBinsize != nOldBinsize_)
    {
      nBinsize_ = nBinsize;
      nOldBinsize_ = nBinsize;
    }
  nBinsizeSat_ = nBinsize;
  nOldBinsizeSat_ = nBinsize;

  const int nWidth = img.GetWidth();
  const int nHeight = img.GetHeight();
  unsigned char *pData = img.GetImageData(); //the image data
  const int nStep = 256/nBinsize_;  //the steps in pData
  const int nStep2=nStep*nStep;
  const int nNumIndex=nStep2+nStep;
  unsigned char hue=0, sat=0, light=0;
  int tmp = nBinsize_; //here calculate the shift
  nShift_=-1;
  while (tmp!=0) { 
    tmp=tmp>>1; nShift_++; 
  }


  //resize the hists array to prevent overflows and init to 0
  if(mHists_.size() != nWidth*nHeight){
    Vector<int> tmp(nNumIndex+1);
    mHists_.newsize(nHeight, nWidth);
    mHists_=tmp;
  }

 
  //now go through picture and generate the pixelvalues
  Vector<int> *uData, *lData, *RowEnd, *end;
  uData = lData = &(mHists_[0][0]);
  RowEnd = uData + nWidth;
  end = uData + nWidth*nHeight;
 
  //set  mHists_ [0][0]
  lData->Fill(0);
  hue = (*pData++)>>nShift_;
  sat = (*pData++);
  light = (*pData++)>>nShift_;
  if (sat>nMinSat_){
    sat = sat>>nShift_;
    (*lData)[hue*nStep+sat]++;
  } else {
    (*lData)[nStep2+light]++;
  }
  (*lData)[nNumIndex]=1; // set bincount at [0][0] = 1 (first element
  lData++;
  // row 0 (setting elements mHists_[0][x]
  while (lData < RowEnd){
    CopyNoCheck(lData[-1], *lData); // copy mHists_[0][0] to mHists_[0][1]
    hue = (*pData++)>>nShift_;
    sat = (*pData++);
    light = (*pData++)>>nShift_;
    if (sat>nMinSat_){
      sat = sat>>nShift_;
      (*lData)[hue*nStep+sat]++;
    } else {
      (*lData)[nStep2+light]++;
    }
    (*lData)[nNumIndex]++; //increase bincount for mHists_[0][x]
    lData++;
    
  }
  RowEnd += nWidth;
  //  usertime.Reset();
  //  usertime.Start();
  while (lData < end){
    // first column
    *lData = *uData;
    hue = (*pData++)>>nShift_;
    sat = (*pData++);
    light = (*pData++)>>nShift_;
    if (sat>nMinSat_){
    sat = sat>>nShift_;
      (*lData)[hue*nStep+sat]++;
    } else {
      (*lData)[nStep2+light]++;
    }
    (*lData)[nNumIndex]++;  //increase bincount 
    lData++;
    uData++;
    while (lData < RowEnd){
      *lData = *uData;
      *lData -= uData[-1];
      *lData += lData[-1];
      hue = (*pData++)>>nShift_;
      sat = (*pData++);
      light = (*pData++)>>nShift_;
      if (sat>nMinSat_){
        sat = sat>>nShift_;
        (*lData)[hue*nStep+sat]++;
      } else {
        (*lData)[nStep2+light]++;
      }
      (*lData)[nNumIndex]++;
      lData++;
      uData++;
    }
    RowEnd+=nWidth;
  }
  //  usertime.Stop();
  //  cout<<"Building the Int Histo takes \n:";
  //  usertime.PrintRealTime();
  // cout <<"usecs"<<endl;
  //  usertime.Reset();
  //Dump(cout, 239, 319);
  //assert(false);
  return 0;
  
}


int IntegralHistogram::GenerateIntegralHistDiffBin(Image<unsigned char>& img,
                                            int nBinsizeHue, int nBinsizeSat)
{
  // cout<<"CreateIntergralHistDiffBin binSizeHue :"<<nBinsizeHue<<endl;
  //cout<<"CreateIntergralHistDiffBin binSizeSat :"<<nBinsizeSat<<endl;
  
  //at first test if image is in HSL
  if (img.GetColorModel()!=ImageBase::CM_HSL){
    BIASERR("GenerateIntegralHist() needs an image with HSL Colormodel!");
    return -1;
  }
  //only resize if changed;
  if(nBinsizeHue != nOldBinsize_)
    {
      nBinsize_ = nBinsizeHue;
      nOldBinsize_ = nBinsizeHue;
    }
  if(nBinsizeSat != nOldBinsizeSat_)
    {
      nBinsizeSat_ = nBinsizeSat;
      nOldBinsizeSat_ = nBinsizeSat;
    }

  const int nWidth = img.GetWidth();
  const int nHeight = img.GetHeight();
  unsigned char *pData = img.GetImageData(); //the image data
  const int nStep = 256/nBinsize_;  //the steps in pData
  const int nStepSat = 256/nBinsizeSat_;
  const int nStep2Sat = nStepSat*nStep; 
  const int nNumIndex=nStep2Sat+nStep;
  unsigned char hue=0, sat=0, light=0;
  int tmp = nBinsize_; //here calculate the shift
  nShift_=-1;
  while (tmp!=0) { 
    tmp=tmp>>1; nShift_++; 
  }

  tmp = nBinsizeSat_;
  nShiftSat_=-1;
  while (tmp!=0) { 
    tmp=tmp>>1; nShiftSat_++;
  }
 


  //resize the hists array to prevent overflows and init to 0
  if(mHists_.size() != nWidth*nHeight){
    Vector<int> tmp(nNumIndex+1);
    mHists_.newsize(nHeight, nWidth);
    mHists_=tmp;
  }


  //now go through picture and generate the pixelvalues
  Vector<int> *uData, *lData, *RowEnd, *end;
  uData = lData = &(mHists_[0][0]);
  RowEnd = uData + nWidth;
  end = uData + nWidth*nHeight;
 
  //set  mHists_[0][0]
  lData->Fill(0);
  hue = (*pData++)>>nShift_;
  sat = (*pData++);
  light = (*pData++)>>nShift_;
  if (sat>nMinSat_){
    sat = sat>>nShiftSat_;
    (*lData)[hue*nStepSat+sat]++;
  } else {  
    (*lData)[nStep2Sat+light]++;
  }
  (*lData)[nNumIndex]=1; // set bincount at [0][0] = 1 (first element
  lData++;
  // row 0 (setting elements mHists_[0][x]
  while (lData < RowEnd){
    CopyNoCheck(lData[-1], *lData); // copy mHists_[0][0] to mHists_[0][1]
    hue = (*pData++)>>nShift_;
    sat = (*pData++);
    light = (*pData++)>>nShift_;
    if (sat>nMinSat_){
      sat = sat>>nShiftSat_; 
      (*lData)[hue*nStepSat+sat]++;
    } else {
      (*lData)[nStep2Sat+light]++;
    }
    (*lData)[nNumIndex]++; //increase bincount for mHists_[0][x]
    lData++;
    
  }
  RowEnd += nWidth;
  
  usertime.Reset();
  usertime.Start();
  while (lData < end){
    // first column
    *lData = *uData;
    hue = (*pData++)>>nShift_;
    sat = (*pData++);
    light = (*pData++)>>nShift_;
    if (sat>nMinSat_){
      sat = sat>>nShiftSat_;
      (*lData)[hue*nStepSat+sat]++;
    } else {
      (*lData)[nStep2Sat+light]++;
    }
    (*lData)[nNumIndex]++;  //increase bincount 
    lData++;
    uData++;
    while (lData < RowEnd){
      *lData = *uData;
      *lData -= uData[-1];
      *lData += lData[-1];
      hue = (*pData++)>>nShift_;
      sat = (*pData++);
      light = (*pData++)>>nShift_;
      if (sat>nMinSat_){  
        sat = sat>>nShiftSat_; 
    
        (*lData)[hue*nStepSat+sat]++;
      } else {
        (*lData)[nStep2Sat+light]++;
      }
      (*lData)[nNumIndex]++;
      lData++;
      uData++;
    }
    RowEnd+=nWidth;
  }
  usertime.Stop();
  //  cout<<"Building the Int Histo takes \n:";
  //usertime.PrintRealTime();
  // cout <<"usecs"<<endl;
  usertime.Reset();
  //Dump(cout, 239, 319);
  //assert(false);
  return 0;
  
}


double IntegralHistogram::
CalcSimilarity(ColorHistogram<unsigned char> &chist,
               int ux,int uy,int lx, int ly)
{
  double dist=0; //indicator for similarity

  //  BIAS::TimeMeasure time;
  //  time.Reset();
  //  time.Start();

  //D = D-C-B+A 
  Vector<int> histo = mHists_[ly-1][lx-1];
  if(ux>0 && uy>0)
    {
      histo -= mHists_[ly-1][ux-1]; //lowery - (upperx-1) 
      histo -= mHists_[uy-1][lx-1]; // (uppery-1) - lowerx 
      histo += mHists_[uy-1][ux-1]; // (uppery-1) - (uppery-1) 
    }
  else if(uy==0 && ux >0)
    {
      histo -= mHists_[ly-1][ux-1];  
    }
  else if(ux==0 && uy >0){
    histo -= mHists_[uy-1][lx-1];
  }  
  int nBinwidth =  256/nBinsize_; //width of one bin
  int nBinwidthSat= 256/nBinsizeSat_; // width of one bin for Sat
  //  int nBinwidthsq = nBinwidth*nBinwidth
  int nBinwidthsq = nBinwidth*nBinwidthSat;
  int nBinCount = histo[nBinwidthsq+nBinwidth];  //Number of bins
  //prevent division by o
  if(nBinCount == 0)
    nBinCount = 1;

  m_cHSData = chist.GetHist();
  m_cLData = chist.GetHistL();

  //save normalised values (divided by bincount) in normHisto
  Vector<double> normHisto(histo.size());
  for(int x=0;x<histo.size()-1;x++){
    normHisto[x]= (double)((double) histo[x] / nBinCount); 
  }
  //but leave last element as the bincount
  normHisto[histo.size()-1] = histo[histo.size()-1];
 
  //sum up the sqrt of the product of q*(n)and q(n;x) over all bins
  for(int i=0;i<nBinwidth;i++){
    //lValue
    dist += sqrt(((double)normHisto[nBinwidthsq+i])*m_cLData[i]);
    //hsValues
    //  for(int j=0;j<nBinwidth;j++){
    for(int j=0;j<nBinwidthSat;j++){
      dist += sqrt((normHisto[i*nBinwidthSat+j])*m_cHSData[i][j]);
    }
  }
  //new 
  if(!bExpDist_)
    dist = sqrt(1-dist);
  else 
    dist = exp(-nLambda_*(1-dist));  

  /* old
     if(bExpDist_) 
     dist = exp(-nLambda_*(1-dist));  
  */
  //  time.Stop();
  // cout<<"CalculateSimilarity took:"<<endl;
  // time.Print();
  return dist;
}

double 
IntegralHistogram::CalcSimilarity(ColorHistogram<unsigned char> &chist1,
                                  int ux1,int uy1,int lx1, int ly1,
                                  ColorHistogram<unsigned char> &chist2,
                                  int ux2,int uy2,int lx2, int ly2,
                                  double weightSecond)
{
  double sim1 = CalcSimilarity(chist1,ux1,uy1,lx1,ly1);
  double sim2 = CalcSimilarity(chist2,ux2,uy2,lx2,ly2);
  return (1.0-weightSecond)*sim1 + weightSecond*sim2;
}

/**
   dumps information on histogram to ostream
   format:
   binsize tab ux tab uy tab lx tab ly
   hs histogram, as a matrix with binsize x binsize (tab separated)
   l values as a line with lenght binsize (tab separated)
*/
void IntegralHistogram::Dump(std::ostream& os, int y, int x)
{
  os << " D = B + C - A\n";
  os << "A: "<<mHists_[y-1][x-1]<<endl
     << "B: "<<mHists_[y-1][x]<<endl
     << "C: "<<mHists_[y][x-1]<<endl
     << "D: "<<mHists_[y][x]<<endl;
  os << " D = B + C - A\n";
}