PositionIterator.h

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/rp/PositionIterator.h
 
  \brief It implements several ways to retrieve positions inside a band with spatial restriction,
         e.g. through a line, inside a bounding box or polygon, etc.
*/
 
#ifndef __TERRALIB_RASTER_INTERNAL_POSITIONITERATOR_H
#define __TERRALIB_RASTER_INTERNAL_POSITIONITERATOR_H
 
// Terralib
#include "../common/STLUtils.h"
#include "../common/MathUtils.h"
#include "../core/logger/Logger.h"
#include "../geometry.h"
#include "Band.h"
#include "BandProperty.h"
#include "BlockUtils.h"
#include "Grid.h"
#include "Exception.h"
 
// STL
#include <iostream>
#include <cmath>
 
namespace te
{
  namespace rst
  {
// Forward declaration.
    class Band;
 
    inline bool StdSortPointPointerComparison( te::gm::Point *p1, te::gm::Point *p2 )
    {
      return *p1 < *p2;
    }
 
    /*!
      \class AbstractPositionIterator
 
      \brief This class is the base for implementing ways to navigate over the band with spatial restriction,
             e.g. through a line, inside a bounding box or polygon, etc.
 
      \ingroup rst
    */
    template<typename T> class AbstractPositionIterator
    {
      public:
 
        /*! \brief Destructor. */
        virtual ~AbstractPositionIterator();
 
        /*! \brief Returns a vector of the values in current position (column, row) from iterator. */
        virtual const std::vector<T> operator*() const = 0;
 
        /*!
          \brief Returns the real value in current position (column, row, band) from iterator.
 
          \param i The band index.
 
          \return The pixel real value in current position.
 
          \note For complex values use operator()
        */
        virtual T operator[](const unsigned int i) const = 0;
 
        /*!
          \brief Returns the complex value in current position (column, row, band) from iterator.
 
          \param i The band index.
 
          \return The pixel comples value in current position.
 
          \note For real values use operator[]
        */
        virtual std::complex< T > operator()(const unsigned int i) const = 0;
 
        /*! \brief Returns the current row in iterator. */
        virtual unsigned int getRow() const = 0;
 
        /*! \brief Returns the current column in iterator. */
        virtual unsigned int getColumn() const = 0;
 
        /*! \brief Advances to the next position. */
        virtual void operator++() = 0;
 
        /*! \brief Returns to the previous position. */
        virtual void operator--() = 0;
 
        /*!
          \brief Assignment operator.
 
          \param rhs The right-hand-side copy used to copy from.
 
          \return A reference to this object.
        */
        virtual AbstractPositionIterator<T>& operator=(const AbstractPositionIterator<T>& rhs);
 
        /*!
          \brief Difference operator.
 
          \param rhs The right-hand side to compare.
 
          \return Returns true if the iterators are at different positions, or false otherwise.
        */
        virtual bool operator!=(const AbstractPositionIterator<T>& rhs) const = 0;
 
        /*!
          \brief Sets the iterator position to the end of the current band.
 
          \param b The band to retrieve the end information.
        */
        virtual void setEnd() = 0;
 
      protected:
 
        /*! \brief Constructor. */
        AbstractPositionIterator();
 
        /*!
          \brief Copy constructor.
 
          \param rhs The right-hand-side copy used to copy from.
        */
        AbstractPositionIterator(const AbstractPositionIterator& rhs );
 
    };
 
