src/terralib/raster/Band.cpp

Go to the documentation of this file.

/*  Copyright (C) 2008 National Institute For Space Research (INPE) - Brazil.

    This file is part of the TerraLib - a Framework for building GIS enabled applications.

    TerraLib 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 3 of the License,
    or (at your option) any later version.

    TerraLib 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 TerraLib. See COPYING. If not, write to
    TerraLib Team at <terralib-team@terralib.org>.
 */

/*!
  \file terralib/raster/Band.cpp

  \brief It gives access to values in one band (dimension) of a raster.
*/

// TerraLib

#include "Band.h"
#include "BandProperty.h"
#include "PositionIterator.h"
#include "Utils.h"

// STL
#include <cassert>
#include <limits>

te::rst::Band::Band(BandProperty* p, std::size_t idx)
  : m_property(p),
    m_idx(idx)
{
}

te::rst::Band::~Band()
{
  delete m_property;
}

te::rst::Band& te::rst::Band::operator=(const Band& rhs)
{
  if(this != &rhs)
  {
    delete m_property;

    m_property = rhs.m_property ? new BandProperty(*rhs.m_property) : nullptr;

    m_idx = rhs.m_idx;
  }

  return *this;
}

void te::rst::Band::getValue(unsigned int c, unsigned int r, std::complex<double>& value) const
{
  double vr = 0.0;

  double vi = 0.0;

  getValue(c, r, vr);

  getIValue(c, r, vi);

  value = std::complex<double>(vr, vi);
}

void te::rst::Band::setValue(unsigned int c, unsigned int r, const std::complex<double>& value)
{
  setValue(c, r, value.real());

  setIValue(c, r, value.imag());
}

std::complex<double> te::rst::Band::getMinValue(bool readall, unsigned int rs, unsigned int cs, unsigned int rf, unsigned int cf) const
{
  std::complex<double> pixel;
  double min_img = std::numeric_limits<double>::max();
  double min_real = std::numeric_limits<double>::max();
  double no_data = getProperty()->m_noDataValue;

  if ((rf == cf) && (rf == 0))
  {
    rf = getRaster()->getNumberOfRows() - 1;
    cf = getRaster()->getNumberOfColumns() - 1;

// read up to 1000 pixels in random positions
    const unsigned int maxInputPoints = 1000;
    if (readall == false && (rf * cf) > maxInputPoints)
    {
      std::vector<te::gm::Point*> randomPoints = te::rst::GetRandomPointsInRaster(*getRaster(), maxInputPoints);
      te::rst::PointSetIterator<double> pit = te::rst::PointSetIterator<double>::begin(getRaster(), randomPoints);
      te::rst::PointSetIterator<double> pitend = te::rst::PointSetIterator<double>::end(getRaster(), randomPoints);
      while (pit != pitend)
      {
        getValue(pit.getColumn(), pit.getRow(), pixel);
        ++pit;

        if(pixel.real() == no_data)
          continue;

        if (pixel.real() < min_real)
          min_real = pixel.real();

        if (pixel.imag() < min_img)
          min_img = pixel.imag();
      }
      
      while (!randomPoints.empty()) {
        delete randomPoints.back();
        randomPoints.pop_back();
      }

      return std::complex<double>(min_real, min_img);
    }
  }
  
// read all pixels in range
  for (unsigned r = rs; r <= rf; r++)
    for (unsigned c = cs; c <= cf; c++)
    {
      getValue(c, r, pixel);

      if(pixel.real() == no_data)
        continue;

      if (pixel.real() < min_real)
        min_real = pixel.real();

      if (pixel.imag() < min_img)
        min_img = pixel.imag();
    }

  return std::complex<double>(min_real, min_img);
}

