ImgReaderGdal.h

/**********************************************************************
ImgReaderGdal.h: class to read raster files using GDAL API library
Copyright (C) 2008-2016 Pieter Kempeneers

This file is part of pktools

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

pktools 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 General Public License for more details.

You should have received a copy of the GNU General Public License
along with pktools.  If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/
#ifndef _IMGREADERGDAL_H_
#define _IMGREADERGDAL_H_

#include "ImgRasterGdal.h"
#include <assert.h>
#include <fstream>
#include <string>
#include <sstream>
#include "gdal_priv.h"
#include "base/Vector2d.h"

class ImgReaderGdal : public virtual ImgRasterGdal
{
public:
  ImgReaderGdal(void);
  ImgReaderGdal(const std::string& filename, const GDALAccess& readMode=GA_ReadOnly){open(filename, readMode);};
  ~ImgReaderGdal(void);
  void open(const std::string& filename, const GDALAccess& readMode=GA_ReadOnly);

  void close(void);
  template<typename T> void readData(T& value, int col, int row, int band=0);
  template<typename T> void readData(std::vector<T>& buffer, int minCol, int maxCol, int row, int band=0);
  template<typename T> void readData(std::vector<T>& buffer, int minCol, int maxCol, double row, int band=0, RESAMPLE resample=NEAR);
  template<typename T> void readDataBlock(Vector2d<T>& buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band=0);
  template<typename T> void readDataBlock(std::vector<T>& buffer , int minCol, int maxCol, int minRow, int maxRow, int band=0);
  template<typename T> void readData(std::vector<T>& buffer, int row, int band=0);
  template<typename T> void readData(std::vector<T>& buffer, double row, int band=0, RESAMPLE resample=NEAR);
  void getMinMax(int startCol, int endCol, int startRow, int endRow, int band, double& minValue, double& maxValue);
  void getMinMax(double& minValue, double& maxValue, int band=0);
  double getMin(int& col, int& row, int band=0);
  double getMax(int& col, int& row, int band=0);
  double getHistogram(std::vector<double>& histvector, double& min, double& max,unsigned int& nbin, int theBand=0, bool kde=false);
  void getRefPix(double& refX, double &refY, int band=0);
  void getRange(std::vector<short>& range, int Band=0);
  unsigned long int getNvalid(int band);

protected:
  void setCodec(const GDALAccess& readMode=GA_ReadOnly);
private:
};

template<typename T> void ImgReaderGdal::readData(T& value, int col, int row, int band)
{
  assert(band<nrOfBand()+1);
  assert(col<nrOfCol());
  assert(col>=0);
  assert(row<nrOfRow());
  assert(row>=0);
  double dvalue=0;
  double theScale=1;
  double theOffset=0;
  if(m_scale.size()>band||m_offset.size()>band){
    if(m_scale.size()>band)
      theScale=m_scale[band];
    if(m_offset.size()>band)
      theOffset=m_offset[band];
  }
  //fetch raster band
  GDALRasterBand  *poBand;
  poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
  poBand->RasterIO(GF_Read,col,row,1,1,&value,1,1,getGDALDataType<T>(),0,0);
  dvalue=theScale*value+theOffset;
  value=static_cast<T>(dvalue);
}

template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, int minCol, int maxCol, int row, int band)
{
  assert(band<nrOfBand()+1);
  assert(minCol<nrOfCol());
  assert(minCol>=0);
  assert(maxCol<nrOfCol());
  assert(minCol<=maxCol);
  assert(row<nrOfRow());
  assert(row>=0);
  double theScale=1;
  double theOffset=0;
  if(m_scale.size()>band||m_offset.size()>band){
    if(m_scale.size()>band)
      theScale=m_scale[band];
    if(m_offset.size()>band)
      theOffset=m_offset[band];
  }
  //fetch raster band
  GDALRasterBand  *poBand;
  poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
  if(buffer.size()!=maxCol-minCol+1)
    buffer.resize(maxCol-minCol+1);
  poBand->RasterIO(GF_Read,minCol,row,buffer.size(),1,&(buffer[0]),buffer.size(),1,getGDALDataType<T>(),0,0);
  if(m_scale.size()>band||m_offset.size()>band){
    for(int index=0;index<buffer.size();++index)
      buffer[index]=theScale*static_cast<double>(buffer[index])+theOffset;
  }
}

template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, int minCol, int maxCol, double row, int band, RESAMPLE resample)
{
  double eps=0.00001;
  std::vector<T> readBuffer_upper;
  std::vector<T> readBuffer_lower;
  double upperRow, lowerRow;
  if(buffer.size()!=maxCol-minCol+1)
    buffer.resize(maxCol-minCol+1);
  /* double upperRow=row-0.5; */
  /* upperRow=static_cast<int>(upperRow); */
  /* double lowerRow=row+0.5; */
  /* lowerRow=static_cast<int>(lowerRow); */
  switch(resample){
  case(BILINEAR):
    /* if(lowerRow>=nrOfRow()) */
    /*   lowerRow=nrOfRow()-1; */
    upperRow=static_cast<int>(row-0.5+eps);
    if(upperRow<0)
      upperRow=0;
    if(upperRow>=nrOfRow())
      upperRow=nrOfRow()-1;
    lowerRow=upperRow+1.0;
    if(lowerRow>=nrOfRow())
      lowerRow=nrOfRow()-1;

    readData(readBuffer_upper,minCol,maxCol,static_cast<int>(upperRow),band);
    readData(readBuffer_lower,minCol,maxCol,static_cast<int>(lowerRow),band);
    //do interpolation in y
    for(int icol=0;icol<maxCol-minCol+1;++icol){
      /* buffer[icol]=(lowerRow-row+0.5)*readBuffer_upper[icol]+(1-lowerRow+row-0.5)*readBuffer_lower[icol]; */
      if(!isNoData(readBuffer_upper[icol])){
        if(!isNoData(readBuffer_lower[icol])){
          buffer[icol]=(lowerRow-row+0.5)*readBuffer_upper[icol]+(1-lowerRow+row-0.5)*readBuffer_lower[icol];
        }
        else{
          buffer[icol]=readBuffer_upper[icol];
        }
      }
      else{
        buffer[icol]=readBuffer_lower[icol];
      }
    }
    break;
  default:
    readData(buffer,minCol,maxCol,static_cast<int>(row),band);
    break;
  }
}