    /*!
      \class PolygonIterator
 
      \brief This class implements the strategy to iterate with spatial restriction,
             the iteration occurs inside a polygon.
 
      \ingroup rst
    */
    template<typename T> class PolygonIterator: public AbstractPositionIterator<T>
    {
      public:
 
        /*! \enum BlendMethod Pixel Blend methods. */
        enum IterationType
        {
          ScanLineIterationT, //!< Scan line iteration type.
          NDisjointBBOXIterationT, //!< Itaration over pixels whose bounding boxes are not disjoint in relation with the polygon.
          DisjointBBOXIterationT //!< Itaration over pixels whose bounding boxes are disjoint in relation with the polygon.
        };
 
        PolygonIterator();
 
        /*!
          \brief Constructor.
 
          \param b The band to iterate.
          \param p The polygon from where the iteration will navigate.
          \param iterationType Iteration type.
          \note Both raster and polygon must have the same SRID.
        */
        PolygonIterator(const te::rst::Raster* r, const te::gm::Polygon* p,
          const IterationType iterationType );
 
        /*!
          \brief Constructor.
 
          \param b The band to iterate.
          \param p The polygon from where the iteration will navigate.
          \note Both raster and polygon must have the same SRID.
          \note Iteration type: ScanLineIterationT.
        */
        PolygonIterator(const te::rst::Raster* r, const te::gm::Polygon* p);
 
        /*!
          \brief Copy constructor.
 
          \param rhs The right-hand-side copy used to copy from.
        */
        PolygonIterator(const PolygonIterator& rhs);
 
        ~PolygonIterator();
 
        const std::vector<T> operator*() const;
 
        T operator[](const unsigned int i) const;
 
        std::complex< T > operator()(const unsigned int i) const;
 
        unsigned int getRow() const;
 
        unsigned int getColumn() const;
 
        void operator++();
 
        void operator--();
 
        PolygonIterator& operator=(const PolygonIterator& rhs);
 
        AbstractPositionIterator<T>& operator=(const AbstractPositionIterator<T>& rhs);
 
        void setEnd();
 
        /*! \brief Returns an iterator referring to the first value of the band.*/
        static PolygonIterator begin(const te::rst::Raster* r, const te::gm::Polygon* p,
          const IterationType iterationType );
 
        /*! \brief Returns an iterator referring to the first value of the band.*/
        static PolygonIterator begin(const te::rst::Raster* r, const te::gm::Polygon* p);
 
        /*! \brief Returns an iterator referring to after the end of the iterator. */
        static PolygonIterator end(const te::rst::Raster* r, const te::gm::Polygon* p);
 
        /*! \brief Returns an iterator referring to after the end of the iterator. */
        static PolygonIterator end(const te::rst::Raster* r, const te::gm::Polygon* p,
          const IterationType iterationType );
 
        bool operator!=(const PolygonIterator<T>& rhs) const;
 
        bool operator!=(const AbstractPositionIterator<T>& rhs) const;
 
      protected:
 
        const te::rst::Raster* m_raster;  //!< The band from where to get the values (default 0).
        const te::gm::Polygon* m_polygon; //!< The spatial restriction to be applied in the iterator (default 0).
        IterationType m_iterationType;  //!< Current iteration type (default ScanLineIterationT).
        int m_totalRasterColumns; //!< The number of columns in band (default: 0).
        int m_totalRasterRows; //!< The number of rows in band (default: 0).
        int m_polygonStartingColumn; //!< Polygon first raster column (default: -1).
        int m_polygonEndingColumn; //!< Polygon last raster column (default: -1).
        int m_polygonStartingRow; //!< Polygon first raster row (default: -1).
        int m_polygonEndingRow;//!< Polygon last raster row (default: -1).
 
        int m_column; //!< The current column of the iterator (default -1).
        int m_row;  //!< The current row of the iterator (default -1).
        int m_currentRangeStartingColumn; //!< The current line iteration range starting column (default: -1).
        int m_currentRangeEndingColumn; //!< The current line iteration range column to finalize the iteration (default: -1).
        int m_currentLineIntersectionRangesIndex; //!< The actual line of the iterator (default -1).
        std::unique_ptr<te::rst::TileIndexer> m_tileIndexerPtr; //!< Tile indexer used to optimize the geometric operations
        std::vector<std::pair<int, int> > m_currentLineIntersectionRanges; //!< Coordinates of the columns to be transversed
 
        // Variables used by the operator[] method.
        mutable double m_operatorBrackets_value;
        mutable std::complex< double > m_operatorParenthesis_value;
 
        /*! \brief Initialize all internal variables to initial states. */
        void initiateVariables();
 
        /*! \brief Clear all internal allocated objects and reset back to the initial state. */
        void clear();
 
        /*!
            \brief Returns the given line intersection ranges (Coordinates of the columns to be transversed) using the current iteration type.
            \param lineIndex Line index.
            \param intersectionRanges Line intersection ranges (first:initial column, second:ending column) or an empty vector if no intersections exist.
        */
        void getIntersectionRanges( int lineIndex,
          std::vector<std::pair<int, int> >& intersectionRanges ) const;
 
        /*!
            \brief Returns the given line intersection ranges (Coordinates of the columns to be transversed) using the ScanLineIterationT mode.
            \param lineIndex Line index.
            \param intersectionRanges Line intersection ranges (first:initial column, second:ending column) or an empty vector if no intersections exist.
        */
        void getScanLineIntersectionRanges( int lineIndex,
          std::vector<std::pair<int, int> >& intersectionRanges ) const;
 
        /*!
            \brief Returns the given line intersection ranges using the pixel bounding box as reference.
            \param lineIndex Line index.
            \param areDisjoint If true, returns only thos pixels disjoint with the polygon, otherwise returns those pixels not disjoint with the polygon.
            \param intersectionRanges Line intersection ranges (first:initial column, second:ending column) or an empty vector if no intersections exist.
        */
        void getBBOXIntersectionRanges( int lineIndex, const bool areDisjoint,
          std::vector<std::pair<int, int> >& intersectionRanges ) const;
 
 
        /*!
          \brief Initialize this instance.
 
          \param r The raster to iterate.
          \param p The polygon from where the iteration will navigate.
          \param iterationType The iteration type
          \note Both raster and polygon must have the same SRID.
        */
        bool initialize(const te::rst::Raster* r, const te::gm::Polygon* p,
          const IterationType iterationType );
        
        /*!
        \brief It will get a GeometryCollection and distribute in a vector.
        \param g      Input GeometryCollection.
        \param geoms  Output Geometry Vector.
        \note If the geomSource is not a te::gm::GeometryCollectionType it will return vector with the input geometry.
        \note Do not delete de returned pointers since they point to parts of the original input geometry.
        */
        void breakIntoSingleGeoms( te::gm::Geometry const * const g, 
          std::vector< te::gm::Geometry const * >& geoms) const;        
 
    };
 
    /*!
      \class LineIterator
 
      \brief This class implements the strategy to iterate with spatial restriction,
             the iteration occurs inside a line.
 
      \ingroup rst
    */
    template<typename T> class LineIterator: public AbstractPositionIterator<T>
    {
      public:
 
        LineIterator();
 
        /*!
          \brief Constructor.
 
          \param b The band to iterate.
          \param l The line from where the iteration will navigate.
        */
        LineIterator(const te::rst::Raster* r, const te::gm::Line* l);
 
        /*!
          \brief Copy constructor.
 
          \param rhs The right-hand-side copy used to copy from.
        */
        LineIterator(const LineIterator& rhs);
 
        ~LineIterator();
 
        const std::vector<T> operator*() const;
 
        T operator[](const unsigned int i) const;
 
        std::complex< T > operator()(const unsigned int i) const;
 
        unsigned int getRow() const;
 
        unsigned int getColumn() const;
 
        void operator++();
 
        void operator--();
 
        LineIterator& operator=(const LineIterator& rhs);
 
        AbstractPositionIterator<T>& operator=(const AbstractPositionIterator<T>& rhs);
 
        void setEnd();
 
        /*! \brief Returns an iterator referring to the first value of the band.*/
        static LineIterator begin(const te::rst::Raster* r, const te::gm::Line* l);
 
        /*! \brief Returns an iterator referring to after the end of the iterator. */
        static LineIterator end(const te::rst::Raster* r, const te::gm::Line* l);
 
        bool operator!=(const LineIterator<T>& rhs) const;
 
        bool operator!=(const AbstractPositionIterator<T>& rhs) const;
 
      protected:
 
        const te::rst::Raster* m_raster;              //!< The band from where to get the values.
        const te::gm::Line* m_line;                   //!< The spatial restriction to be applied in the iterator.
        int m_currentpixelindex;                      //!< The index of the current pixel location.
        std::vector<te::gm::Point*> m_pixelsinline;   //!< A vector of pixel locations that intersects the line.
        mutable double m_operatorBrackets_value;                   //!< Used by the operator[] method.
        mutable std::complex< double > m_operatorParenthesis_value;     //!< Used by the operator() method.
 
 
    };
 