std::complex<double> te::rst::Band::getMaxValue(bool readall, unsigned int rs, unsigned int cs, unsigned int rf, unsigned int cf) const
{
  std::complex<double> pixel;
  double max_img = std::numeric_limits<double>::min();
  double max_real = std::numeric_limits<double>::min();
  double no_data = getProperty()->m_noDataValue;

  if ((rf == cf) && (rf == 0))
  {
    rf = getRaster()->getNumberOfRows() - 1;
    cf = getRaster()->getNumberOfColumns() - 1;

// read up to 1000 pixels in random positions
    const unsigned int maxInputPoints = 1000;
    if (readall == false && (rf * cf) > maxInputPoints)
    {
      std::vector<te::gm::Point*> randomPoints = te::rst::GetRandomPointsInRaster(*getRaster(), maxInputPoints);
      te::rst::PointSetIterator<double> pit = te::rst::PointSetIterator<double>::begin(getRaster(), randomPoints);
      te::rst::PointSetIterator<double> pitend = te::rst::PointSetIterator<double>::end(getRaster(), randomPoints);
      while (pit != pitend)
      {
        getValue(pit.getColumn(), pit.getRow(), pixel);
        ++pit;

        if(pixel.real() == no_data)
          continue;

        if (pixel.real() > max_real)
          max_real = pixel.real();

        if (pixel.imag() > max_img)
          max_img = pixel.imag();
      }
      
      while (!randomPoints.empty()) {
        delete randomPoints.back();
        randomPoints.pop_back();
      }
      
      return std::complex<double>(max_real, max_img);
    }
  }
  
// read all pixels in range
  for (unsigned r = rs; r <= rf; r++)
    for (unsigned c = cs; c <= cf; c++)
    {
      getValue(c, r, pixel);

      if(pixel.real() == no_data)
        continue;

      if (pixel.real() > max_real)
        max_real = pixel.real();

      if (pixel.imag() > max_img)
        max_img = pixel.imag();
    }

  return std::complex<double>(max_real, max_img);
}

std::complex<double> te::rst::Band::getStdValue(unsigned int rs, unsigned int cs, unsigned int rf, unsigned int cf) const
{
  if ((rf == cf) && (rf  == 0))
  {
    rf = getRaster()->getNumberOfRows() - 1;

    cf = getRaster()->getNumberOfColumns() - 1;
  }

  std::complex<double> pixel;

  std::complex<double> mean = getMeanValue(rs, cs, rf, cf);

  std::complex<double> diff;

  std::complex<double> sumDiffs(0.0, 0.0);

  unsigned n = (rf-rs+1) * (cf-cs+1) - 1;

  if (n == 0)
    return std::complex<double> (1.0, 1.0);

  for (unsigned r = rs; r <= rf; r++)
    for (unsigned c = cs; c <= cf; c++)
    {
      getValue(c, r, pixel);

      diff = pixel - mean;

      sumDiffs += diff * diff;
    }

  return std::complex<double> (std::sqrt(sumDiffs.real() / n), std::sqrt(sumDiffs.imag() / n));
}

std::complex<double> te::rst::Band::getMeanValue(unsigned int rs, unsigned int cs, unsigned int rf, unsigned int cf) const
{
  if ((rf == cf) && (rf  == 0))
  {
    rf = getRaster()->getNumberOfRows() - 1;

    cf = getRaster()->getNumberOfColumns() - 1;
  }

  std::complex<double> pixel;

  std::complex<double> sumValue(0.0, 0.0);

  unsigned int n = (rf-rs+1) * (cf-cs+1);

  for (unsigned r = rs; r <= rf; r++)
    for (unsigned c = cs; c <= cf; c++)
    {
      getValue(c, r, pixel);

      sumValue += pixel;
    }

  if (n == 0)
    return std::complex<double> (0.0, 0.0);

  return std::complex<double> (sumValue.real() / n, sumValue.imag() / n);
}

std::map<double, unsigned> te::rst::Band::getHistogramR(unsigned int rs, unsigned int cs, unsigned int rf, unsigned int cf, unsigned int b) const
{
  std::map<double, unsigned> histogram;
  getHistogramR( rs, cs,rf, cf, b, 1, histogram );
  return histogram;
}

