RasterBand.hpp

Go to the documentation of this file.

#ifndef LGL_RASTERBAND_HPP
#define LGL_RASTERBAND_HPP

#include <limits>
#include <iostream>
#include <boost/static_assert.hpp>

#include <gdal_priv.h>

#include <lgl/Raster.hpp>

namespace jafar { // namespace jafar
  namespace lgl { // namespace lgl
    
    template<typename T>
    struct DataType {
      static const bool value=false;
    };
    
    template<>
    struct DataType<double> {
      static const bool value=true;
      static const GDALDataType dataType = GDT_Float64;
      static const double min() { return 0.0; };
      static const double max() { return 1.0; };
      /*
      const static double MIN = 0.0;
      const static double MAX = 1.0;        
      */
    };

    template<>
    struct DataType<float> {
      static const bool value=true;
      static const GDALDataType dataType = GDT_Float32;
      static const double min() { return 0.0; };
      static const double max() { return 1.0; };
      /*
      const static double MIN = 0.0;
      const static double MAX = 1.0;        
      */
    };
    
    template<>
    struct DataType<unsigned char> {
      static const bool value=true;
      static const GDALDataType dataType = GDT_Byte;
      static const unsigned char min() { return 0;  };
      static const unsigned char max() { return 255;};        
    };
    
    template<>
    struct DataType<unsigned int> {
      static const bool value=true;
      static const GDALDataType dataType = GDT_UInt32;
      static const int min() { return 0;};
      static const int max() { return UINT_MAX;};   
    };
    
    enum FUSION_METHOD {MAX, MIN, MEAN, BAYES};
    enum HOMOGENEOUS_TYPE {HT_NONE=-1, SIMPLE, NEARLY, SIGMA};
      
    template <typename T>
    class GDRasterBand : public Raster { //{{{
        
      public:
        GDRasterBand(GDALDataset* _dataset, GDALRasterBand *_band, 
              double _weight=1.0,
              double _decompositionLimit=8.0,
              FUSION_METHOD _fmethod=MEAN,
              HOMOGENEOUS_TYPE _homogType=SIGMA);
        
        ~GDRasterBand();
        
        void setGDALRasterBand(GDALRasterBand* _band);      

        //-------- Getters --------//
        GDALRasterBand* getGDALRasterBand();

        double getWeight();
        
        int getTypeSize()const;
        GDALDataType getType()const;
              
        int getxroot() const;     
        int getyroot() const;     

        int getxsize() const { return xsize; }
        int getysize() const { return ysize; };

        double getxsize_m() const;      
        double getysize_m() const;      

        double getxsize_m(const RasterRect& rect) const;      
        double getysize_m(const RasterRect& rect) const;      

        double getScale() const;

        void reset(T _v); 

        // a raster object should be able to return the data for a cluster
        // or a raster cell
        //T getData(const RasterCellIndex& cell) const;   
        //T getData(const RasterRect& rect) const;  
        T getValue(const RasterPoint& rpt) const; 

        T getValue(int i) const;  

        T* getValuePtr() const; 

        bool getValue(const int& x, const int& y, void* value) const;

        T get(const RasterCellIndex& cell) const;

        bool get(const RasterRect& rect, std::vector<T>& v) const;

        //void* getBlock(const RasterRect&) const;
        double getDecompositionLimit();
        
        //-------- Setters --------//

        void setValue(const RasterCellIndex& cell, T value) const;
        void setValue(const int index, T value) const;

        void setScale(double _scale);

        void setDecompositionLimit(double _dl);

        //-------- Calculation/Operations --------//

        bool sameValue(const RasterRect& rect) const;     

        bool nearlySameValue(const RasterRect& rect) const;

        bool sigmaSameValue(const RasterRect& rect) const;

        void setProbabilityLag(T newProbaLag);

        void setSigmaFactor(T newSigmaFactor);

        T getProbabilityLag() const;

        T getSigmaFactor() const;

        double mean(const RasterRect& rect) const;      
        double sd(const RasterRect& rect) const;
        
        // given a rect ((top,left),(right,bottom)), tell if this is homogeneous
        virtual bool isHomogeneous(const RasterRect& rect) const=0;

