HistogramEqualization.cpp

#include "HistogramEqualization.hh"
#include <Image/HistogramImage.hh>
#include <Base/Image/ImageIO.hh>
#include <limits>

using namespace BIAS;
using namespace std;

template <class InputStorageType, class OutputStorageType>
int HistogramEqualization<InputStorageType,OutputStorageType>::
Filter(const Image<InputStorageType>& src,Image<OutputStorageType>& dst){

  if(src.GetChannelCount()!=1){
    BIASERR("HistogramEqualization: not for multiple channel images");
    return -1;
  }
  if (!dst.SamePixelAndChannelCount(src)){
    dst.Release();
    dst.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }

  //1. compute histogram
  HistogramImage hsrc;
  hsrc.AboveThresholdToValue(minCutOff_, minCutOff_);
  hsrc.BelowThresholdToValue(maxCutOff_, maxCutOff_);
  hsrc.Init(256,256,1);//we don't really need the drawing functions but we heave to initialize the class...
  unsigned bincount =(unsigned) (1.0 + std::numeric_limits<InputStorageType>::max() + 
                      fabs((double) std::numeric_limits<InputStorageType>::min()));
  hsrc.InitHist(bincount);
  hsrc.AddHist<InputStorageType>(src);
  
  unsigned binnum = hsrc.GetBinCount(); 
  BIASASSERT(binnum>0);
  unsigned pixnum = 0;
  map<unsigned,unsigned> filledbins;   
  //we allow only integer values. thus for every possible grey value a bin exists.   
  for(unsigned b=0; b<binnum; b++){
    unsigned val = (unsigned) hsrc.GetBin(b);
    if(val>0){
      filledbins[b]=val; 
      pixnum += val;
    }       
  }

  map<unsigned,unsigned>::iterator it;
  //unsigned numoffilledbins = 0; numoffilledbins = filledbins.size();
  map<InputStorageType,InputStorageType> NewGreyVal;   
  double sum = 0.0; 
  
  double maxval = (double)std::numeric_limits<InputStorageType>::max();
  double minval = (double)std::numeric_limits<InputStorageType>::min();
  double valrange = (maxval-minval);
  for(it=filledbins.begin(); it != filledbins.end(); it++){
    sum+= (double) it->second;
    InputStorageType oldgreyvalue = (InputStorageType)
                                    (it->first -  fabs((double) std::numeric_limits<InputStorageType>::min()));
    NewGreyVal[oldgreyvalue] = (InputStorageType)
                               (((sum/(double)pixnum)*valrange)+minval);
  }
  
  const InputStorageType** pDat = src.GetImageDataArray();
  OutputStorageType** pDstDat = dst.GetImageDataArray(); 
  const ROI *roi = src.GetROI();
  
  unsigned wid = src.GetWidth();
  unsigned hei = src.GetHeight();
  for(unsigned y=0; y<hei; y++)
  for(unsigned x=0; x<wid; x++){    
    if (roi->IsInROI(x,y)){
      InputStorageType oldgreyval = pDat[y][x];  
      //get new value:
      pDstDat[y][x]= (OutputStorageType) NewGreyVal[oldgreyval];
    }
  }

  return 0;
}


