ColorHistogram.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 "ColorHistogram.hh"
#include <Base/ImageUtils/BresenhamCircle.hh>
#include <Base/Image/ImageBase.hh>
#include <Base/Image/ImageConvert.hh>

#include <Base/Common/BIASpragma.hh>

using namespace BIAS;
using namespace std;

namespace BIAS {

const float maxShutter = 0.034f;
const float minShutter = 0.01f;
const float maxGain = 29.0f;
const float minGain = 0.1f;



/**
   This class creates a Histogram in the HSL color space.
   @author ingo schiller
   @date feb 2004
 */
template<class StorageType>
ColorHistogram<StorageType>::ColorHistogram()
{
  image_=NULL;
  histCount_ = 0;
  binCount_ = 0;
  expDist_=true;
  minSat_=10;
  lambda_=20;
  binSize_=32;
  binSizeSat_=64;
  ux_=uy_=lx_=ly_=0;
  SetBinSize(binSize_,binSizeSat_);

}

template<class StorageType>
ColorHistogram<StorageType>::~ColorHistogram()
{

}


//set the region of interest
template<class StorageType>
int ColorHistogram<StorageType>::SetROI(int ux, int uy, int lx, int ly){

  ux_=ux;
  uy_=uy;
  lx_=lx;
  ly_=ly;

    circle_=false;

  return 0;
}

//set the region of interest
template<class StorageType>
int ColorHistogram<StorageType>::SetROI(int x, int y, int r){

  mx_=x;
  my_=y;
  r_=r;

    circle_=true;

  return 0;
}

template<class StorageType>
int ColorHistogram<StorageType>::SetBinSize(int binsize){
  binSize_ = binsize;
  binSizeSat_=binsize;
  unsigned int nBins = (unsigned int)(256/binSize_);
  hist_.resize(nBins);
  for(unsigned int k=0;k<nBins;k++){
    hist_[k].resize(nBins);
  }
  histl_.resize(nBins);
#ifdef BIAS_DEBUG
  if (binSize_!=4 && binSize_!=8 && binSize_!=16 && binSize_!=32 && binSize_!=64 && binSize_!=128){
    BIASERR("slow, not supported");
    BIASABORT
  }
#endif
  int tmp = binSize_;
  shift_=-1;
  while (tmp!=0) {
    tmp=tmp>>1; shift_++;
  }

  shiftSat_ = shift_;
  //cerr << "binSize_: "<<binSize_<<" shift: "<<shift_<<endl;

  luminance_.resize(8); // declare fixed binsize for luminace of 32 -> 8 bins
  shiftLum_=5;
  return 0;
}


template<class StorageType>
int ColorHistogram<StorageType>::SetBinSize(int binSizeHue,
                                            int binSizeSat)
{
  binSize_ = binSizeHue;
  binSizeSat_ = binSizeSat;
  unsigned int nBinsHue = (unsigned int)(256/binSize_);
  unsigned int nBinsSat = (unsigned int)(256/binSizeSat_);

  hist_.resize(nBinsHue);
  for(unsigned int k=0;k<nBinsHue;k++){
    hist_[k].resize(nBinsSat);
  }
  histl_.resize(nBinsHue);
#ifdef BIAS_DEBUG
  if (binSize_!=4 && binSize_!=8 && binSize_!=16 && binSize_!=32 && binSize_!=64 && binSize_!=128){
    BIASERR("slow, not supported");
    BIASABORT;
  }
#endif
  int tmp = binSize_;
  shift_=-1;
  while (tmp!=0) {
    tmp=tmp>>1; shift_++;
  }
  //set shift for the sat binsize
  tmp = binSizeSat_;
  shiftSat_ = -1;
  while (tmp!=0) {
    tmp=tmp>>1; shiftSat_++;
  }

  //cerr << "binSize_: "<<binSize_<<" binSizeSat: "<<binSizeSat_<<endl;
  //cerr << "Shift: "<<shift_<<" shiftSat: "<<shiftSat_<<endl;

 luminance_.resize(8); // declare fixed binsize for luminace of 32 -> 8 bins
 shiftLum_=5;

 return 0;

}


template<class StorageType>
std::vector<std::vector<double> > ColorHistogram<StorageType>::GetHist() const
{
  return hist_;
}


template<class StorageType>
std::vector<double> ColorHistogram<StorageType>::GetHistL() const
{
  return histl_;
}

template<class StorageType>
int ColorHistogram<StorageType>::GetHeightROI() const {

  return ly_-uy_;
}

template<class StorageType>
int ColorHistogram<StorageType>::GetWidthROI() const {

  return lx_-ux_;
}


template<class StorageType>
void ColorHistogram<StorageType>::GetROI(std::vector<int>& roi) const
{
  roi.clear();
    if(circle_){
        roi.push_back(mx_);
        roi.push_back(my_);
        roi.push_back(r_);
    }else{
      roi.push_back(ux_);
    roi.push_back(uy_);
    roi.push_back(lx_);
    roi.push_back(ly_);
    }
}




template<class StorageType>
int ColorHistogram<StorageType>::Clear()
{
  expDist_=true;
  normalised_=false;
  binCount_=0;
  histCount_=0;
  ux_=0;
  uy_=0;
  lx_=0;
  ly_=0;
    mx_=0;
    my_=0;
    r_=0;
    circle_=false;

  return 0;
}

template<class StorageType>
int
ColorHistogram<StorageType>::GenerateHist(BIAS::Image<StorageType>&image,
                                          int ux, int uy, int lx, int ly)
{
  ux_=ux;
  uy_=uy;
  lx_=lx;
  ly_=ly;

    circle_=false;

  GenerateHist(image);
  return 0;
}

template<class StorageType>
int
ColorHistogram<StorageType>::GenerateCircleHist(BIAS::Image<StorageType>&image,
                                          int x, int y, int r){
  mx_=x;
  my_=y;
  r_=r;

    circle_ = true;

  GenerateCircleHist(image);
  return 0;
}



template<class StorageType>
int
ColorHistogram<StorageType>::GenerateHist(BIAS::Image<StorageType>& img)
{

  StorageType *data;

  if  ( (abs(ux_- lx_)==0) |  (abs(uy_-ly_)==0) )
    return 1;

  //at first test if image is in HSL
  if (img.GetColorModel()==ImageBase::CM_HSL){
    image_=&img;
  } else {
    BIASERR("GenerateHist() needs an image with HSL Colormodel!");
    return -1;
  }

  DeleteHist_();


  //go through the picturearray and collect the value and saturation
  //and set the l-values that are missing

  StorageType hue=0, sat=0, light=0;
  StorageType lum=0;
  data = image_->GetImageData();

  //go through imagedata and collect infos on colors and lightness
  int width = image_->GetWidth(); //width of image
  int step = 3*width-3*lx_+3*ux_; //right part left of image
  StorageType *p = data+3*ux_+3*width*uy_; //pointer to beginning of ROI
  StorageType *e = data+3*lx_+3*width*(ly_-1); //pointer to end of ROI
  StorageType *le = data+3*lx_+3*width*uy_; //pointer to line end


  //go through the picturearray and collect the value and saturation
  //and set the l-values that are missing

  double invbin=1.0/(double)binSize_;
  double invbinSat=1.0/(double)binSizeSat_;
  while(p<e){
    //cerr<<"ALERT!!! IN WHILE p<e!!"<<endl;
    while(p<le){
      hue = *p++; //get values and increase pointers
      sat = *p++;

      lum = *p;//take the lum but don't increase
      lum = lum*float(1.0/(double)32);//and then sort into luminacne vector
      luminance_[(int)floor((float)lum)]+=1;

      if(sat>minSat_)
        {
          //if SAT is >minSat_ put value in hist_
          sat = sat*float(invbinSat);
          hue = hue*float(invbin);
          hist_[(int)floor((float)hue)][(int)floor((float)sat)]+=1;
          binCount_++;//increase bincount for the normalisation
          p++;//increase pointer because it's not done in else case
        }
      else {
        //build the array for the L-values here
        light = *p++; //get value
        light = light*float(invbin);
        histl_[(int)floor((float)light)]+=1;
        //also increase bincount if in histl_
        binCount_++;

      }
    }
    le+=3*width;
    p+=step;
  }

  Normalise_();

  return 0;
}

template<>
int
ColorHistogram<unsigned char>::GenerateHist(Image<unsigned char>& img)
{
  unsigned char *data;
 // cout<<"ColorHistogram GenerateHist binsizeHue :"<<binSize_<<endl;
 // cout<<"ColorHistogram GenerateHist binsizeSat :"<<binSizeSat_<<endl;

  if ( (abs(ux_- lx_)==0) |  (abs(uy_-ly_)==0) )
    return 1;

  //at first test if image is in HSL
  if (img.GetColorModel()==ImageBase::CM_HSL){
    image_=&img;
  } else {
    BIASERR("GenerateHist() needs an image with HSL Colormodel!");
    return -1;
  }

  DeleteHist_();

  data = image_->GetImageData();

  //go through imagedata and collect infos on colors and lightness
  int width = image_->GetWidth(); //width of image
  int step = 3*width-3*lx_+3*ux_; //right part left of image
  unsigned char *p = data+3*ux_+3*width*uy_; //pointer to beginning of ROI
  unsigned char *e = data+3*lx_+3*width*(ly_-1); //pointer to end of ROI
  unsigned char *le = data+3*lx_+3*width*uy_; //pointer to line end


  //go through the picturearray and collect the value and saturation
  //and set the l-values that are missing

  unsigned char hue=0, sat=0, light=0;
  unsigned char lum=0;

  while(p<e){
    //cerr<<"ALERT!!! IN WHILE p<e!!"<<endl;
    while(p<le){
      hue = *p++; //get values and increase pointers
      sat = *p++;
      lum = *p;//take the lum but don't increase
      lum = lum>>shiftLum_;;//and then sort into luminacne vector
      luminance_[lum]+=1;

      if(sat>minSat_)
        {
          //if SAT is >minSat_ put value in hist_
          sat = sat>>shiftSat_;
          hue = hue>>shift_;
          hist_[hue][sat]+=1;
          binCount_++;//increase bincount for the normalisation
          p++;//increase pointer because it's not done in else case
        }
      else {
        //build the array for the L-values here
        light = *p++; //get value
        light = light>>shift_;
        histl_[light]+=1;
        //also increase bincount if in histl_
        binCount_++;

      }
    }
    le+=3*width;
    p+=step;
  }

  Normalise_();

  return 0;
}

template<class StorageType>
int
ColorHistogram<StorageType>::GenerateCircleHist(BIAS::Image<StorageType>& img){
  // Clear the histogram
  DeleteHist_();


    // Initialize the Bresenham Algorithm
    int height, width;
    int mp[2] = {mx_,my_};
  int next[2] = {0,0};
    BresenhamCircle circle = BresenhamCircle(mp,r_);

  //at first test if image is in HSL
  if (img.GetColorModel()==ImageBase::CM_HSL){
    image_=&img;
  } else {
    BIASERR("GenerateHist() needs an image with HSL Colormodel!");
    return -1;
  }

    // Initialize field for setting the histogram
  StorageType hue=0, sat=0, light=0;
  StorageType lum=0;
  StorageType **data = image_->GetImageDataArray();
    unsigned int cc = image_->GetChannelCount();
    height = image_->GetHeight();
    width = image_->GetWidth();


    // check, if the range of the circle is completly in the image area
    if( (mx_-r_<0) && (my_-r_<0) && (mx_+r_>width) && (my_+r_>height)){
    BIASERR ("The range of the circle needs to be completly in the image.");
        return -1;
  }



  //go through the picturearray and collect the value and saturation
  //and set the l-values that are missing

  double invbin=1.0/(double)binSize_;
  double invbinSat=1.0/(double)binSizeSat_;

  while (circle.GetNext(next)){
        lum = data[next[1]][next[0]*cc+2];
        lum = lum*float(1.0/(double)32);//and then sort into luminacne vector
        luminance_[(int)floor((float)lum)]+=1;
        if (sat>minSat_){
      hue = data[next[1]][next[0]*cc+0];
      sat = data[next[1]][next[0]*cc+1];
            hue = hue*float(invbin);
            sat = sat*float(invbinSat);
            hist_[(int)floor((float)hue)][(int)floor((float)sat)]+=1;
            binCount_++;//increase bincount for the normalisation
        }else{
        //build the array for the L-values here
            light = data[next[1]][next[0]*cc+2]; //get value
            light = light*float(invbin);
            histl_[(int)floor((float)light)]+=1;
            binCount_++;//also increase bincount if in histl_
        }
  }

  Normalise_();

  return 0;
}

template<>
int
ColorHistogram<unsigned char>::GenerateCircleHist(Image<unsigned char>& img){
  // Clear the histogram
  DeleteHist_();

  // Initialize the Bresenham Algorithm
  int height, width;
  int mp[2] = {mx_,my_};
  int next[2] = {0,0};
  BresenhamCircle circle = BresenhamCircle(mp,r_);

  //at first test if image is in HSL
  if (img.GetColorModel()==ImageBase::CM_HSL){
    image_=&img;
  } else {
    BIASERR("GenerateHist() needs an image with HSL Colormodel!");
    return -1;
  }

  // Initialize field for setting the histogram
  unsigned char hue=0, sat=0, light=0;
  unsigned char lum=0;
  unsigned char **data = image_->GetImageDataArray();
  unsigned int cc = image_->GetChannelCount();
  height = image_->GetHeight();
  width = image_->GetWidth();


  // check, if the range of the circle is completly in the image area
  if( (mx_-r_<0) && (my_-r_<0) && (mx_+r_>width) && (my_+r_>height)){
    BIASERR ("The range of the circle needs to be completly in the image.");
    return -1;
  }

  //go through the picturearray and collect the value and saturation
  //and set the l-values that are missing

  double invbin=1.0/(double)binSize_;
  double invbinSat=1.0/(double)binSizeSat_;

  while (circle.GetNext(next)){
    lum = data[next[1]][next[0]*cc+2];
    //and then sort into luminacne vector
    lum = (unsigned char)(float(lum)*(1.0f/32.0f));
    luminance_[(int)floor((float)lum)]+=1;
    if (sat>minSat_){
      hue = data[next[1]][next[0]*cc+0];
      sat = data[next[1]][next[0]*cc+1];
      hue = (unsigned char)(double(hue)*invbin);
      sat = (unsigned char)(double(sat)*invbinSat);
      hist_[(int)floor((float)hue)][(int)floor((float)sat)]+=1;
      binCount_++;//increase bincount for the normalisation
    }else{
      //build the array for the L-values here
      light = data[next[1]][next[0]*cc+2]; //get value
      light = (unsigned char)(double(light)*invbin);
      histl_[(int)floor((float)light)]+=1;
      binCount_++;//also increase bincount if in histl_
    }
  }

  Normalise_();

  return 0;
}


template<class StorageType>
double ColorHistogram<StorageType>::CalcSimilarity(const BIAS::ColorHistogram<StorageType> chist) const
{
  double dist=0; //indicator for similarity
  int nBins = 256/binSize_;
  int nBinsSat = 256/binSizeSat_;
  //sum up the sqrt of the product of q*(n)and q(n;x) over all bins
  for(int i=0;i<nBins;i++){
    dist += sqrt(histl_[i]*chist.histl_[i]);
    for(int j=0;j<nBinsSat;j++){
      dist += sqrt(hist_[i][j]*chist.hist_[i][j]);
    }
  }
  //new
  if(!expDist_)
    dist = sqrt(1-dist);
  else
    dist = exp(-lambda_*(1-dist));

  /* old
  if(expDist_)
    dist = exp(-lambda_*(1-dist));
  */
  //  cout<<"CalcSimilarity, dist: "<<dist<<endl;
  return dist;
}

template<class StorageType>
double ColorHistogram<StorageType>::CalcSimilarity(const BIAS::ColorHistogram<StorageType> chist, double &weight) const
{
  double dist=0; //indicator for similarity
  weight=0;
  int nBins = 256/binSize_;
  int nBinsSat = 256/binSizeSat_;
  //sum up the sqrt of the product of q*(n)and q(n;x) over all bins
  for(int i=0;i<nBins;i++){
    if (histl_[i]!=0) {
      weight++;
      if (chist.histl_[i]!=0)
        dist += sqrt(histl_[i]*chist.histl_[i]);
    }
    for(int j=0;j<nBinsSat;j++){
      if (hist_[i][j]!=0) {
        weight++;
        if (chist.hist_[i][j]!=0)
          dist += sqrt(hist_[i][j]*chist.hist_[i][j]);
      }
    }
  }

 //new
  if(!expDist_)
    dist = sqrt(1-dist);
  else
    dist = exp(-lambda_*(1-dist));

  /* old
  if(expDist_)
    dist = exp(-lambda_*(1-dist));
  */


  weight /= (256.0/(double)binSize_) * (256.0/(double)binSize_);
  return dist;
}

template<class StorageType>
int
ColorHistogram<StorageType>::Update(const ColorHistogram<StorageType> &hist,
                                    double factor)
{
#ifdef BIAS_DEBUG
  if (hist.binSize_!=binSize_) {
    BIASERR("can't update , beause histogramms have different binSize");
    return -1;
  }
#endif

  //update the size
  ux_= int((double)ux_*(1.0-factor)+factor*(double)hist.ux_ + 0.5);
  uy_= int((double)uy_*(1.0-factor)+factor*(double)hist.uy_ + 0.5);
  lx_= int((double)lx_*(1.0-factor)+factor*(double)hist.lx_ + 0.5);
  ly_= int((double)ly_*(1.0-factor)+factor*(double)hist.ly_ + 0.5);

  int nBins = 256/binSize_;
  int nBinsSat = 256/binSizeSat_;

  // update bins
  for(int i=0;i<nBins;i++){
    histl_[i]= histl_[i]*(1.0-factor) + factor * hist.histl_[i];
    for(int j=0;j<nBinsSat;j++){
      hist_[i][j] = hist_[i][j]*(1.0-factor) + factor * hist.hist_[i][j];
    }
  }
  return 0;
}


template<class StorageType> int ColorHistogram<StorageType>::
CheckImageQualityRef(ColorHistogram<StorageType>& cHist,
             float percent,
             float &shutter,
             float &gain)
{
  int ret =0;
  float shutterstep = 0.001f;
  float gainstep = 0.1f;
  double lowerLum =luminance_[0]+luminance_[1];
  double higherLum = luminance_[6]+luminance_[7];

  double lowerLumRef = cHist.luminance_[0]+cHist.luminance_[1];
  double higherLumRef = cHist.luminance_[6]+cHist.luminance_[7];
  //  cout<<"LowerLum :"<<lowerLum<<" LowerLumRef :"<<lowerLumRef<<endl;
  // cout<<"HigherLum :"<<higherLum<<" HigherLumRef :"<<higherLumRef<<endl;

  if((lowerLum - lowerLumRef) > percent)
    {
      if(shutter + shutterstep < maxShutter)
    {
      shutter += shutterstep;
      //  cout<<"Increasing shutter, lum too low"<<endl;
    }
      else
    {
      if(gain + gainstep < maxGain )
        {
          gain += gainstep;
          //  cout<<"Increasing gain, lum too low and shutter already max"<<endl;
        }
      else{}
        //cout<<"All max it is to dark to play outside"<<endl;
    }
      ret = 1;
    }
  else if((higherLum - higherLumRef) > percent)
    {
      if(gain - gainstep > minGain + gainstep)
    {
      gain -= gainstep;
      //cout<<"Decreasing gain, lum too high and gain not min"<<endl;
    }
      else if(shutter - shutterstep > minShutter)
    {
      shutter -= shutterstep;
      //  cout<<"Decreasing shutter, lum too high"<<endl;
    }
      else
    {
      //cout<<"All min, high risk of sunburn"<<endl;
    }
      ret = -1;
    }
  else {
    //  cout<<"Light is good, no need to change"<<endl;
  }

  return ret;
}




template<class StorageType>
int ColorHistogram<StorageType>::CheckImageQualityComplete(float percent,
                               float &shutter,
                               float &gain)
{
  double lowerLum =luminance_[0]+luminance_[1];
  double higherLum = luminance_[6]+luminance_[7];
  int ret =0;
  if(lowerLum > percent)
    {
      //cout<<"Luminance is low: increase!"<<endl;
      shutter = 1;
      ret = 1;
    }
  else if(higherLum > percent)
    {
      //cout<<"Luminance is high: decrease!"<<endl;
      shutter = -1;
      ret = -1;
    }
  else
    {
      //cout<<"Luminace is good, no need to modify"<<endl;
    }
  return ret;
}


template<class StorageType>
int ColorHistogram<StorageType>::CheckImageQuality(float percent,
                           float &shutter,
                           float &gain)
{
  std::vector<double> vHistL = GetHistL();

     if(GetBinSize() == 32)
       {

         double lowerLum =vHistL[0]+vHistL[1]+vHistL[2]+vHistL[3];
         double higherLum = vHistL[5]+vHistL[6]+vHistL[7];

         if(lowerLum > percent)
           {
             cout<<"Luminance is low: increase!"<<endl;
             shutter = 1;
         return 1;
           }
         else if(higherLum > percent)
          {
             cout<<"Luminance is high: decrease!"<<endl;
             shutter = -1;
         return -1;
           }
     else
       {
         cout<<"Luminace is good, no need to modify"<<endl;
         return 0;
       }
       }
     else
       {
         BIASERR("Only implemented for binSizeHue==8");
     return -1;
       }

  //return 0;
}//end CheckImageQuality

/**
   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)
*/
template<class StorageType>
void ColorHistogram<StorageType>::DumpCSV(std::ostream& os)
{
  int nBins = 256/binSize_;
  int nBinsSat = 256/binSizeSat_;

  for(int i=0;i<nBins;i++){
    for(int j=0;j<nBinsSat;j++){
     os <<hist_[i][j]<<",";
    }
    os <<endl;
  }
  os<<endl; //create empty line
  //luminance
  for(int i=0;i<nBins;i++){
    os <<histl_[i]<<",";
  }

  //complete luminance
  os <<endl<<endl;
  for(int i=0;i<8;i++){
    os <<luminance_[i]<<",";
  }
}

/**
   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)
*/
template<class StorageType>
void ColorHistogram<StorageType>::Dump(std::ostream& os)
{
  int nBins = 256/binSize_;
  int nBinsSat = 256/binSizeSat_;

    if(circle_){
    os<<binSize_<<"\t"<<binSizeSat_<<"\t"<<mx_<<
      "\t"<<my_<<"\t"<<r_<<"\t"<<(-1)<<endl;
    }else{
    os<<binSize_<<"\t"<<binSizeSat_<<"\t"<<ux_<<
      "\t"<<uy_<<"\t"<<lx_<<"\t"<<ly_<<endl;
  }

  for(int i=0;i<nBins;i++){
    for(int j=0;j<nBinsSat;j++){
     os <<hist_[i][j]<<"\t";
    }
    os <<endl;
  }
  os<<endl; //create empty line

  for(int i=0;i<nBins;i++){
    os <<histl_[i]<<"\t";
  }
  //create empty line
  os <<endl<<endl;
   for(int i=0;i<8;i++){
    os <<luminance_[i]<<"\t";
  }

  os <<endl<<endl;
  os <<binCount_<<endl;
}

template<class StorageType>
int ColorHistogram<StorageType>::GetNumberOfEntries()
{

  return 6 + (256/binSize_)*(256/binSizeSat_) + 2*(256/binSize_) +1;
}




//reads the information of a histogram from disk
template<class StorageType>
int ColorHistogram<StorageType>::Read(const char *filename)
{
  cout<<"Reading the histogram"<<endl;
  ifstream histofile(filename);
  if (!histofile){
    cerr << "unable to open " << filename<< endl;
    perror("");
    return -1;
  }
  int binSize,binSizeSat,ux,uy,lx,ly;

  histofile >> binSize;
  histofile >> binSizeSat;
  histofile >> ux >> uy >> lx >> ly;
  SetBinSize(binSize,binSizeSat);
    if(ly <0){
        SetROI(ux,uy,lx);
    }else{
    SetROI(ux,uy,lx,ly);
    }
  int nBins = 256/binSize_;
  int nBinsSat = 256/binSizeSat_;
  for(int i=0;i<nBins;i++){
    for(int j=0;j<nBinsSat;j++){
      histofile >> hist_[i][j];
    }
  }

  //And now the L-vector
  for(int i=0;i<nBins;i++){
    histofile >> histl_[i];

  }

 for(int i=0;i<8;i++){
    histofile >> luminance_[i];

  }


  return 0;
}


template<class StorageType>
void ColorHistogram<StorageType>::Normalise_()
{
  //go through whole vector and divide by binCount_
  const int end=256>>shift_;
  const int endSat=256>>shiftSat_;
  const int endLum=256>>shiftLum_;


  for(int i=0;i<end;i++){
    histl_[i]/=binCount_;
    for(int j=0;j<endSat;j++){
     hist_[i][j] /=binCount_;
    }
  }
  //normalise luminance as well
  for(int k=0;k<endLum;k++)
    {
      luminance_[k]/=binCount_;
    }

  normalised_=true;
}


template<class StorageType>
int ColorHistogram<StorageType>::DeleteHist_()
{
  binCount_=0;
  histCount_=0;
  histl_.assign(histl_.size(),0);
  for (unsigned int i=0;i<hist_.size();i++) {
    hist_[i].assign(hist_.size(),0);
  }

  return 0;
}

template<class StorageType>
float ColorHistogram<StorageType>::GetSaturatedBins() const
{
  const int end=256>>shift_;
  const int endSat=256>>shiftSat_;
  float count=0;
  for(int i=0;i<end;i++){
    for(int j=0;j<endSat;j++){
      if (hist_[i][j]!=0) count++;
    }
  }
  count/=float(end*endSat);
  return count;
}


template<class StorageType>
char* ColorHistogram<StorageType>::GetBinaryData(int &length){
  int dataLength= sizeof(double) * (hist_.size()*hist_[0].size() + histl_.size()*2)
    +  7*sizeof(int);

  char *ptr = new char[dataLength];
  int *tmp_int_ptr=(int *)ptr;
  *tmp_int_ptr++ = binSize_;
  *tmp_int_ptr++ = binSizeSat_;
  *tmp_int_ptr++ = ux_;
  *tmp_int_ptr++ = uy_;
  *tmp_int_ptr++ = lx_;
  *tmp_int_ptr++ = ly_;
  *tmp_int_ptr++ = binCount_;

  length += 7*sizeof(int);
  int numSat = 256/binSizeSat_;
  int numHue = 256/binSize_;


  double *tmp_ptr=(double*)tmp_int_ptr;
  for (int h=0; h<numHue; h++) {
    for (int s=0; s<numSat; s++){
      *tmp_ptr++ = hist_[h][s];
      length+=sizeof(double);
    }
  }
  for(int l=0;l<numHue;l++){
    *tmp_ptr++ = histl_[l];  length+=sizeof(double);}


  for(int t=0;t<numHue;t++){
    *tmp_ptr++ = luminance_[t];  length+=sizeof(double);}

  return ptr;
}

template<class StorageType>
void ColorHistogram<StorageType>::ReadBinaryData(char *data)
{

  int *tmp_int_ptr =(int *)data;
  binSize_ = *tmp_int_ptr++;
  binSizeSat_ = *tmp_int_ptr++;
  ux_ = *tmp_int_ptr++;
  uy_ = *tmp_int_ptr++;
  lx_ = *tmp_int_ptr++;
  ly_ = *tmp_int_ptr++;
  binCount_ = *tmp_int_ptr++;

  double *tmp_ptr=(double*)tmp_int_ptr;

  SetBinSize(binSize_,binSizeSat_);
  SetROI(ux_,uy_,lx_,ly_);
  int nBins = 256/binSize_;
  int nBinsSat = 256/binSizeSat_;
  for(int i=0;i<nBins;i++){
    for(int j=0;j<nBinsSat;j++){
      hist_[i][j] = *tmp_ptr++;
    }
  }

  //And now the L-vector
  for(int i=0;i<nBins;i++){
    histl_[i] = *tmp_ptr++;
  }

 for(int i=0;i<nBins;i++){
   luminance_[i] = *tmp_ptr++;
 }

}


#define INST(type) template class ColorHistogram<type>;

INST(float)
INST(unsigned char)

  // disabled untested instances, myriards of warnings

// #ifdef BUILD_IMAGE_INT
// INST(int);
// #endif
// #ifdef BUILD_IMAGE_CHAR
// INST(char);
// #endif
// #ifdef BUILD_IMAGE_SHORT
// INST(short);
// #endif
// #ifdef BUILD_IMAGE_USHORT
// INST(unsigned short);
// #endif
// #ifdef BUILD_IMAGE_USHORT
// INST(unsigned int);
// #endif
// #ifdef BUILD_IMAGE_DOUBLE
// INST(double);
// #endif

} // namespace BIAS