    /*!
      \class PointSetIterator
 
      \brief This class implements the strategy to iterate with spatial restriction,
             the iteration occurs inside a vector of points.
 
      \ingroup rst
    */
    template<typename T> class PointSetIterator: public AbstractPositionIterator<T>
    {
      public:
 
        PointSetIterator();
 
        /*!
          \brief Constructor.
 
          \param b The band to iterate.
          \param p The vector of points where the iteration will navigate.
        */
        PointSetIterator(const te::rst::Raster* r, const std::vector<te::gm::Point*> p);
 
        /*!
          \brief Copy constructor.
 
          \param rhs The right-hand-side copy used to copy from.
        */
        PointSetIterator(const PointSetIterator& rhs);
 
        ~PointSetIterator();
 
        const std::vector<T> operator*() const;
 
        T operator[](const unsigned int i) const;
 
        std::complex< T > operator()(const unsigned int i) const;
 
        unsigned int getRow() const;
 
        unsigned int getColumn() const;
 
        void operator++();
 
        void operator--();
 
        PointSetIterator& operator=(const PointSetIterator& rhs);
 
        AbstractPositionIterator<T>& operator=(const AbstractPositionIterator<T>& rhs);
 
        void setEnd();
 
        /*! \brief Returns an iterator referring to the first value of the band.*/
        static PointSetIterator begin(const te::rst::Raster* r, const std::vector<te::gm::Point*> p);
 
        /*! \brief Returns an iterator referring to after the end of the iterator. */
        static PointSetIterator end(const te::rst::Raster* r, const std::vector<te::gm::Point*> p);
 
        bool operator!=(const PointSetIterator<T>& rhs) const;
 
        bool operator!=(const AbstractPositionIterator<T>& rhs) const;
 
      protected:
 
        const te::rst::Raster* m_raster;              //!< The band from where to get the values.
        std::vector<te::gm::Point*> m_pixelsinpointset;    //!< The spatial restriction to be applied in the iterator.
        int m_currentpixelindex;                           //!< The index of the current pixel location.
        mutable double m_operatorBrackets_value;                   //!< Used by the operator[] method.
        mutable std::complex< double > m_operatorParenthesis_value;     //!< Used by the operator() method.
 
 
    };
 
    // ------------------------------------------------------------------------
 
// implementation of abstract position iterator
 
    template<typename T> te::rst::AbstractPositionIterator<T>::AbstractPositionIterator()
    {
    }
 
    template<typename T> te::rst::AbstractPositionIterator<T>::AbstractPositionIterator(const AbstractPositionIterator& )
    {
    }
 
    template<typename T> te::rst::AbstractPositionIterator<T>::~AbstractPositionIterator()
    {
    }
 
    template<typename T>
    AbstractPositionIterator<T>& te::rst::AbstractPositionIterator<T>::operator=(
      const AbstractPositionIterator<T>&)
    {
      return *this;
    }
 
    // ------------------------------------------------------------------------
 
// implementation of iteration strategy bounded by a polygon
    template<typename T> te::rst::PolygonIterator<T>::PolygonIterator()
    {
      initiateVariables();
    }
 
    template<typename T> te::rst::PolygonIterator<T>::PolygonIterator(
      const te::rst::Raster* r, const te::gm::Polygon* p,
      const IterationType iterationType)
    {
      initiateVariables();
 
      if( ! initialize( r, p, iterationType ) )
      {
        throw te::rst::Exception( TE_TR("Iterator initialization error") );
      }
    }
 
    template<typename T> te::rst::PolygonIterator<T>::PolygonIterator(const te::rst::Raster* r, const te::gm::Polygon* p)
    {
      initiateVariables();
 
      if( ! initialize( r, p, ScanLineIterationT ) )
      {
        throw te::rst::Exception( TE_TR("Iterator initialization error") );
      }
    }
 
    template<typename T> te::rst::PolygonIterator<T>::PolygonIterator(const PolygonIterator<T>& rhs)
    {
      initiateVariables();
      operator=( rhs );
    }
 
    template<typename T> te::rst::PolygonIterator<T>::~PolygonIterator()
    {
      clear();
    }
 
    template<typename T> const std::vector<T> te::rst::PolygonIterator<T>::operator*() const
    {
      std::vector<T> values( m_raster->getNumberOfBands());
      double value;
 
      for (unsigned int b = 0; b < this->m_raster->getNumberOfBands(); b++)
      {
        m_raster->getValue( m_column, m_row, value, b);
        values[b] = ((T) value);
      }
 
      return values;
    }
 
    template<typename T> T te::rst::PolygonIterator<T>::operator[](const unsigned int i) const
    {
      this->m_raster->getValue(m_column, m_row, m_operatorBrackets_value, i);
 
      return (T) m_operatorBrackets_value;
    }
 
    template<typename T> std::complex< T > te::rst::PolygonIterator<T>::operator()(const unsigned int i) const
    {
      this->m_raster->getValue(m_column, m_row, m_operatorParenthesis_value, i);
 
      return (std::complex< T >) m_operatorParenthesis_value;
    }
 
    template<typename T> unsigned int te::rst::PolygonIterator<T>::getRow() const
    {
      return m_row;
    }
 
    template<typename T> unsigned int te::rst::PolygonIterator<T>::getColumn() const
    {
      return m_column;
    }
 