        HOMOGENEOUS_TYPE getHomogeneousType() const;

        void setHomogeneousType(HOMOGENEOUS_TYPE hhh); 

        bool  bufferToMem();

        bool  memToBuffer();

      protected:
        GDALDataset *dataset;

        int xsize;

        int ysize;
        
        double absWEResolution;

        double absNSResolution;

        GDALRasterBand* band;

        GDALDataType dataType;

        double weight;

        double decompositionLimit;

        FUSION_METHOD fmethod;
        
        T* values;

        T probabilityLag;

        T sigmaFactor;

        HOMOGENEOUS_TYPE homogType;

      private :
        BOOST_STATIC_ASSERT(DataType<T>::value);

    }; //}}} end of class GDRasterBand 


     template<typename T>
     GDRasterBand<T>::GDRasterBand(GDALDataset* _dataset, GDALRasterBand *_band, double _weight, double _decompositionLimit, FUSION_METHOD _fmethod, HOMOGENEOUS_TYPE _homogType):
      dataset(_dataset),
      band(_band),
      weight(_weight),
      decompositionLimit(_decompositionLimit),
      fmethod(_fmethod),
      homogType(_homogType)
    { //{{{
      
      if ((dataset==NULL)||(band==NULL)) {
        std::cout << "#WWW# [RasterBand] Create raster band with NULL dataset and gdal band " << std::endl;
        return;
      }     
      
      //std::cout << "#III# [RasterBand] Load raster band of type : " << GDALGetDataTypeName(band->GetRasterDataType()) << std::endl;
      
      if ((band->GetRasterDataType()!=GDT_Byte)
          &&(band->GetRasterDataType()!=GDT_UInt16)) {
        std::cout << "#EEE# [RasterBand] Error, the band type should be Byte or Uint16 " << std::endl;
        return;
      }

      // NULL pointing values array
      values = NULL;

      // Nb of columns and lines
      xsize=dataset->GetRasterXSize();
      ysize=dataset->GetRasterYSize();

      // Getting the resolution
      double thisGeoTransform[6];
      dataset->GetGeoTransform( thisGeoTransform );
      absWEResolution = fabs(thisGeoTransform[1]);
      absNSResolution = fabs(thisGeoTransform[5]);

      // buffer the data to memory and make the conversion from input
      // type (uint8 or uint16) to double value between [0,1]
      bufferToMem();

      dataType = GDT_Float64;

      setProbabilityLag((T)0.05);
      setSigmaFactor((T)1.5);
    } //}}}
     
    template<typename T>
    GDRasterBand<T>::~GDRasterBand() {
      free(values);
    }
    
    template<typename T>
    bool GDRasterBand<T>::bufferToMem() 
    {//{{{
      if (values != NULL)
      {
        free(values);
      }
      values = (T*)malloc(sizeof(T)*xsize*ysize);  
      
      // convert all to double
      double* tmp = (double*)malloc(sizeof(double)*xsize*ysize);    
      band->RasterIO( GF_Read, 0, 0, xsize , ysize, tmp, xsize, ysize, GDT_Float64, 0, 0 );
      
      T inputMin;
      T inputMax;
      if (band->GetRasterDataType()==GDT_Byte) {
        inputMin = (T)(std::numeric_limits<unsigned char>::min());
        inputMax = (T)(std::numeric_limits<unsigned char>::max());
      } else if (band->GetRasterDataType()==GDT_UInt16) {
        inputMin = (T)(std::numeric_limits<unsigned short>::min());
        inputMax = (T)(std::numeric_limits<unsigned short>::max());
      } else {
        return false;
      }
      
      T constFactor = ((DataType<T>::max() - DataType<T>::min())/(inputMax-inputMin));
      for (int j=0;j<ysize;j++) {
        for (int i=0;i<xsize;i++) {
          // saving data line by line (X axis first)
          values[i+xsize*j] = (T)((tmp[i+xsize*j]-inputMin)*constFactor);
        }
      }
      
      free(tmp);
      return true;
    }//}}}
    