void te::rst::Band::getHistogramR( const unsigned int rowStart, const unsigned int colStart,
  const unsigned int finalRow, const unsigned int finalCol, const unsigned int histoBins,
  const unsigned int sampleStep, std::map<double, unsigned>& histogram ) const
{
  if( finalRow < rowStart )
  {
    throw te::rst::Exception( "Invalid Row start/final" );
  }
  if( finalCol < colStart )
  {
    throw te::rst::Exception( "Invalid Col start/final" );
  }
  if( sampleStep < 1 )
  {
    throw te::rst::Exception( "Invalid sample step" );
  }

  histogram.clear();

  // globals

  const bool processAllPixels = ( ( rowStart == 0 ) && ( colStart == 0 ) &&
    ( finalRow == 0 ) && ( finalCol == 0 ) );
  const long int nBlocksX = (long int)m_property->m_nblocksx;
  const long int nBlocksY = (long int)m_property->m_nblocksy;
  const long int blocksWidth = (long int)m_property->m_blkw;
  const long int blocksHeight = (long int)m_property->m_blkh;
  const long int rowsBound = processAllPixels ?  getRaster()->getNumberOfRows() :
    std::min( finalRow + 1, getRaster()->getNumberOfRows() );
  const long int colsBound = processAllPixels ? getRaster()->getNumberOfColumns() :
    std::min( finalCol + 1, getRaster()->getNumberOfColumns() );
  const double& noDataValue = m_property->m_noDataValue;
  long int blockXIndex = 0;
  long int blockYIndex = 0;
  long int blockColStart = 0;
  long int blockRowStart = 0;
  long int blockColsBound = 0;
  long int blockRowsBound = 0;
  double pixelValue = 0;
  long int samplesRowStart = 0;
  long int samplesColStart = 0;
  long int sampleRow = 0;
  long int sampleCol = 0;

  // find limits to divide into bins

  double pmin = std::numeric_limits<double>::max();
  double pmax = -1.0 * std::numeric_limits<double>::max();
  unsigned int pixelsCount = 0;

  for( blockYIndex = 0 ; blockYIndex < nBlocksY ; ++blockYIndex )
  {
    blockRowStart = blockYIndex * blocksHeight;
    blockRowsBound = std::min( rowsBound, ( blockYIndex + 1 ) * blocksHeight );

    samplesRowStart = blockRowStart % sampleStep;
    if( samplesRowStart )
    {
      samplesRowStart = blockRowStart + ( sampleStep - samplesRowStart );
    }
    else
    {
      samplesRowStart = blockRowStart;
    }


    for( blockXIndex = 0 ; blockXIndex < nBlocksX ; ++blockXIndex )
    {
      blockColStart = blockXIndex * blocksWidth;
      blockColsBound = std::min( colsBound, ( blockXIndex + 1 ) * blocksWidth );

      samplesColStart = blockColStart % sampleStep;
      if( samplesColStart )
      {
        samplesColStart = blockColStart + ( sampleStep - samplesColStart );
      }
      else
      {
        samplesColStart = blockColStart;
      }

      for( sampleRow = samplesRowStart ; sampleRow < blockRowsBound ; sampleRow += sampleStep )
      {
        for( sampleCol = samplesColStart ; sampleCol < blockColsBound; sampleCol += sampleStep )
        {
          getValue(sampleCol, sampleRow, pixelValue);

          if (pixelValue == noDataValue)
            continue;

          if (pixelValue > pmax)
            pmax = pixelValue;

          if (pixelValue < pmin)
            pmin = pixelValue;

          ++pixelsCount;
        }
      }
    }
  }

  if( pixelsCount == 0 )
  {
    return;
  }

  // create histogram with bins

  std::vector< unsigned int > histoVector;
  double delta = 0.0;

  if( histoBins == 0 )
  {
    if(
        ( m_property->m_type == te::dt::FLOAT_TYPE )
        ||
        ( m_property->m_type == te::dt::DOUBLE_TYPE )
      )
    {
      delta = ( (pmax - pmin) / ( (double)( 254 ) ) );
      histoVector.resize( 255, 0 );
    }
    else
    {
      delta = 1.0;
      histoVector.resize( (std::size_t)( (pmax - pmin) + 1 ), 0 );
    }
  }
  else
  {
    delta = ( (pmax - pmin) / ( (double)( histoBins - 1 ) ) );
    histoVector.resize( (std::size_t)histoBins, 0 );
  }

  if( delta == 0 )
  {
    if( pixelsCount )
    {
      histogram[ pmax ] = pixelsCount;
    }

    return;
  }

  // fill histogram

  std::size_t histoVectorIdx = 0;

  for( blockYIndex = 0 ; blockYIndex < nBlocksY ; ++blockYIndex )
  {
    blockRowStart = blockYIndex * blocksHeight;
    blockRowsBound = std::min( rowsBound, ( blockYIndex + 1 ) * blocksHeight );

    samplesRowStart = blockRowStart % sampleStep;
    if( samplesRowStart )
    {
      samplesRowStart = blockRowStart + ( sampleStep - samplesRowStart );
    }
    else
    {
      samplesRowStart = blockRowStart;
    }

    for( blockXIndex = 0 ; blockXIndex < nBlocksX ; ++blockXIndex )
    {
      blockColStart = blockXIndex * blocksWidth;
      blockColsBound = std::min( colsBound, ( blockXIndex + 1 ) * blocksWidth );

      samplesColStart = blockColStart % sampleStep;
      if( samplesColStart )
      {
        samplesColStart = blockColStart + ( sampleStep - samplesColStart );
      }
      else
      {
        samplesColStart = blockColStart;
      }

      for( sampleRow = samplesRowStart ; sampleRow < blockRowsBound ; sampleRow += sampleStep )
      {
        for( sampleCol = samplesColStart ; sampleCol < blockColsBound; sampleCol += sampleStep )
        {
          getValue(sampleCol, sampleRow, pixelValue);

          if (pixelValue != noDataValue)
          {
            histoVectorIdx = static_cast< size_t >( std::floor( ( pixelValue - pmin ) / delta ) );
            assert( histoVectorIdx < histoVector.size() );
            ++( histoVector[ histoVectorIdx ] );
          }
        }
      }
    }
  }

  const std::size_t histoVectorSize = histoVector.size();

  for ( histoVectorIdx = 0 ; histoVectorIdx < histoVectorSize ; ++histoVectorIdx )
  {
    histogram[ pmin + ( ((double)histoVectorIdx) * delta ) ] = histoVector[ histoVectorIdx ];
  }
}