    template<typename T> void te::rst::PolygonIterator<T>::operator++()
    {
      if( m_column < m_currentRangeEndingColumn )
      {
        ++m_column;
      }
      else
      {
        if ( ( m_currentLineIntersectionRangesIndex + 1 ) < (int)( m_currentLineIntersectionRanges.size() ) )
        {
          ++m_currentLineIntersectionRangesIndex;
          m_column = m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].first;
          m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].second;
        }
        else
        {
          if( m_row < m_polygonEndingRow )
          {
            m_currentLineIntersectionRanges.clear();
 
            while( m_row <= m_polygonEndingRow )
            {
              ++m_row;
 
              getIntersectionRanges( m_row, m_currentLineIntersectionRanges );
 
              if( ! m_currentLineIntersectionRanges.empty() )
              {
                break;
              }
            }
 
            if( m_currentLineIntersectionRanges.empty() )
            {
              setEnd();
            }
            else
            {
              m_column = m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ 0 ].first;
              m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ 0 ].second;
              m_currentLineIntersectionRangesIndex = 0;
            }
          }
          else
          {
            setEnd();
          }
        }
      }
    }
 
    template<typename T> void te::rst::PolygonIterator<T>::operator--()
    {
      if(m_column > m_currentRangeStartingColumn )
      {
        --m_column;
      }
      else
      {
        if( m_currentLineIntersectionRangesIndex > 0 )
        {
          --m_currentLineIntersectionRangesIndex;
          m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].first;
          m_column = m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].second;
        }
        else
        {
          if( m_row > m_polygonStartingRow )
          {
            m_currentLineIntersectionRanges.clear();
 
            while( m_row >= m_polygonStartingRow )
            {
              --m_row;
 
              getIntersectionRanges( m_row, m_currentLineIntersectionRanges );
 
              if( ! m_currentLineIntersectionRanges.empty() )
              {
                break;
              }
            }
 
            if( m_currentLineIntersectionRanges.empty() )
            {
              setEnd();
            }
            else
            {
              m_currentLineIntersectionRangesIndex = static_cast<int>(m_currentLineIntersectionRanges.size()) - 1;
              m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].first;
              m_column = m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].second;
            }
          }
          else
          {
            setEnd();
          }
        }
      }
    }
 
    template<typename T>
    te::rst::PolygonIterator<T>& te::rst::PolygonIterator<T>::operator=(
      const te::rst::PolygonIterator<T>& rhs)
    {
      clear();
 
      if (this != &rhs)
      {
        te::rst::AbstractPositionIterator<T>::operator=(rhs);
 
        m_raster = rhs.m_raster;
        m_polygon = rhs.m_polygon;
        m_iterationType = rhs.m_iterationType;
        m_totalRasterColumns = rhs.m_totalRasterColumns;
        m_totalRasterRows = rhs.m_totalRasterRows;
        m_polygonStartingColumn = rhs.m_polygonStartingColumn;
        m_polygonEndingColumn = rhs.m_polygonEndingColumn;
        m_polygonStartingRow = rhs.m_polygonStartingRow;
        m_polygonEndingRow = rhs.m_polygonEndingRow;
 
        m_column = rhs.m_column;
        m_row = rhs.m_row;
        m_currentRangeStartingColumn = rhs.m_currentRangeStartingColumn;
        m_currentRangeEndingColumn = rhs.m_currentRangeEndingColumn;
        m_currentLineIntersectionRangesIndex = rhs.m_currentLineIntersectionRangesIndex;
 
        if( rhs.m_tileIndexerPtr.get() )
        {
          m_tileIndexerPtr.reset(rhs.m_tileIndexerPtr->clone());
        }
 
        m_currentLineIntersectionRanges = rhs.m_currentLineIntersectionRanges;
      }
 
      return *this;
    }
 
    template<typename T>
    AbstractPositionIterator< T >& te::rst::PolygonIterator<T>::operator=(
      const AbstractPositionIterator< T >& rhs)
    {
      return operator=( *((te::rst::PolygonIterator<T>*)&rhs) );
    }
 
    template<typename T> void te::rst::PolygonIterator<T>::setEnd()
    {
      m_column = -1;
      m_row = -1;
    }
 
    template<typename T> te::rst::PolygonIterator<T> te::rst::PolygonIterator<T>::begin(
      const te::rst::Raster* r, const te::gm::Polygon* p, const IterationType iterationType )
    {
      return te::rst::PolygonIterator<T>(r, p, iterationType);
    }
 
    template<typename T> te::rst::PolygonIterator<T> te::rst::PolygonIterator<T>::begin(const te::rst::Raster* r, const te::gm::Polygon* p)
    {
      return te::rst::PolygonIterator<T>(r, p);
    }
 
    template<typename T> te::rst::PolygonIterator<T> te::rst::PolygonIterator<T>::end(
      const te::rst::Raster* r, const te::gm::Polygon* p, const IterationType iterationType)
    {
      te::rst::PolygonIterator<T> it(r, p, iterationType);
 
      it.setEnd();
 
      return it;
    }
 
    template<typename T> te::rst::PolygonIterator<T> te::rst::PolygonIterator<T>::end(const te::rst::Raster* r, const te::gm::Polygon* p)
    {
      te::rst::PolygonIterator<T> it(r, p);
 
      it.setEnd();
 
      return it;
    }
 
    template<typename T> bool te::rst::PolygonIterator<T>::operator!=(const te::rst::PolygonIterator<T>& rhs) const
    {
      return ( (this->m_row != rhs.m_row) && (this->m_column != rhs.m_column));
    }
 
    template<typename T>
    bool te::rst::PolygonIterator<T>::operator!=(
      const AbstractPositionIterator<T>& rhs) const
    {
      return operator!=( *((te::rst::PolygonIterator<T>*)&rhs) );
    }
 
    template<typename T> void te::rst::PolygonIterator<T>::initiateVariables()
    {
      m_raster = 0;
      m_polygon = 0;
      m_iterationType = ScanLineIterationT;
      m_totalRasterColumns = 0;
      m_totalRasterRows = 0;
      m_polygonStartingColumn = -1;
      m_polygonEndingColumn = -1;
      m_polygonStartingRow = -1;
      m_polygonEndingRow = -1;
 
      m_column = -1;
      m_row = -1;
      m_currentRangeStartingColumn = -1;
      m_currentRangeEndingColumn = -1;
      m_currentLineIntersectionRangesIndex = -1;
    }
 
    template<typename T> void te::rst::PolygonIterator<T>::clear()
    {
      m_tileIndexerPtr.reset();
      m_currentLineIntersectionRanges.clear();
 
      initiateVariables();
    }
 
    template<typename T>
    void te::rst::PolygonIterator<T>::getIntersectionRanges( int lineIndex,
       std::vector<std::pair<int, int> >& intersectionRanges ) const
    {
      switch( m_iterationType )
      {
        case ScanLineIterationT :
        {
          getScanLineIntersectionRanges( lineIndex, intersectionRanges );
          break;
        }
        case NDisjointBBOXIterationT :
        {
          getBBOXIntersectionRanges( lineIndex, false, intersectionRanges );
          break;
        }
        case DisjointBBOXIterationT :
        {
          getBBOXIntersectionRanges( lineIndex, true, intersectionRanges );
          break;
        }
        default :
        {
          throw te::rst::Exception( TE_TR("Invalid iteration type") );
          break;
        }
      }
    }
 
    template<typename T>
    void te::rst::PolygonIterator<T>::getScanLineIntersectionRanges( int lineIndex,
       std::vector<std::pair<int, int> >& intersectionRanges ) const
    {
      intersectionRanges.clear();
 
      // Globals
 
      const double& rasterEnvelopeLLX =
        m_raster->getGrid()->getExtent()->m_llx;
      const double& rasterEnvelopeURX =
        m_raster->getGrid()->getExtent()->m_urx;
 
      // current line
 
      const double lineY = this->m_raster->getGrid()->gridToGeo(0, (double)lineIndex).y;
 
      te::gm::Line currline( te::gm::Point( m_polygon->getSRID() ),
        te::gm::Point( m_polygon->getSRID() ), te::gm::LineStringType,
        m_polygon->getSRID() );
 
      currline.setX(0, std::min(
        rasterEnvelopeLLX,
        m_polygon->getMBR()->getLowerLeft().x) );
      currline.setY(0, lineY);
 
      currline.setX(1, std::max(
        rasterEnvelopeURX,
        m_polygon->getMBR()->getUpperRight().x ) );
      currline.setY(1, lineY);
 
      // Vector to store the points where the current line intersects the current tile of the polygon
      std::vector<te::gm::Point*> intersectionPoints;
 
      // Retrieve the tile corresponding to the current line and store the points where the line
      // intersects the tile
 
      te::rst::TileIndexer::TileSegIndex* currentTilePtr = 0;
 
      if(
          m_tileIndexerPtr->getTile(currline.getY(0), &currentTilePtr)
          &&
          ( currentTilePtr != 0 )
        )
      {
        te::gm::LinearRing const* ringTile;
        std::unique_ptr<te::gm::Line> tileSeg;
        std::unique_ptr< te::gm::Geometry > interResultPtr;
        std::vector< te::gm::Geometry const* > singleGeomsPtrs;
        std::size_t singleGeomsPtrsIdx = 0;
 
        try
        {
          // in some cases the intersection presents an unhandled exception, in this case we do not paint the current line
          // Transverse the segments of the tile retrieving the intersection between them and the current line
 
          te::gm::LineString* lineStrPtr = 0;
 
          for (std::size_t i = 0; i < currentTilePtr->size(); i++)
          {
            assert((*currentTilePtr)[i].first < m_polygon->getNumRings());
 
            // Retrieve the current ring of the current tile
            assert(dynamic_cast<te::gm::LinearRing const*>((*m_polygon)[(*currentTilePtr)[i].first]));
            ringTile = (te::gm::LinearRing const*)(*m_polygon)[(*currentTilePtr)[i].first];
 
            assert((*currentTilePtr)[i].second < m_polygon->getNPoints());
 
            // Retrieve the current segment of the current tile
            tileSeg.reset(new te::gm::Line(te::gm::Point(ringTile->getX((*currentTilePtr)[i].second),
                                                          ringTile->getY((*currentTilePtr)[i].second),
                                                          ringTile->getSRID()),
                                            te::gm::Point(ringTile->getX((*currentTilePtr)[i].second + 1),
                                                          ringTile->getY((*currentTilePtr)[i].second + 1),
                                                          ringTile->getSRID()),
                                            te::gm::LineStringType, ringTile->getSRID()));
 
            // Computes the intersection point between the current segment and the current line
 
            interResultPtr.reset( tileSeg->intersection( &currline ) );
 
            if( interResultPtr.get() != 0 )
            {
              singleGeomsPtrs.clear();
              breakIntoSingleGeoms( interResultPtr.get(), singleGeomsPtrs );
 
              for( singleGeomsPtrsIdx = 0 ; singleGeomsPtrsIdx < singleGeomsPtrs.size() ;
                ++singleGeomsPtrsIdx )
              {
                if( singleGeomsPtrs[ singleGeomsPtrsIdx ]->getGeomTypeId() ==
                  te::gm::PointType )
                {
                  intersectionPoints.push_back( (te::gm::Point*)(
                    singleGeomsPtrs[ singleGeomsPtrsIdx ]->clone() ) );
                }
                else if( singleGeomsPtrs[ singleGeomsPtrsIdx ]->getGeomTypeId() ==
                  te::gm::LineStringType )
                {
                  lineStrPtr = ((te::gm::LineString*) singleGeomsPtrs[
                    singleGeomsPtrsIdx ]);
                  
                  if( lineStrPtr->size() > 1 )
                  {
                    intersectionPoints.push_back( lineStrPtr->getStartPoint().release() );
                    intersectionPoints.push_back( lineStrPtr->getEndPoint().release() );
                  }
                }
              }
            }
          }
        }
        catch(const std::exception& e)
        {
          TE_LOG_WARN( "Geometry intersection error:" + e.what() );
 
          // deleting the intersections points
 
          for( std::size_t intersectionPointsIdx = 0; intersectionPointsIdx <
            intersectionPoints.size() ; ++intersectionPointsIdx )
          {
            delete intersectionPoints[intersectionPointsIdx];
          }
 
          intersectionPoints.clear();
        }
        catch(...)
        {
          TE_LOG_WARN( "Geometry intersection error" );
 
          // deleting the intersections points
 
          for( std::size_t intersectionPointsIdx = 0; intersectionPointsIdx <
            intersectionPoints.size() ; ++intersectionPointsIdx )
          {
            delete intersectionPoints[intersectionPointsIdx];
          }
 
          intersectionPoints.clear();
        }
      }
 
      // Removing duplicated points
 
      {
        std::size_t positionBegin = 0;
        std::size_t positionEnd = 0;
        std::vector<te::gm::Point*> intersectionPointsAux = intersectionPoints;
 
        intersectionPoints.clear();
 
        while( positionBegin < intersectionPointsAux.size() )
        {
          if( intersectionPointsAux[positionBegin] )
          {
            positionEnd = positionBegin + 1;
 
            while( positionEnd < intersectionPointsAux.size() )
            {
              if(
                  ( intersectionPointsAux[positionEnd] != 0 )
                  &&
                  intersectionPointsAux[positionBegin]->equals(
                    intersectionPointsAux[positionEnd], true)
                )
              {
                delete intersectionPointsAux[positionEnd];
                intersectionPointsAux[positionEnd] = 0;
              }
 
              positionEnd++;
            }
 
            intersectionPoints.push_back( intersectionPointsAux[positionBegin] );
          }
 
          ++positionBegin;
        }
      }
 
      // Sort the intersection points through its Y coordinates (column)
      std::sort(intersectionPoints.begin(), intersectionPoints.end(),
                te::rst::StdSortPointPointerComparison);
 
      // Using the intersection points, build a vector of coordinates (columns) with the start and
      // end column of the current line for each stretch
 
      std::size_t positionBegin = 0;
      std::size_t positionEnd = 0;
      int startingCol = 0;
      int endingCol = 0;
      double startingX = 0;
      double startingY = 0;
      double endingX = 0;
      double endingY = 0;
 
      while( ( positionBegin + 1 ) < intersectionPoints.size() )
      {
        positionEnd = positionBegin + 1;
 
        startingX = intersectionPoints[positionBegin]->getX();
        startingX = std::max( startingX, rasterEnvelopeLLX );
 
        startingY = intersectionPoints[positionBegin]->getY();
 
        endingX = intersectionPoints[positionEnd]->getX();
        endingX = std::min( endingX, rasterEnvelopeURX );
 
        endingY = intersectionPoints[positionEnd]->getY();
 
        // Build the vector of coordinates with the folowing structure:
        // vector<pair<startcolumn, endcolumn>>;
        // where each pair represents a stretch to be transversed
 
        if(
            ( startingX != endingX )
            &&
            (
                m_tileIndexerPtr->withinOrTouches(te::gm::Point(
                (startingX + (endingX - startingX) / 2.0 ),
                (startingY +(endingY - startingY) / 2.0 ),
                m_polygon->getSRID()))
            )
          )
        {
          // The middle-point (between the start end end point)
          // inside the current polygon
 
          startingCol = (int)std::ceil( m_raster->getGrid()->geoToGrid(startingX, startingY).x );
          startingCol = std::max( 0, startingCol );
          startingCol = std::min( m_totalRasterColumns - 1, startingCol );
 
          endingCol = (int)std::floor( m_raster->getGrid()->geoToGrid(endingX, endingY).x );
          endingCol = std::max( 0, endingCol );
          endingCol = std::min( m_totalRasterColumns - 1, endingCol );
 
          if( endingCol >= startingCol )
          {
            intersectionRanges.push_back(std::pair<int, int>(startingCol, endingCol));
          }
        }
 
        positionBegin = positionEnd;
      }
 
      // deleting the intersections points
 
      positionBegin = 0;
 
      while ( positionBegin < intersectionPoints.size() )
      {
        delete intersectionPoints[positionBegin];
        ++positionBegin;
      }
    }
 
    template<typename T>
    void te::rst::PolygonIterator<T>::getBBOXIntersectionRanges( int lineIndex,
       const bool areDisjoint, std::vector<std::pair<int, int> >& intersectionRanges ) const
    {
      intersectionRanges.clear();
 
      double ulXCoord = 0;
      double ulYCoord = 0;
      m_raster->getGrid()->gridToGeo( ((double)m_polygonStartingColumn) - 0.5,
        ((double)lineIndex) - 0.5, ulXCoord, ulYCoord);
 
      const double yRes = this->m_raster->getGrid()->getResolutionY();
      const double xRes = this->m_raster->getGrid()->getResolutionX();
 
      te::gm::LinearRing currBBOX( te::gm::LineStringType,
        m_polygon->getSRID(), 0 );
      currBBOX.setNumCoordinates( 5 );
      te::gm::Coord2D* currBBOXCoordsPtr = currBBOX.getCoordinates();
 
      bool rangeStarted = false;
      std::pair<int, int> range;
 
      for(  int col = m_polygonStartingColumn ; col <= m_polygonEndingColumn ; ++col )
      {
        // UL
        currBBOXCoordsPtr[ 0 ].x = currBBOXCoordsPtr[ 4 ].x =
          ulXCoord + ( xRes * (double)( col - m_polygonStartingColumn ) );
        currBBOXCoordsPtr[ 0 ].y = currBBOXCoordsPtr[ 4 ].y = ulYCoord;
 
        // UR
        currBBOXCoordsPtr[ 1 ].x = currBBOXCoordsPtr[ 0 ].x + xRes;
        currBBOXCoordsPtr[ 1 ].y = ulYCoord;
 
        // LR
        currBBOXCoordsPtr[ 2 ].x = currBBOXCoordsPtr[ 1 ].x;
        currBBOXCoordsPtr[ 2 ].y = ulYCoord - yRes;
 
        // LL
        currBBOXCoordsPtr[ 3 ].x = currBBOXCoordsPtr[ 0 ].x;
        currBBOXCoordsPtr[ 3 ].y = ulYCoord - yRes;
 
        if( areDisjoint == currBBOX.disjoint( m_polygon ) )
        { // this pixel must be considered
          if( rangeStarted )
          {
            range.second = col;
          }
          else
          {
            rangeStarted = true;
            range.first = range.second = col;
          }
        }
        else
        { // this pixel must NOT be considered
          if( rangeStarted )
          {
            rangeStarted = false;
            intersectionRanges.push_back( range );
          }
        }
      }
 
      // finalizing the last range
      if( rangeStarted )
      {
        intersectionRanges.push_back( range );
      }
    }
 
    template<typename T>
    bool te::rst::PolygonIterator<T>::initialize(const te::rst::Raster* r,
      const te::gm::Polygon* p, const IterationType iterationType )
    {
      clear();
 
      m_raster = r;
      m_polygon = p;
      m_totalRasterColumns = r->getNumberOfColumns();
      m_totalRasterRows = r->getNumberOfRows();
      m_iterationType = iterationType;
 
      if( r->getSRID() != p->getSRID() )
      {
        return false;
      }
 
      const te::rst::Grid& rasterGrid = *r->getGrid();
      const te::gm::Envelope& rasterEnvelope = *rasterGrid.getExtent();
 
      te::gm::Coord2D ll = m_polygon->getMBR()->getLowerLeft();
      te::gm::Coord2D ur = m_polygon->getMBR()->getUpperRight();
 
      // Is the polygon extent is outside the current raster area
 
      if(
          ( ll.x > rasterEnvelope.m_urx )
          ||
          ( ur.x < rasterEnvelope.m_llx )
          ||
          ( ur.y < rasterEnvelope.m_lly )
          ||
          ( ll.y > rasterEnvelope.m_ury )
        )
      {
        setEnd();
      }
      else
      {
        // defining starting/ending rows
 
        te::gm::Coord2D llIndexed;
        rasterGrid.geoToGrid(ll.x, ll.y, llIndexed.x, llIndexed.y);
        te::gm::Coord2D urIndexed;
        rasterGrid.geoToGrid(ur.x, ur.y, urIndexed.x, urIndexed.y);
 
        m_polygonStartingRow = (int)std::floor( urIndexed.y );
        m_polygonStartingRow = std::max( 0, m_polygonStartingRow );
        m_polygonStartingRow = std::min( m_polygonStartingRow, ((int)rasterGrid.getNumberOfRows()) - 1 );
 
        m_polygonEndingRow = (int)std::ceil( llIndexed.y );
        m_polygonEndingRow = std::max( 0, m_polygonEndingRow );
        m_polygonEndingRow = std::min( m_polygonEndingRow, ((int)rasterGrid.getNumberOfRows()) - 1 );
 
        m_polygonStartingColumn = (int)std::floor( llIndexed.x );
        m_polygonStartingColumn = std::max( 0, m_polygonStartingColumn );
        m_polygonStartingColumn = std::min( m_polygonStartingColumn, ((int)rasterGrid.getNumberOfColumns()) - 1 );
 
        m_polygonEndingColumn = (int)std::ceil( urIndexed.x );
        m_polygonEndingColumn = std::max( 0, m_polygonEndingColumn );
        m_polygonEndingColumn = std::min( m_polygonEndingColumn, ((int)rasterGrid.getNumberOfColumns()) - 1 );
 
        if(
            ( m_polygonEndingRow < m_polygonStartingRow )
            ||
            ( m_polygonEndingColumn < m_polygonStartingColumn )
          )
        {
          setEnd();
        }
        else
        {
          // initialize the TileIndexer
 
          if( m_iterationType == ScanLineIterationT )
          {
            m_tileIndexerPtr = std::unique_ptr<te::rst::TileIndexer>(new te::rst::TileIndexer(*m_polygon, rasterGrid.getResolutionY()));
          }
 
          // defining initial state
 
          m_row = m_polygonStartingRow;
 
          while( m_row <= m_polygonEndingRow )
          {
            getIntersectionRanges( m_row, m_currentLineIntersectionRanges );
 
            if( ! m_currentLineIntersectionRanges.empty() )
            {
              break;
            }
 
            ++m_row;
          }
 
          if( m_currentLineIntersectionRanges.empty() )
          {
            setEnd();
          }
          else
          {
            m_column = m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ 0 ].first;
            m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ 0 ].second;
            m_currentLineIntersectionRangesIndex = 0;
          }
        }
      }
 
      return true;
    }
    
    template<typename T>
    void te::rst::PolygonIterator<T>::breakIntoSingleGeoms(
      te::gm::Geometry const * const g, 
      std::vector< te::gm::Geometry const * >& geoms) const
    {
      te::gm::GeometryCollection const * const gc = 
        dynamic_cast< te::gm::GeometryCollection const * >(g);
        
      if (gc)
      {
        te::gm::Geometry const* currentGeom = 0;
        const std::size_t nGeoms = gc->getNumGeometries();
        
        for (std::size_t i = 0; i < nGeoms ; ++i)
        {
          currentGeom = gc->getGeometryN(i);
          breakIntoSingleGeoms(currentGeom, geoms);
        }
      }
      else
      {
        geoms.push_back( g );
      }
    }     
 
    // ------------------------------------------------------------------------
 
 
