27 #ifndef __TERRALIB_RASTER_INTERNAL_POSITIONITERATOR_H
28 #define __TERRALIB_RASTER_INTERNAL_POSITIONITERATOR_H
31 #include "../common/STLUtils.h"
32 #include "../common/MathUtils.h"
33 #include "../core/logger/Logger.h"
34 #include "../geometry.h"
95 virtual std::complex< T >
operator()(
const unsigned int i)
const = 0;
98 virtual unsigned int getRow()
const = 0;
204 std::complex< T >
operator()(
const unsigned int i)
const;
206 unsigned int getRow()
const;
274 std::vector<std::pair<int, int> >& intersectionRanges )
const;
282 std::vector<std::pair<int, int> >& intersectionRanges )
const;
291 std::vector<std::pair<int, int> >& intersectionRanges )
const;
313 std::vector< te::gm::Geometry const * >& geoms)
const;
352 std::complex< T >
operator()(
const unsigned int i)
const;
354 unsigned int getRow()
const;
425 std::complex< T >
operator()(
const unsigned int i)
const;
427 unsigned int getRow()
const;
499 if( ! initialize( r, p, iterationType ) )
501 throw te::rst::Exception(
TE_TR(
"Iterator initialization error") );
509 if( ! initialize( r, p, ScanLineIterationT ) )
511 throw te::rst::Exception(
TE_TR(
"Iterator initialization error") );
528 std::vector<T> values( m_raster->getNumberOfBands());
531 for (
unsigned int b = 0; b < this->m_raster->getNumberOfBands(); b++)
533 m_raster->getValue( m_column, m_row, value, b);
534 values[b] = ((T) value);
542 this->m_raster->getValue(m_column, m_row, m_operatorBrackets_value, i);
544 return (T) m_operatorBrackets_value;
549 this->m_raster->getValue(m_column, m_row, m_operatorParenthesis_value, i);
551 return (std::complex< T >) m_operatorParenthesis_value;
566 if( m_column < m_currentRangeEndingColumn )
572 if ( ( m_currentLineIntersectionRangesIndex + 1 ) < (
int)( m_currentLineIntersectionRanges.size() ) )
574 ++m_currentLineIntersectionRangesIndex;
575 m_column = m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].first;
576 m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].second;
580 if( m_row < m_polygonEndingRow )
582 m_currentLineIntersectionRanges.clear();
584 while( m_row <= m_polygonEndingRow )
588 getIntersectionRanges( m_row, m_currentLineIntersectionRanges );
590 if( ! m_currentLineIntersectionRanges.empty() )
596 if( m_currentLineIntersectionRanges.empty() )
602 m_column = m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ 0 ].first;
603 m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ 0 ].second;
604 m_currentLineIntersectionRangesIndex = 0;
617 if(m_column > m_currentRangeStartingColumn )
623 if( m_currentLineIntersectionRangesIndex > 0 )
625 --m_currentLineIntersectionRangesIndex;
626 m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].first;
627 m_column = m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].second;
631 if( m_row > m_polygonStartingRow )
633 m_currentLineIntersectionRanges.clear();
635 while( m_row >= m_polygonStartingRow )
639 getIntersectionRanges( m_row, m_currentLineIntersectionRanges );
641 if( ! m_currentLineIntersectionRanges.empty() )
647 if( m_currentLineIntersectionRanges.empty() )
653 m_currentLineIntersectionRangesIndex =
static_cast<int>(m_currentLineIntersectionRanges.size()) - 1;
654 m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].first;
655 m_column = m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ m_currentLineIntersectionRangesIndex ].second;
748 return ( (this->m_row != rhs.
m_row) && (this->m_column != rhs.
m_column));
762 m_iterationType = ScanLineIterationT;
763 m_totalRasterColumns = 0;
764 m_totalRasterRows = 0;
765 m_polygonStartingColumn = -1;
766 m_polygonEndingColumn = -1;
767 m_polygonStartingRow = -1;
768 m_polygonEndingRow = -1;
772 m_currentRangeStartingColumn = -1;
773 m_currentRangeEndingColumn = -1;
774 m_currentLineIntersectionRangesIndex = -1;
779 m_tileIndexerPtr.reset();
780 m_currentLineIntersectionRanges.clear();
787 std::vector<std::pair<int, int> >& intersectionRanges )
const
789 switch( m_iterationType )
791 case ScanLineIterationT :
793 getScanLineIntersectionRanges( lineIndex, intersectionRanges );
796 case NDisjointBBOXIterationT :
798 getBBOXIntersectionRanges( lineIndex,
false, intersectionRanges );
801 case DisjointBBOXIterationT :
803 getBBOXIntersectionRanges( lineIndex,
true, intersectionRanges );
808 throw te::rst::Exception(
TE_TR(
"Invalid iteration type") );
816 std::vector<std::pair<int, int> >& intersectionRanges )
const
818 intersectionRanges.clear();
822 const double& rasterEnvelopeLLX =
823 m_raster->getGrid()->getExtent()->m_llx;
824 const double& rasterEnvelopeURX =
825 m_raster->getGrid()->getExtent()->m_urx;
829 const double lineY = this->m_raster->getGrid()->gridToGeo(0, (
double)lineIndex).y;
833 m_polygon->getSRID() );
835 currline.
setX(0, std::min(
837 m_polygon->getMBR()->getLowerLeft().x) );
838 currline.
setY(0, lineY);
840 currline.
setX(1, std::max(
842 m_polygon->getMBR()->getUpperRight().x ) );
843 currline.
setY(1, lineY);
846 std::vector<te::gm::Point*> intersectionPoints;
854 m_tileIndexerPtr->getTile(currline.
getY(0), ¤tTilePtr)
856 ( currentTilePtr != 0 )
860 std::unique_ptr<te::gm::Line> tileSeg;
861 std::unique_ptr< te::gm::Geometry > interResultPtr;
862 std::vector< te::gm::Geometry const* > singleGeomsPtrs;
863 std::size_t singleGeomsPtrsIdx = 0;
872 for (std::size_t i = 0; i < currentTilePtr->size(); i++)
874 assert((*currentTilePtr)[i].first < m_polygon->getNumRings());
877 assert(
dynamic_cast<te::gm::LinearRing const*
>((*m_polygon)[(*currentTilePtr)[i].first]));
880 assert((*currentTilePtr)[i].second < m_polygon->getNPoints());
884 ringTile->
getY((*currentTilePtr)[i].second),
887 ringTile->
getY((*currentTilePtr)[i].second + 1),
893 interResultPtr.reset( tileSeg->intersection( &currline ) );
895 if( interResultPtr.get() != 0 )
897 singleGeomsPtrs.clear();
898 breakIntoSingleGeoms( interResultPtr.get(), singleGeomsPtrs );
900 for( singleGeomsPtrsIdx = 0 ; singleGeomsPtrsIdx < singleGeomsPtrs.size() ;
901 ++singleGeomsPtrsIdx )
903 if( singleGeomsPtrs[ singleGeomsPtrsIdx ]->getGeomTypeId() ==
907 singleGeomsPtrs[ singleGeomsPtrsIdx ]->clone() ) );
909 else if( singleGeomsPtrs[ singleGeomsPtrsIdx ]->getGeomTypeId() ==
913 singleGeomsPtrsIdx ]);
915 if( lineStrPtr->
size() > 1 )
917 intersectionPoints.push_back( lineStrPtr->
getStartPoint().release() );
918 intersectionPoints.push_back( lineStrPtr->
getEndPoint().release() );
925 catch(
const std::exception& e)
927 TE_LOG_WARN(
"Geometry intersection error:" + e.what() );
931 for( std::size_t intersectionPointsIdx = 0; intersectionPointsIdx <
932 intersectionPoints.size() ; ++intersectionPointsIdx )
934 delete intersectionPoints[intersectionPointsIdx];
937 intersectionPoints.clear();
945 for( std::size_t intersectionPointsIdx = 0; intersectionPointsIdx <
946 intersectionPoints.size() ; ++intersectionPointsIdx )
948 delete intersectionPoints[intersectionPointsIdx];
951 intersectionPoints.clear();
958 std::size_t positionBegin = 0;
959 std::size_t positionEnd = 0;
960 std::vector<te::gm::Point*> intersectionPointsAux = intersectionPoints;
962 intersectionPoints.clear();
964 while( positionBegin < intersectionPointsAux.size() )
966 if( intersectionPointsAux[positionBegin] )
968 positionEnd = positionBegin + 1;
970 while( positionEnd < intersectionPointsAux.size() )
973 ( intersectionPointsAux[positionEnd] != 0 )
975 intersectionPointsAux[positionBegin]->equals(
976 intersectionPointsAux[positionEnd],
true)
979 delete intersectionPointsAux[positionEnd];
980 intersectionPointsAux[positionEnd] = 0;
986 intersectionPoints.push_back( intersectionPointsAux[positionBegin] );
994 std::sort(intersectionPoints.begin(), intersectionPoints.end(),
1000 std::size_t positionBegin = 0;
1001 std::size_t positionEnd = 0;
1002 int startingCol = 0;
1004 double startingX = 0;
1005 double startingY = 0;
1009 while( ( positionBegin + 1 ) < intersectionPoints.size() )
1011 positionEnd = positionBegin + 1;
1013 startingX = intersectionPoints[positionBegin]->getX();
1014 startingX = std::max( startingX, rasterEnvelopeLLX );
1016 startingY = intersectionPoints[positionBegin]->getY();
1018 endingX = intersectionPoints[positionEnd]->getX();
1019 endingX = std::min( endingX, rasterEnvelopeURX );
1021 endingY = intersectionPoints[positionEnd]->getY();
1028 ( startingX != endingX )
1032 (startingX + (endingX - startingX) / 2.0 ),
1033 (startingY +(endingY - startingY) / 2.0 ),
1034 m_polygon->getSRID()))
1041 startingCol = (int)std::ceil( m_raster->getGrid()->geoToGrid(startingX, startingY).x );
1042 startingCol = std::max( 0, startingCol );
1043 startingCol = std::min( m_totalRasterColumns - 1, startingCol );
1045 endingCol = (int)std::floor( m_raster->getGrid()->geoToGrid(endingX, endingY).x );
1046 endingCol = std::max( 0, endingCol );
1047 endingCol = std::min( m_totalRasterColumns - 1, endingCol );
1049 if( endingCol >= startingCol )
1051 intersectionRanges.push_back(std::pair<int, int>(startingCol, endingCol));
1055 positionBegin = positionEnd;
1062 while ( positionBegin < intersectionPoints.size() )
1064 delete intersectionPoints[positionBegin];
1069 template<
typename T>
1071 const bool areDisjoint, std::vector<std::pair<int, int> >& intersectionRanges )
const
1073 intersectionRanges.clear();
1075 double ulXCoord = 0;
1076 double ulYCoord = 0;
1077 m_raster->getGrid()->gridToGeo( ((
double)m_polygonStartingColumn) - 0.5,
1078 ((
double)lineIndex) - 0.5, ulXCoord, ulYCoord);
1080 const double yRes = this->m_raster->getGrid()->getResolutionY();
1081 const double xRes = this->m_raster->getGrid()->getResolutionX();
1084 m_polygon->getSRID(), 0 );
1088 bool rangeStarted =
false;
1089 std::pair<int, int> range;
1091 for(
int col = m_polygonStartingColumn ; col <= m_polygonEndingColumn ; ++col )
1094 currBBOXCoordsPtr[ 0 ].
x = currBBOXCoordsPtr[ 4 ].
x =
1095 ulXCoord + ( xRes * (double)( col - m_polygonStartingColumn ) );
1096 currBBOXCoordsPtr[ 0 ].
y = currBBOXCoordsPtr[ 4 ].
y = ulYCoord;
1099 currBBOXCoordsPtr[ 1 ].
x = currBBOXCoordsPtr[ 0 ].
x + xRes;
1100 currBBOXCoordsPtr[ 1 ].
y = ulYCoord;
1103 currBBOXCoordsPtr[ 2 ].
x = currBBOXCoordsPtr[ 1 ].
x;
1104 currBBOXCoordsPtr[ 2 ].
y = ulYCoord - yRes;
1107 currBBOXCoordsPtr[ 3 ].
x = currBBOXCoordsPtr[ 0 ].
x;
1108 currBBOXCoordsPtr[ 3 ].
y = ulYCoord - yRes;
1110 if( areDisjoint == currBBOX.
disjoint( m_polygon ) )
1118 rangeStarted =
true;
1119 range.first = range.second = col;
1126 rangeStarted =
false;
1127 intersectionRanges.push_back( range );
1135 intersectionRanges.push_back( range );
1139 template<
typename T>
1149 m_iterationType = iterationType;
1165 ( ll.
x > rasterEnvelope.
m_urx )
1167 ( ur.
x < rasterEnvelope.
m_llx )
1169 ( ur.
y < rasterEnvelope.
m_lly )
1171 ( ll.
y > rasterEnvelope.
m_ury )
1185 m_polygonStartingRow = (int)std::floor( urIndexed.
y );
1186 m_polygonStartingRow = std::max( 0, m_polygonStartingRow );
1187 m_polygonStartingRow = std::min( m_polygonStartingRow, ((
int)rasterGrid.
getNumberOfRows()) - 1 );
1189 m_polygonEndingRow = (int)std::ceil( llIndexed.
y );
1190 m_polygonEndingRow = std::max( 0, m_polygonEndingRow );
1191 m_polygonEndingRow = std::min( m_polygonEndingRow, ((
int)rasterGrid.
getNumberOfRows()) - 1 );
1193 m_polygonStartingColumn = (int)std::floor( llIndexed.
x );
1194 m_polygonStartingColumn = std::max( 0, m_polygonStartingColumn );
1195 m_polygonStartingColumn = std::min( m_polygonStartingColumn, ((
int)rasterGrid.
getNumberOfColumns()) - 1 );
1197 m_polygonEndingColumn = (int)std::ceil( urIndexed.
x );
1198 m_polygonEndingColumn = std::max( 0, m_polygonEndingColumn );
1199 m_polygonEndingColumn = std::min( m_polygonEndingColumn, ((
int)rasterGrid.
getNumberOfColumns()) - 1 );
1202 ( m_polygonEndingRow < m_polygonStartingRow )
1204 ( m_polygonEndingColumn < m_polygonStartingColumn )
1213 if( m_iterationType == ScanLineIterationT )
1220 m_row = m_polygonStartingRow;
1222 while( m_row <= m_polygonEndingRow )
1224 getIntersectionRanges( m_row, m_currentLineIntersectionRanges );
1226 if( ! m_currentLineIntersectionRanges.empty() )
1234 if( m_currentLineIntersectionRanges.empty() )
1240 m_column = m_currentRangeStartingColumn = m_currentLineIntersectionRanges[ 0 ].first;
1241 m_currentRangeEndingColumn = m_currentLineIntersectionRanges[ 0 ].second;
1242 m_currentLineIntersectionRangesIndex = 0;
1250 template<
typename T>
1253 std::vector< te::gm::Geometry const * >& geoms)
const
1263 for (std::size_t i = 0; i < nGeoms ; ++i)
1266 breakIntoSingleGeoms(currentGeom, geoms);
1271 geoms.push_back( g );
1282 m_currentpixelindex(0),
1290 m_currentpixelindex(0),
1295 throw te::rst::Exception(
TE_TR(
"Invalid line SRID") );
1315 double startingcolumn;
1317 std::unique_ptr<te::gm::Point> startpoint = inrasterline->
getStartPoint();
1319 startingcolumn, startingrow);
1321 double endingcolumn;
1323 std::unique_ptr<te::gm::Point> endpoint = inrasterline->
getEndPoint();
1325 endingcolumn, endingrow);
1330 double x1, x2, y1, y2, geoX, geoY;
1331 for(
int r = (
int)startingrow; r <= (int)endingrow; r++)
1332 for(
int c = (
int)startingcolumn; c <= (int)endingcolumn; c++)
1336 x1 = geoX - resXdiv2; y1 = geoY - resYdiv2;
1337 x2 = geoX + resXdiv2; y2 = geoY + resYdiv2;
1351 : m_raster(rhs.m_raster),
1352 m_currentpixelindex(rhs.m_currentpixelindex),
1353 m_pixelsinline(rhs.m_pixelsinline)
1359 m_pixelsinline.clear();
1364 std::vector<T> values(this->m_raster->getNumberOfBands());
1367 for (
unsigned int b = 0; b < this->m_raster->getNumberOfBands(); b++)
1369 this->m_raster->getValue(getColumn(), getRow(), value, b);
1370 values[b] = ((T) value);
1378 this->m_raster->getValue(getColumn(), getRow(), m_operatorBrackets_value, i);
1380 return (T) m_operatorBrackets_value;
1385 this->m_raster->getValue(getColumn(), getRow(), m_operatorParenthesis_value, i);
1387 return (std::complex< T >) m_operatorParenthesis_value;
1392 return (
unsigned int)(m_pixelsinline[m_currentpixelindex]->getY());
1397 return (
unsigned int)(m_pixelsinline[m_currentpixelindex]->getX());
1402 m_currentpixelindex++;
1404 if (m_currentpixelindex >= (
int)(m_pixelsinline.size()))
1410 m_currentpixelindex--;
1412 if (m_currentpixelindex < 0)
1430 template<
typename T>
1439 this->m_currentpixelindex = -1;
1461 template<
typename T>
1471 m_pixelsinpointset(0),
1472 m_currentpixelindex(0)
1478 m_pixelsinpointset(p),
1479 m_currentpixelindex(0)
1487 std::vector<te::gm::Point*> inside_points;
1494 throw te::rst::Exception(
TE_TR(
"Invalid point SRID") );
1513 : m_raster(rhs.m_raster),
1514 m_pixelsinpointset(rhs.m_pixelsinpointset),
1515 m_currentpixelindex(rhs.m_currentpixelindex)
1521 m_pixelsinpointset.clear();
1526 std::vector<T> values(this->m_raster->getNumberOfBands());
1529 for (
unsigned int b = 0; b < this->m_raster->getNumberOfBands(); b++)
1531 this->m_raster->getValue(getColumn(), getRow(), value, b);
1532 values[b] = ((T) value);
1540 this->m_raster->getValue(getColumn(), getRow(), m_operatorBrackets_value, i);
1542 return (T) m_operatorBrackets_value;
1547 this->m_raster->getValue(getColumn(), getRow(), m_operatorParenthesis_value, i);
1549 return (std::complex< T >) m_operatorParenthesis_value;
1554 return (
unsigned int)(m_pixelsinpointset[m_currentpixelindex]->getY());
1559 return (
unsigned int)(m_pixelsinpointset[m_currentpixelindex]->getX());
1564 m_currentpixelindex++;
1566 if (m_currentpixelindex >= (
int) m_pixelsinpointset.size())
1572 m_currentpixelindex--;
1574 if (m_currentpixelindex < 0)
1591 template<
typename T>
1600 this->m_currentpixelindex = -1;
1622 template<
typename T>
1632 #endif // __TERRALIB_RASTER_INTERNAL_POSITIONITERATOR_H