    template<typename T>
    bool GDRasterBand<T>::memToBuffer()
    { //{{{
      // convert all to double
      double* tmp = (double*)malloc(sizeof(double)*xsize*ysize);    

      T inputMin;
      T inputMax;
      if (band->GetRasterDataType()==GDT_Byte) {
        inputMin = (T)(std::numeric_limits<unsigned char>::min());
        inputMax = (T)(std::numeric_limits<unsigned char>::max());
      } else if (band->GetRasterDataType()==GDT_UInt16) {
        inputMin = (T)(std::numeric_limits<unsigned short>::min());
        inputMax = (T)(std::numeric_limits<unsigned short>::max());
      } else {
        return false;
      }

      for (int j=0;j<ysize;j++) {
        for (int i=0;i<xsize;i++) {
          // Buffering data line by line (X axis first)
          tmp[i+xsize*j] = (T)(values[i+xsize*j]*((inputMax-inputMin)/(DataType<T>::max() - DataType<T>::min())));
        }
      }

      band->RasterIO( GF_Write, 0, 0, xsize , ysize, tmp, xsize, ysize, GDT_Float64, 0, 0 );    

      free(tmp);

      return true;
    } //}}}

    template<typename T>
    void GDRasterBand<T>::setGDALRasterBand(GDALRasterBand* _band) {
      band=_band;  
    }

    template<typename T>
    GDALRasterBand* GDRasterBand<T>::getGDALRasterBand() {
      return band;
    }


    template<typename T>
    int GDRasterBand<T>::getxroot() const {
      return 0;
    }

    template<typename T>
    int GDRasterBand<T>::getyroot() const {
      return 0;
    }

    template<typename T>
    double GDRasterBand<T>::getxsize_m() const {
      return ((double)xsize)*absWEResolution;
    }

    template<typename T>
    double GDRasterBand<T>::getysize_m() const {
      return ((double)ysize)*absNSResolution;
    }

    template<typename T>
    double GDRasterBand<T>::getxsize_m(const RasterRect& rect) const {
      return ((double)(rect.xsize()))*absWEResolution;
    }

    template<typename T>
    double GDRasterBand<T>::getysize_m(const RasterRect& rect) const {
      return ((double)(rect.ysize()))*absNSResolution;
    }

    template<typename T>
    double GDRasterBand<T>::getScale() const {
      return (double)band->GetScale();
    }

    template<typename T>
    void GDRasterBand<T>::setScale(double _scale) {
      band->SetScale((double)_scale);
    }

    // get the weight of this band  
    template<typename T>
    double GDRasterBand<T>::getWeight() {
      return weight;
    }

    template<typename T>
    int GDRasterBand<T>::getTypeSize() const {
      return GDALGetDataTypeSize(dataType);
    }
   
    template<typename T>
    GDALDataType GDRasterBand<T>::getType() const {
      return dataType;
    }
     
    template<typename T>
    T GDRasterBand<T>::get(const RasterCellIndex& cell) const {
      return values[cell.i+xsize*cell.j];
    }

    template<typename T>
    void GDRasterBand<T>::setValue(const RasterCellIndex& cell, T value) const {
      values[cell.i+xsize*cell.j]=value;
    }

    template<typename T>
    void GDRasterBand<T>::setValue(const int index, T value) const {
      values[index]=value;
    }

    template<typename T>
    T* GDRasterBand<T>::getValuePtr() const {
      return values;
    }

    template<typename T>
    T GDRasterBand<T>::getValue(const RasterPoint& rpt) const {
      return values[rpt.x+xsize*rpt.y];
    }

    template<typename T>
    T GDRasterBand<T>::getValue(const int i) const {
      return values[i];
    }

    template<typename T>
    bool GDRasterBand<T>::getValue(const int& x, const int& y, void* value) const 
    {//{{{
      // read at (x,y) pixel with RasterIO
      // automatic conversion from unsigned byte to double 
      if ( band->RasterIO( GF_Read, x, y, 1, 1, value, 1, 1, dataType, 0, 0 ) == CE_Failure )
      {
        return false;
      }
      return true;
    }//}}}

    template<typename T>
    void GDRasterBand<T>::reset(T _v)
    {//{{{
      for (int j=0;j<ysize;j++) {
        for (int i=0;i<xsize;i++) {
          // saving data line by line (X axis first)
          values[i+xsize*j] = _v;
        }
      }
    }//}}}