std::map<double, unsigned> te::rst::Band::getHistogramI(unsigned int rs, unsigned int cs, unsigned int rf, unsigned int cf, unsigned int b) const
{
  if ((rf == cf) && (rf  == 0))
  {
    rf = getRaster()->getNumberOfRows() - 1;

    cf = getRaster()->getNumberOfColumns() - 1;
  }

  std::map<double, unsigned> hist;

  double pixel;

  if (b == 0)
    for (unsigned r = rs; r <= rf; r++)
      for (unsigned c = cs; c <= cf; c++)
      {
        getIValue(c, r, pixel);

        hist[pixel]++;
      }
  else
  {
// find limits to divide into bins
    double pmin = std::numeric_limits<double>::max();
    double pmax = -1.0 * std::numeric_limits<double>::max();

    for (unsigned r = rs; r <= rf; r++)
      for (unsigned c = cs; c <= cf; c++)
      {
        getIValue(c, r, pixel);

        if (pixel > pmax)
          pmax = pixel;

        if (pixel < pmin)
          pmin = pixel;
      }

// create histogram with bins
    double delta = (pmax * 1.000001 - pmin) / b;

    std::map<std::size_t, double> binsLocations;

    std::size_t location = 0;

    for (double bins = pmin; bins < pmax; bins += delta)
    {
      hist[bins] = 0;

      binsLocations[location++] = bins;
    }

// fill histogram
    for (unsigned r = rs; r <= rf; r++)
      for (unsigned c = cs; c <= cf; c++)
      {
        getIValue(c, r, pixel);

        location = (std::size_t) ((pixel - pmin) / delta);

        hist[binsLocations[location]]++;
      }

// removing empty bins from histogram
    for (double bins = pmin; bins < pmax; bins += delta)
      if(hist[bins] == 0)
        hist.erase(bins);

  }

  return hist;
}