template<typename T> void ImgReaderGdal::readDataBlock(Vector2d<T>& buffer2d, int minCol, int maxCol, int minRow, int maxRow, int band)
{
  buffer2d.resize(maxRow-minRow+1);
  typename std::vector<T> buffer;
  readDataBlock(buffer,minCol,maxCol,minRow,maxRow,band);
  typename std::vector<T>::const_iterator startit=buffer.begin();
  typename std::vector<T>::const_iterator endit=startit;
  for(int irow=minRow;irow<=maxRow;++irow){
    buffer2d[irow-minRow].resize(maxCol-minCol+1);
    endit+=maxCol-minCol+1;
    buffer2d[irow-minRow].assign(startit,endit);
    startit+=maxCol-minCol+1;
  }
}

template<typename T> void ImgReaderGdal::readDataBlock(std::vector<T>& buffer, int minCol, int maxCol, int minRow, int maxRow, int band)
{
  double theScale=1;
  double theOffset=0;
  if(m_scale.size()>band)
    theScale=m_scale[band];
  if(m_offset.size()>band)
    theOffset=m_offset[band];
  if(minCol>=nrOfCol() ||
     (minCol<0) ||
     (maxCol>=nrOfCol()) ||
     (minCol>maxCol) ||
     (minRow>=nrOfRow()) ||
     (minRow<0) ||
     (maxRow>=nrOfRow()) ||
     (minRow>maxRow)){
    std::string errorString="block not within image boundaries";
    throw(errorString);
  }
  if(buffer.size()!=(maxRow-minRow+1)*(maxCol-minCol+1))
    buffer.resize((maxRow-minRow+1)*(maxCol-minCol+1));
  //fetch raster band
  GDALRasterBand  *poBand;
  assert(band<nrOfBand()+1);
  poBand = m_gds->GetRasterBand(band+1);//GDAL uses 1 based index
  poBand->RasterIO(GF_Read,minCol,minRow,maxCol-minCol+1,maxRow-minRow+1,&(buffer[0]),(maxCol-minCol+1),(maxRow-minRow+1),getGDALDataType<T>(),0,0);
  if(m_scale.size()>band||m_offset.size()>band){
    for(int index=0;index<buffer.size();++index)
      buffer[index]=theScale*buffer[index]+theOffset;
  }
}

template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, int row, int band)
{
  readData(buffer,0,nrOfCol()-1,row,band);
}

template<typename T> void ImgReaderGdal::readData(std::vector<T>& buffer, double row, int band, RESAMPLE resample)
{
  readData(buffer,0,nrOfCol()-1,row,band,resample);
}

#endif // _IMGREADERGDAL_H_