    template<typename T>
    bool GDRasterBand<T>::get(const RasterRect& rect, std::vector<T>& v) const
    {//{{{
      int i,j;
      // modify vector capacity 
      v.reserve(rect.xsize()*rect.ysize());
      // Fill vector with raster values in the ROI specified by rect
      for (j = rect.ytopleft; j <= rect.yend(); j++) {
        for (i = rect.xtopleft; i <= rect.xend(); i++) {
          v.push_back(values[i+xsize*j]);
        }
      }
      return true;
    }//}}}
   
    template<typename T>
    bool GDRasterBand<T>::sameValue(const RasterRect& rect) const
    {//{{{  
      int i,j;
      T previous_value, cur_value = 0;
      previous_value = get(RasterCellIndex(rect.xtopleft,rect.ytopleft));
      for (i=rect.xtopleft; i < rect.xtopleft+rect.xsize(); i++) {
      for (j=rect.ytopleft; j < rect.ytopleft+rect.ysize(); j++) {
          cur_value = get(RasterCellIndex(i,j));
          if (cur_value != previous_value) {
            return false;
          }
          previous_value=cur_value;
        }
      }
      return true;
    }//}}}

    template<typename T>
    T GDRasterBand<T>::getProbabilityLag() const
    {//{{{
      return probabilityLag;
    }//}}}

    template<typename T>
    T GDRasterBand<T>::getSigmaFactor() const
    {//{{{
      return sigmaFactor;
    }//}}}

    template<typename T>
    void GDRasterBand<T>::setHomogeneousType(HOMOGENEOUS_TYPE hhh) 
    {//{{{
      homogType = hhh;
    }//}}}

    template<typename T>
    HOMOGENEOUS_TYPE GDRasterBand<T>::getHomogeneousType() const
    {//{{{
      return homogType;
    }//}}}

    template<typename T>
    void GDRasterBand<T>::setProbabilityLag(T newProbaLag) 
    {//{{{
      probabilityLag = newProbaLag;
    }//}}}

    template<typename T>
    void GDRasterBand<T>::setSigmaFactor(T newSigmaFactor) 
    {//{{{
      sigmaFactor = newSigmaFactor;
    }//}}}

    template<typename T>
    bool GDRasterBand<T>::nearlySameValue(const RasterRect& rect) const
    {//{{{  
      int i,j;
      int Nmeans = 0;
      T tmp = 0;
      // Mean and value
      //T meanVariation = 0;
      T previousMeanValue = 0;
      T currentMeanValue = 0;
      T curValue = 0;
      T smallestValue = (T)1;
      T highestValue = 0;
      // The first value will be the mean and 
      T upBorderLimit   = (T)probabilityLag;
      T downBorderLimit = (T)probabilityLag;
      T zeroValue = (T)0;
      //--- Scanning the rectangle zone
      for (i=rect.xtopleft; i < rect.xtopleft+rect.xsize(); i++)
      {
        for (j=rect.ytopleft; j < rect.ytopleft+rect.ysize(); j++) 
        {
          //--- Getting value and average at i j
          curValue = get(RasterCellIndex(i,j));
          if (curValue < smallestValue)
          { smallestValue = curValue; }
          if (curValue > highestValue)
          { highestValue = curValue; }
          // Running mean of the raster rect zone
          if (Nmeans > 0)
          {
            currentMeanValue = (currentMeanValue * ((T)Nmeans/(T)(Nmeans+1))) + (curValue / ((T)(Nmeans+1)));
          } else { // Initial Case
            currentMeanValue  = curValue;
            previousMeanValue = currentMeanValue;
          }
          Nmeans++;

          //--- Checking if up and down limits are respected else zone is not homogeneous
          //meanVariation = previousMeanValue - currentMeanValue;
          // Zone not nearly homogeneous
          //if ((downBorderLimit - meanVariation) < zeroValue)
          //{ return false; }
          // Zone not nearly homogeneous
          //if ((upBorderLimit + meanVariation)   < zeroValue)
          //{ return false; }
          // Setting down-limit
          tmp = (currentMeanValue - (T)probabilityLag);
          if ( tmp < zeroValue)
          { tmp = zeroValue; }
          downBorderLimit = smallestValue - tmp;
          // smallest value out of zone segmentation
          if (downBorderLimit < zeroValue) { return false; }
          // Setting up-limit
          tmp = (currentMeanValue + (T)probabilityLag);
          if ( tmp > (T)1)
          { tmp = (T)1; }
          upBorderLimit = tmp - highestValue;
          // highest value out of zone segmentation
          if (upBorderLimit < zeroValue) { return false; }
          // adjusting for next step
          previousMeanValue = currentMeanValue;
        }
      }
      return true;
    }//}}}