std::complex<double> te::rst::Band::getScaleValue() const
{
  return m_property->m_valuesScale;
}

void te::rst::Band::setScaleValue(const std::complex<double> s)
{
  m_property->m_valuesScale = s;
}

std::complex<double> te::rst::Band::getOffsetValue() const
{
  return m_property->m_valuesOffset;
}

void te::rst::Band::setOffsetValue(const std::complex<double> o)
{
  m_property->m_valuesOffset = o;
}

te::rst::BandProperty* te::rst::Band::getProperty()
{
  return m_property;
}

const te::rst::BandProperty* te::rst::Band::getProperty() const
{
  return m_property;
}

te::rst::Band& te::rst::Band::callOperator(std::complex<double>(*f)(std::complex<double>, std::complex<double>), te::rst::Band& rhs)
{
  assert(getRaster()->getNumberOfRows() == rhs.getRaster()->getNumberOfRows());
  assert(getRaster()->getNumberOfColumns() == rhs.getRaster()->getNumberOfColumns());
  assert(getRaster()->getAccessPolicy() == te::common::RWAccess || getRaster()->getAccessPolicy() == te::common::WAccess);

  std::complex<double> lhsv, rhsv;

  std::vector<std::pair<unsigned, unsigned> > rcStPos, rcFPos;

  unsigned last_y;
  unsigned last_x;
  for (unsigned x = 0; x < (unsigned) getProperty()->m_nblocksx; x++)
  {
    for (unsigned y = 0; y < (unsigned) getProperty()->m_nblocksy; y++)
    {
      rcStPos.push_back(std::pair<unsigned, unsigned> (y * getProperty()->m_blkh, x * getProperty()->m_blkw));
      last_y = (y + 1) * getProperty()->m_blkh;
      last_x = (x + 1) * getProperty()->m_blkw;
      if (last_y > getRaster()->getNumberOfRows())
        last_y = getRaster()->getNumberOfRows();
      if (last_x > getRaster()->getNumberOfColumns())
        last_x = getRaster()->getNumberOfColumns();
      rcFPos.push_back(std::pair<unsigned, unsigned> (last_y, last_x));
    }
  }

// rasters without no data values
  if (getProperty()->m_noDataValue == std::numeric_limits<double>::max())
  {
    for (unsigned i = 0; i < rcStPos.size(); i++)
    {
      for (unsigned r = rcStPos[i].first; r < rcFPos[i].first; r++)
        for (unsigned c = rcStPos[i].second; c < rcFPos[i].second; c++)
        {
          getValue(c, r, lhsv);
          rhs.getValue(c, r, rhsv);

          lhsv = f(lhsv, rhsv);
          setValue(c, r, lhsv);
        }
      }
  }
// rasters with no data values
  else
  {
    for (unsigned i = 0; i < rcStPos.size(); i++)
    {
      for (unsigned r = rcStPos[i].first; r < rcFPos[i].first; r++)
        for (unsigned c = rcStPos[i].second; c < rcFPos[i].second; c++)
        {
// lhs data
          getValue(c, r, lhsv);
          if (lhsv.real() == getProperty()->m_noDataValue)
            continue;

// rhs data
          rhs.getValue(c, r, rhsv);
          if (rhsv.real() == rhs.getProperty()->m_noDataValue)
            continue;

          lhsv = f(lhsv, rhsv);
          setValue(c, r, lhsv);
        }
    }
  }

  return *this;
}