// implementation of iteration strategy bounded by a line
    template<typename T> te::rst::LineIterator<T>::LineIterator()
      : m_raster(0),
        m_line(0),
        m_currentpixelindex(0),
        m_pixelsinline(0)
    {
    }
 
    template<typename T> te::rst::LineIterator<T>::LineIterator(const te::rst::Raster* r, const te::gm::Line* l)
      : m_raster(r),
        m_line(l),
        m_currentpixelindex(0),
        m_pixelsinline(0)
    {
      if( r->getSRID() != l->getSRID() )
      {
        throw te::rst::Exception( TE_TR("Invalid line SRID") );
      }
 
      int srid = this->m_raster->getSRID();
 
// make intersection between line and band's envelope
      te::gm::Geometry* bandEnvelope = te::gm::GetGeomFromEnvelope(this->m_raster->getExtent(), srid);
      te::gm::Geometry* inter = bandEnvelope->intersection(m_line);
 
      if (inter->isEmpty())
      {
        setEnd();
 
        return;
      }
 
// create line that intersects only band's envelope
      te::gm::Line* inrasterline = (te::gm::Line*) inter;
 
// find starting and ending points
      double startingcolumn;
      double startingrow;
      std::unique_ptr<te::gm::Point> startpoint = inrasterline->getStartPoint();
      this->m_raster->getGrid()->geoToGrid(startpoint->getX(), startpoint->getY(),
                                           startingcolumn, startingrow);
 
      double endingcolumn;
      double endingrow;
      std::unique_ptr<te::gm::Point> endpoint = inrasterline->getEndPoint();
      this->m_raster->getGrid()->geoToGrid(endpoint->getX(), endpoint->getY(),
                                           endingcolumn, endingrow);
 
// creating one envelope per pixel, and intersects with line
      const double resXdiv2 = this->m_raster->getResolutionX() / 2;
      const double resYdiv2 = this->m_raster->getResolutionY() / 2;
      double x1, x2, y1, y2, geoX, geoY;
      for(int r = (int)startingrow; r <= (int)endingrow; r++)
        for(int c = (int)startingcolumn; c <= (int)endingcolumn; c++)
        {
// define envelope of pixel
          this->m_raster->getGrid()->gridToGeo(c, r, geoX, geoY);
          x1 = geoX - resXdiv2; y1 = geoY - resYdiv2;
          x2 = geoX + resXdiv2; y2 = geoY + resYdiv2;
 
          te::gm::Envelope* pixelbox = new te::gm::Envelope(x1, y1, x2, y2);
          te::gm::Geometry* pixelboxgeometry = GetGeomFromEnvelope(pixelbox, srid);
 
          if (te::gm::SatisfySpatialRelation(inrasterline, pixelboxgeometry, te::gm::INTERSECTS))
            m_pixelsinline.push_back(new te::gm::Point(c, r, srid));
        }
 
        if (m_pixelsinline.empty())
          setEnd();
    }
 
    template<typename T> te::rst::LineIterator<T>::LineIterator(const LineIterator<T>& rhs)
      : m_raster(rhs.m_raster),
        m_currentpixelindex(rhs.m_currentpixelindex),
        m_pixelsinline(rhs.m_pixelsinline)
    {
    }
 
    template<typename T> te::rst::LineIterator<T>::~LineIterator()
    {
      m_pixelsinline.clear();
    }
 
    template<typename T> const std::vector<T> te::rst::LineIterator<T>::operator*() const
    {
      std::vector<T> values(this->m_raster->getNumberOfBands());
      double value;
 
      for (unsigned int b = 0; b < this->m_raster->getNumberOfBands(); b++)
      {
        this->m_raster->getValue(getColumn(), getRow(), value, b);
        values[b] = ((T) value);
      }
 
      return values;
    }
 
    template<typename T> T te::rst::LineIterator<T>::operator[](const unsigned int i) const
    {
      this->m_raster->getValue(getColumn(), getRow(), m_operatorBrackets_value, i);
 
      return (T) m_operatorBrackets_value;
    }
 
    template<typename T> std::complex< T > te::rst::LineIterator<T>::operator()(const unsigned int i) const
    {
      this->m_raster->getValue(getColumn(), getRow(), m_operatorParenthesis_value, i);
 
      return (std::complex< T >) m_operatorParenthesis_value;
    }
 
    template<typename T> unsigned int te::rst::LineIterator<T>::getRow() const
    {
      return (unsigned int)(m_pixelsinline[m_currentpixelindex]->getY());
    }
 
    template<typename T> unsigned int te::rst::LineIterator<T>::getColumn() const
    {
      return (unsigned int)(m_pixelsinline[m_currentpixelindex]->getX());
    }
 
    template<typename T> void te::rst::LineIterator<T>::operator++()
    {
      m_currentpixelindex++;
 
      if (m_currentpixelindex >= (int)(m_pixelsinline.size()))
        setEnd();
    }
 
    template<typename T> void te::rst::LineIterator<T>::operator--()
    {
      m_currentpixelindex--;
 
      if (m_currentpixelindex < 0)
        setEnd();
    }
 
    template<typename T> te::rst::LineIterator<T>& te::rst::LineIterator<T>::operator=(const te::rst::LineIterator<T>& rhs)
    {
      if (this != &rhs)
      {
        te::rst::AbstractPositionIterator<T>::operator=(rhs);
 
        m_line = rhs.m_line;
        m_currentpixelindex = rhs.m_currentpixelindex;
        m_pixelsinline = rhs.m_pixelsinline;
      }
 
      return *this;
    }
 
    template<typename T>
    AbstractPositionIterator<T>& te::rst::LineIterator<T>::operator=(
      const AbstractPositionIterator<T>& rhs)
    {
      return operator=( *((LineIterator<T>*)&rhs) );
    }
 
    template<typename T> void te::rst::LineIterator<T>::setEnd()
    {
      this->m_currentpixelindex = -1;
    }
 
    template<typename T> te::rst::LineIterator<T> te::rst::LineIterator<T>::begin(const te::rst::Raster* r, const te::gm::Line* l)
    {
      return te::rst::LineIterator<T>(r, l);
    }
 
    template<typename T> te::rst::LineIterator<T> te::rst::LineIterator<T>::end(const te::rst::Raster* r, const te::gm::Line* l)
    {
      te::rst::LineIterator<T> it(r, l);
 
      it.setEnd();
 
      return it;
    }
 
    template<typename T> bool te::rst::LineIterator<T>::operator!=(const te::rst::LineIterator<T>& rhs) const
    {
      return ( (this->m_currentpixelindex != rhs.m_currentpixelindex) );
    }
 
    template<typename T>
    bool te::rst::LineIterator<T>::operator!=(
      const AbstractPositionIterator<T>& rhs) const
    {
      return operator!=( *((te::rst::LineIterator<T>*)&rhs) );
    }
 