    template<typename T>
    bool GDRasterBand<T>::sigmaSameValue(const RasterRect& rect) const
    {//{{{
      // All values are considered the same in the next case.
      if (sigmaFactor <= 0.0)
      {
        return false;
      }

      if (sigmaFactor > 3.0)
      {
        return true;
      }
      //-- else calculate and test sigma
      int i,j;
      int Nmeans = 0;
      T tmp = 0;
      // Mean and value
      T currentMeanValue =  0;
      T sigmaVal = 0;
      //--- Scanning the rectangle zone
      for (i=rect.xtopleft; i < rect.xtopleft+rect.xsize(); i++)
      {
        for (j=rect.ytopleft; j < rect.ytopleft+rect.ysize(); j++) 
        {
          //--- Getting value and finding average
          currentMeanValue = currentMeanValue + get(RasterCellIndex(i,j));
          Nmeans++;
        }
      }
      //--- Getting mean
      if (Nmeans > 0)
      {
        currentMeanValue = currentMeanValue / ((T)Nmeans);
      } else {
        return false;
      }

      //--- Calculating Sigma
      for (i=rect.xtopleft; i < rect.xtopleft+rect.xsize(); i++)
      {
        for (j=rect.ytopleft; j < rect.ytopleft+rect.ysize(); j++) 
        {
          //--- Getting value and finding average
          tmp = currentMeanValue - get(RasterCellIndex(i,j));
          sigmaVal = sigmaVal + tmp*tmp;
        }
      }
      sigmaVal = sqrt(sigmaVal/((T)Nmeans));

      //--- Checking for very low sigmas too narrow-probability-spread cells.
      if (sigmaVal <= (probabilityLag/100.0))
      {
        // is homogeneous
        return true;
      }
      //--- Checking for high sigmas too widely-probability-spread cells.
      if (sigmaVal > probabilityLag)
      {
        // is hereterogeneous
        return false;
      }

      //--- Checking for homogeneity
      for (i=rect.xtopleft; i < rect.xtopleft+rect.xsize(); i++)
      {
        for (j=rect.ytopleft; j < rect.ytopleft+rect.ysize(); j++) 
        {
          tmp = fabs(currentMeanValue - get(RasterCellIndex(i,j)));
          if (tmp > sigmaFactor * sigmaVal)
          {
            return false;
          }
        }
      }
      return true;
    }//}}}

    template<typename T>
    double GDRasterBand<T>::getDecompositionLimit() {
      return decompositionLimit;
    }

    template<typename T>
    void GDRasterBand<T>::setDecompositionLimit(double _dl) {
      decompositionLimit = _dl;
    }

   template<typename T>
   double GDRasterBand<T>::mean(const RasterRect& rect) const {
    double result = -1.0;

    // Bad case
      if (rect.size() <= 0)
    {
      return result;
    }

    // Get the data
    std::vector<T> data;
    get(rect,data);
    
    // Calculate the mean of data
    double sum = 0.0;
    for (int i=0; i<rect.size(); i++) {
      sum += (double)data[i];
    }
    result = sum/(double)rect.size();

    return result;
   }