te::rst::Band& te::rst::Band::callOperator(std::complex<double>(*f)(std::complex<double>, std::complex<double>), std::complex<double>& cvalue)
{
  assert(getRaster()->getAccessPolicy() == te::common::RWAccess || getRaster()->getAccessPolicy() == te::common::WAccess);

  std::complex<double> lhsv;

  std::vector<std::pair<unsigned, unsigned> > rcStPos, rcFPos;

  unsigned last_y;
  unsigned last_x;
  for (unsigned x = 0; x < (unsigned) getProperty()->m_nblocksx; x++)
  {
    for (unsigned y = 0; y < (unsigned) getProperty()->m_nblocksy; y++)
    {
      rcStPos.push_back(std::pair<unsigned, unsigned> (y * getProperty()->m_blkh, x * getProperty()->m_blkw));
      last_y = (y + 1) * getProperty()->m_blkh;
      last_x = (x + 1) * getProperty()->m_blkw;
      if (last_y > getRaster()->getNumberOfRows())
        last_y = getRaster()->getNumberOfRows();
      if (last_x > getRaster()->getNumberOfColumns())
        last_x = getRaster()->getNumberOfColumns();
      rcFPos.push_back(std::pair<unsigned, unsigned> (last_y, last_x));
    }
  }

// rasters without no data values
  if (getProperty()->m_noDataValue == std::numeric_limits<double>::max())
  {
    for (unsigned i = 0; i < rcStPos.size(); i++)
    {
      for (unsigned r = rcStPos[i].first; r < rcFPos[i].first; r++)
        for (unsigned c = rcStPos[i].second; c < rcFPos[i].second; c++)
        {
          getValue(c, r, lhsv);

          lhsv = f(lhsv, cvalue);
          setValue(c, r, lhsv);
        }
      }
  }
// rasters with no data values
  else
  {
    for (unsigned i = 0; i < rcStPos.size(); i++)
    {
      for (unsigned r = rcStPos[i].first; r < rcFPos[i].first; r++)
        for (unsigned c = rcStPos[i].second; c < rcFPos[i].second; c++)
        {
// lhs data
          getValue(c, r, lhsv);
          if (lhsv.real() == getProperty()->m_noDataValue)
            continue;

          lhsv = f(lhsv, cvalue);
          setValue(c, r, lhsv);
        }
    }
  }

  return *this;
}

std::complex<double> plus(std::complex<double> lhs, std::complex<double> rhs)
{
  return lhs + rhs;
}

std::complex<double> minus(std::complex<double> lhs, std::complex<double> rhs)
{
  return lhs - rhs;
}

std::complex<double> times(std::complex<double> lhs, std::complex<double> rhs)
{
  return lhs * rhs;
}

std::complex<double> divide(std::complex<double> lhs, std::complex<double> rhs)
{
  return lhs / rhs;
}

te::rst::Band& te::rst::Band::operator+=(te::rst::Band& rhs)
{
  return callOperator(plus, rhs);
}

te::rst::Band& te::rst::Band::operator+=(std::complex<double>& cvalue)
{
  return callOperator(plus, cvalue);
}

te::rst::Band& te::rst::Band::operator-=(te::rst::Band& rhs)
{
  return callOperator(minus, rhs);
}

te::rst::Band& te::rst::Band::operator-=(std::complex<double>& cvalue)
{
  return callOperator(minus, cvalue);
}

te::rst::Band& te::rst::Band::operator*=(std::complex<double>& cvalue)
{
  return callOperator(times, cvalue);
}

te::rst::Band& te::rst::Band::operator*=(te::rst::Band& rhs)
{
  return callOperator(times, rhs);
}

te::rst::Band& te::rst::Band::operator/=(te::rst::Band& rhs)
{
  return callOperator(divide, rhs);
}

te::rst::Band& te::rst::Band::operator/=(std::complex<double>& cvalue)
{
  return callOperator(divide, cvalue);
}

int te::rst::Band::getBlockSize() const
{
  int blkw = m_property->m_blkw;

  int blkh = m_property->m_blkh;

  int pxlsize = GetPixelSize(m_property->getType());

  int blksize = blkw * blkh * pxlsize;

  return blksize;
}

te::rst::Band::Band(const Band& rhs)
  : m_property(nullptr),
    m_idx(rhs.m_idx)
{
  m_property = rhs.m_property ? new BandProperty(*rhs.m_property) : nullptr;
}