// implementation of iteration strategy bounded by a vector of points
    template<typename T> te::rst::PointSetIterator<T>::PointSetIterator()
      : m_raster(0),
        m_pixelsinpointset(0),
        m_currentpixelindex(0)
    {
    }
 
    template<typename T> te::rst::PointSetIterator<T>::PointSetIterator(const te::rst::Raster* r, const std::vector<te::gm::Point*> p)
      : m_raster(r),
        m_pixelsinpointset(p),
        m_currentpixelindex(0)
    {
      const int rasterSRID = this->m_raster->getSRID();
 
      const te::gm::Envelope* rasterbox = r->getExtent();
      te::gm::Geometry* rasterboxgeometry = GetGeomFromEnvelope(rasterbox, rasterSRID);
 
// remove points that are not inside the band's envelope
      std::vector<te::gm::Point*> inside_points;
      double column;
      double row;
      for (unsigned int i = 0; i < m_pixelsinpointset.size(); i++)
      {
        if( rasterSRID != m_pixelsinpointset[i]->getSRID() )
        {
          throw te::rst::Exception( TE_TR("Invalid point SRID") );
        }
 
        if (te::gm::SatisfySpatialRelation(m_pixelsinpointset[i], rasterboxgeometry, te::gm::INTERSECTS))
        {
          this->m_raster->getGrid()->geoToGrid(m_pixelsinpointset[i]->getX(), m_pixelsinpointset[i]->getY(), column, row);
 
          inside_points.push_back(new te::gm::Point(column, row));
        }
      }
 
      m_pixelsinpointset.clear();
      m_pixelsinpointset = inside_points;
 
      if (m_pixelsinpointset.empty())
        setEnd();
    }
 
    template<typename T> te::rst::PointSetIterator<T>::PointSetIterator(const PointSetIterator<T>& rhs)
      : m_raster(rhs.m_raster),
        m_pixelsinpointset(rhs.m_pixelsinpointset),
        m_currentpixelindex(rhs.m_currentpixelindex)
    {
    }
 
    template<typename T> te::rst::PointSetIterator<T>::~PointSetIterator()
    {
      m_pixelsinpointset.clear();
    }
 
    template<typename T> const std::vector<T> te::rst::PointSetIterator<T>::operator*() const
    {
      std::vector<T> values(this->m_raster->getNumberOfBands());
      double value;
 
      for (unsigned int b = 0; b < this->m_raster->getNumberOfBands(); b++)
      {
        this->m_raster->getValue(getColumn(), getRow(), value, b);
        values[b] = ((T) value);
      }
 
      return values;
    }
 
    template<typename T> T te::rst::PointSetIterator<T>::operator[](const unsigned int i) const
    {
      this->m_raster->getValue(getColumn(), getRow(), m_operatorBrackets_value, i);
 
      return (T) m_operatorBrackets_value;
    }
 
    template<typename T> std::complex< T > te::rst::PointSetIterator<T>::operator()(const unsigned int i) const
    {
      this->m_raster->getValue(getColumn(), getRow(), m_operatorParenthesis_value, i);
 
      return (std::complex< T >) m_operatorParenthesis_value;
    }
 
    template<typename T> unsigned int te::rst::PointSetIterator<T>::getRow() const
    {
      return (unsigned int)(m_pixelsinpointset[m_currentpixelindex]->getY());
    }
 
    template<typename T> unsigned int te::rst::PointSetIterator<T>::getColumn() const
    {
      return (unsigned int)(m_pixelsinpointset[m_currentpixelindex]->getX());
    }
 
    template<typename T> void te::rst::PointSetIterator<T>::operator++()
    {
      m_currentpixelindex++;
 
      if (m_currentpixelindex >= (int) m_pixelsinpointset.size())
        setEnd();
    }
 
    template<typename T> void te::rst::PointSetIterator<T>::operator--()
    {
      m_currentpixelindex--;
 
      if (m_currentpixelindex < 0)
        setEnd();
    }
 
    template<typename T> te::rst::PointSetIterator<T>& te::rst::PointSetIterator<T>::operator=(const te::rst::PointSetIterator<T>& rhs)
    {
      if (this != &rhs)
      {
        te::rst::AbstractPositionIterator<T>::operator=(rhs);
 
        m_pixelsinpointset = rhs.m_pixelsinpointset;
        m_currentpixelindex = rhs.m_currentpixelindex;
      }
 
      return *this;
    }
 
    template<typename T>
    AbstractPositionIterator<T>& te::rst::PointSetIterator<T>::operator=(
      const AbstractPositionIterator<T>& rhs)
    {
      return operator=( *((te::rst::PointSetIterator<T>*)&rhs) );
    }
 
    template<typename T> void te::rst::PointSetIterator<T>::setEnd()
    {
      this->m_currentpixelindex = -1;
    }
 
    template<typename T> te::rst::PointSetIterator<T> te::rst::PointSetIterator<T>::begin(const te::rst::Raster* r, const std::vector<te::gm::Point*> p)
    {
      return te::rst::PointSetIterator<T>(r, p);
    }
 
    template<typename T> te::rst::PointSetIterator<T> te::rst::PointSetIterator<T>::end(const te::rst::Raster* r, const std::vector<te::gm::Point*> p)
    {
      te::rst::PointSetIterator<T> it(r, p);
 
      it.setEnd();
 
      return it;
    }
 
    template<typename T> bool te::rst::PointSetIterator<T>::operator!=(const te::rst::PointSetIterator<T>& rhs) const
    {
      return ( (this->m_currentpixelindex != rhs.m_currentpixelindex) );
    }
 
    template<typename T>
    bool te::rst::PointSetIterator<T>::operator!=(
      const AbstractPositionIterator<T>& rhs) const
    {
      return operator!=( *((te::rst::PointSetIterator<T>*)&rhs) );
    }
 
  } // end namespace rst
}   // end namespace te
 
#endif  // __TERRALIB_RASTER_INTERNAL_POSITIONITERATOR_H