   template<typename T>
   double  GDRasterBand<T>::sd(const RasterRect& rect) const {
    double result = -1.0;
    double mean   = GDRasterBand<T>::mean(rect);
    double tmp = 0.0;

    // Bad case
      if (rect.size() <= 0)
    {
      return result;
    }

    // get the data
    std::vector<T> data;
    get(rect,data);

    // Calculate the standard deviation
    result = 0.0;
    for (int i=0;i<rect.size();i++) {
      tmp = mean - (double)data[i];
      tmp = tmp * tmp;
      result += tmp;
    }
    result = result/(double)(rect.size());
  
    return result;
   }

   
//   template<typename T>
//   void* GDRasterBand<T>::getBlock(const RasterRect& rect) const {
//    
//    int i,j,k,m,l;
//    l=0;
//    int nXBlocks, nYBlocks, nXBlockSize, nYBlockSize;
//   
//    // the data size
//    int dataSize = GDALGetDataTypeSize(dataType);
// 
//    // alloc the value
//    void* value = malloc(dataSize*1);
//   
//    void* block = (unsigned char*)malloc(dataSize*rect.size());
//   
//    band->GetBlockSize( &nXBlockSize, &nYBlockSize );
//   
//    // if the wanted block has a small size
//    if (((nXBlockSize == 1) && (nYBlockSize==1)) || (rect.size()<= 4)) {
//      // read pixel by pixel
//      for (i=rect.xtopleft;i<rect.xtopleft+rect.xsize();i++) {
//       for (j=rect.ytopleft;j<rect.ytopleft+rect.ysize();j++) {
//        k = (i-rect.xtopleft)*rect.ysize() + (j-rect.ytopleft);
//        if (getValue(i,j,dataType,value)) {
//          swapPixelValue(block,value,k,0);
//          l++;
//        } 
//       }
//      }
//    } else {
//      // read data by block by block way
//      int        iXBlock, iYBlock;
//      int        iBlock_xbegin, iBlock_ybegin;
//      void* papyData;
//      RasterRect interRect;
// 
//      // read block by block
// 
//      // number of blocks of the raster band
//      nXBlocks = (band->GetXSize() + nXBlockSize - 1) / nXBlockSize;
//      nYBlocks = (band->GetYSize() + nYBlockSize - 1) / nYBlockSize;
//   
//      // buffer contains current read block from band data
//      papyData = CPLMalloc(dataSize * nXBlockSize * nYBlockSize);
//   
//      // iterate all blocks : iXBlock = horizontal block offset, iYBlock = vertical block offset
//      for( iXBlock = 0; iXBlock < nXBlocks; iXBlock++ ) {
//       for( iYBlock = 0; iYBlock < nYBlocks; iYBlock++ ) {
// 
//        // for each block, check if the wanted block data is inside the block
//        // compute the xroot,yroot of the current block
//        iBlock_xbegin = iXBlock*nXBlockSize;
//        iBlock_ybegin = iYBlock*nYBlockSize;
//      
//        // if the current read block contains data for wanted block
//        if (rect.intersects(RasterRect(iBlock_xbegin,iBlock_ybegin,nXBlockSize,nYBlockSize),interRect)) {
//          // cout << "read block " << iXBlock << ":" << iYBlock << "rect :" << interRect << endl; 
//          band->ReadBlock(iXBlock,iYBlock,papyData);
//          // for each pixel in the intersection
//          for (i=interRect.xtopleft;i<interRect.xtopleft+interRect.xsize();i++) {
//           for (j=interRect.ytopleft;j<interRect.ytopleft+interRect.ysize();j++) {
//            // from raster cell index (i,j), compute the index 
//            k = (i-rect.xtopleft)*rect.ysize() + (j-rect.ytopleft);
//            m = (i-iBlock_xbegin)*nYBlockSize + (j-iBlock_ybegin);
//            // cout << "read pixel " << i << ":" << j << " k=" << k<< " m=" << m << endl; 
//            swapPixelValue(block,papyData,k,m); 
//            // increment the read pixel number
//            l++;
//           }
//          }
//        }
//        if (l>=rect.size())
//          continue;     
//       }
//       if (l>=rect.size())
//        continue;
//      }
//   
//      free(papyData);  
//    }
// 
//    if (l!=rect.size()) 
//      cout << "[RasterBand] Read block error, should read " << rect.size() << " values instead of " << l << endl;
//   
//    return block;
//   }
// 
   

  }
}

#endif /* LGL_RASTERBAND_HPP */