template <class InputStorageType, class OutputStorageType>
int HistogramEqualization<InputStorageType,OutputStorageType>::
AdaptiveHistogramEqualization(const Image<InputStorageType>& src,Image<OutputStorageType>& dst,
                              unsigned int numRows, unsigned int numCols){

  if(src.GetChannelCount()!=1){
    BIASERR("HistogramEqualization: not for multiple channel images");
    return -1;
  }
  if (!dst.SamePixelAndChannelCount(src)){
    dst.Release();
    dst.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }

  // init counts
  unsigned int imWidth = src.GetWidth();
  unsigned int imHeight = src.GetHeight();
  unsigned int numPxCellx = imWidth/numRows;
  unsigned int numPxCelly = imHeight/numCols;

  if(src.GetBitDepth() == 0){
    BIASERR("bit depth not set in image ");
    BIASABORT;
  }
  
  unsigned int histoSize = 1;
  for (unsigned int i = 0; i < src.GetBitDepth(); i++){
    histoSize *= 2;
  }
  
  // init histograms
  vector<vector<unsigned int> > histograms(numRows*numCols, vector<unsigned int>(histoSize, 0));
  vector<vector<unsigned int> > newGreyVals(numRows*numCols, vector<unsigned int>(histoSize, 0));
  vector<unsigned int> histoPixNum(numRows*numCols, 0);
  vector<Vector2<unsigned int> > cellCenters(numRows*numCols);
  
  const InputStorageType** srcPtr = src.GetImageDataArray();
  
  unsigned int startInCellX=0, startInCellY=0;
  unsigned int endInCellX=0, endInCellY=0;
  unsigned int histoCount=0;
  InputStorageType val=0;
  for(unsigned int cx=0; cx < numRows; cx++){
    for(unsigned int cy=0; cy < numCols; cy++){
      startInCellX = cx*numPxCellx;
      startInCellY = cy*numPxCelly;
      endInCellX = (cx+1)*numPxCellx;
      endInCellY = (cy+1)*numPxCelly;
      histoCount = numRows*cy+cx;
      
      cellCenters[histoCount] = Vector2<unsigned int>((endInCellX-numPxCellx/2),
                                                      (endInCellY-numPxCelly/2));
      
      for(unsigned int x=startInCellX; x < endInCellX; x++){
        for(unsigned int y=startInCellY; y < endInCellY; y++){
          val = srcPtr[y][x];
          if(val < minCutOff_) val = minCutOff_;
          if(val > maxCutOff_) val = maxCutOff_;
          histograms[histoCount][val]++;
        } // end celly
      } // end cellx
      
      // get pixel count in each cell - might not be equal to pixels in orig image due to
      // thresholding
      for(unsigned int hs=0; hs < histoSize; hs++){
        histoPixNum[histoCount] += histograms[histoCount][hs];
      }

      // compute new grey values
      double sum = 0.0; 
      double maxval = (double)std::numeric_limits<InputStorageType>::max();
      double minval = (double)std::numeric_limits<InputStorageType>::min();
      double valrange = (maxval-minval);
      for(unsigned int hs=0; hs < histoSize; hs++){
        sum+= (double) histograms[histoCount][hs];
        InputStorageType oldgreyvalue = (InputStorageType)
                                        (hs - fabs((double) std::numeric_limits<InputStorageType>::min()));
        newGreyVals[histoCount][oldgreyvalue] = (InputStorageType)
                                                (((sum/(double)histoPixNum[histoCount])*valrange)+minval);
      } // end for each histogram entry, compute new grey value
      
    } // end for all cell cols
  }// end for all cell rows
  
  // compute new grey values in output image
  OutputStorageType** pDstDat = dst.GetImageDataArray(); 
  const ROI *roi = src.GetROI();
 

  double w00=0, w01=0;
  double w10=0, w11=0;
  unsigned int cx1 = 0;
  unsigned int cx2 = 0;
  unsigned int cy1 = 0;
  unsigned int cy2 = 0;
  double distx1=0, distx2=0;
  double disty1=0, disty2=0;
  int xm1=0, xm2=0, ym1=0, ym2=0;
  for(unsigned y=0; y<imHeight; y++){
    for(unsigned x=0; x<imWidth; x++){    
      // determine nearby cells
      xm1 = ((int)x-numPxCellx/2);
      xm2 = ((int)x+numPxCellx/2);
      ym1 = ((int)y-numPxCelly/2);
      ym2 = ((int)y+numPxCelly/2);
      

      
      if(xm1 < 0) xm1=0;
      if(xm2 >= (int)imWidth) xm2 = imWidth-1;
      if(ym1 < 0) ym1=0;
      if(ym2 >= (int)imHeight) ym2 = imHeight-1;
      
      cx1 = xm1 / numPxCellx;
      cx2 = xm2 / numPxCellx;
      cy1 = ym1 / numPxCelly;
      cy2 = ym2 / numPxCelly;
      
      if(cx2 >= numRows) cx2 = numRows-1;
      if(cy2 >= numCols) cy2 = numCols-1;

      distx1 = fabs((double)x-(double)cellCenters[cy1*numRows+cx1][0]);
      distx2 = fabs((double)cellCenters[cy1*numRows+cx2][0]-(double)x);
      disty1 = fabs((double)y-(double)cellCenters[cy1*numRows+cx1][1]);
      disty2 = fabs((double)cellCenters[cy2*numRows+cx2][1]-(double)y);
      
      
      distx1 /= numPxCellx;
      distx2 /= numPxCellx;
      disty1 /= numPxCelly;
      disty2 /= numPxCelly;
      
      if(cx1==cx2){
        distx2 = 1.0 - distx1;
      }
      if(cy1==cy2){
        disty2 = 1.0 - disty1;
      }

      w00 = distx2*disty2;
      w01 = distx1*disty2;
      w10 = distx2*disty1;
      w11 = distx1*disty1;

      
      if (roi->IsInROI(x,y)){
        InputStorageType oldgreyval = srcPtr[y][x];  
        unsigned int oldIndex = (unsigned int)oldgreyval;
        //get new value:
//        cout << "dst acces " << dst.GetWidth() << " " << dst.GetHeight() << "  " << x << " " << y << endl;
//        cout << "newgreyvals " << newGreyVals.size() << " " << cy1*numRows+cx1
//             << " " << cy1*numRows+cx2 << " " << cy2*numRows+cx1 << " " << cy2*numRows+cx2
//             << " " << newGreyVals[cy1*numRows+cx1].size() << " " << newGreyVals[cy1*numRows+cx2].size()
//             << " " << newGreyVals[cy2*numRows+cx1].size() << " " << newGreyVals[cy2*numRows+cx2].size()
//             << " old index  " << oldIndex << endl;
//        cout << " num rows " << numRows << " cy2 " << cy2 << " cx1 " << cx1 << " cx2 " << cx2 << endl;
        pDstDat[y][x]= (OutputStorageType) (w00*newGreyVals[cy1*numRows+cx1][oldIndex] +
                                            w01*newGreyVals[cy1*numRows+cx2][oldIndex] +
                                            w10*newGreyVals[cy2*numRows+cx1][oldIndex] +
                                            w11*newGreyVals[cy2*numRows+cx2][oldIndex]);
      }
    }
  }

  return 0;
}


#define INST(intype,outtype) \
template class BIASImage_EXPORT HistogramEqualization<intype,outtype>;
namespace BIAS{
INST(char,char)
INST(unsigned char, unsigned char)
}
//won't work for float and double instances.
//even int and unsigned int aren't so good ideas.