27 #include "../common/progress/TaskProgress.h"
51 == 1,
"Invalid number of raster 1 bands" );
53 == 1,
"Invalid number of raster 2 bands" );
60 "Invalid m_moravecCorrelationWindowWidth" );
65 "Invalid m_moravecWindowWidth" );
69 "Invalid m_moravecMinAbsCorrelation" );
91 const double raster1ToRaster2TransfDMapError,
94 matchedInterestPoints.clear();
104 double raster1XRescFact = 1.0;
105 double raster1YRescFact = 1.0;
106 double raster2XRescFact = 1.0;
107 double raster2YRescFact = 1.0;
160 const double area1 = width1 * height1;
161 const double area2 = width2 * height2;
165 maxInterestPoints1 = (
unsigned int)( ((
double)maxInterestPoints1)*( area1 / area2 ) );
169 maxInterestPoints2 = (
unsigned int)( ((
double)maxInterestPoints2)*( area2 / area1 ) );
175 std::auto_ptr< te::common::TaskProgress > progressPtr;
179 progressPtr->setTotalSteps( 11 );
180 progressPtr->setMessage(
"Locating tie points" );
181 progressPtr->pulse();
182 if( ! progressPtr->isActive() )
return false;
191 std::vector< boost::shared_ptr< FloatsMatrix > > raster1Data;
193 double achievedRescaleFactorX = 0;
194 double achievedRescaleFactorY = 0;
211 achievedRescaleFactorX,
212 achievedRescaleFactorY ) )
217 raster1XRescFact = achievedRescaleFactorX;
218 raster1YRescFact = achievedRescaleFactorY;
224 progressPtr->pulse();
225 if( ! progressPtr->isActive() )
return false;
232 boost::shared_ptr< FloatsMatrix > tempMatrix(
237 raster1Data[ 0 ]->getMaxTmpFileSize(),
238 raster1Data[ 0 ]->getMaxMemPercentUsage() ),
239 "Cannot allocate image matrix" );
247 raster1Data[ 0 ]->reset();
248 raster1Data[ 0 ] = tempMatrix;
255 progressPtr->pulse();
256 if( ! progressPtr->isActive() )
return false;
265 raster1InterestPoints ) )
279 progressPtr->pulse();
280 if( ! progressPtr->isActive() )
return false;
289 raster1InterestPoints,
294 "Error generating raster 1 features" );
300 raster1InterestPoints.clear();
303 InterestPointsSetT::iterator itB = auxInterestPoints.begin();
304 const InterestPointsSetT::iterator itE = auxInterestPoints.end();
310 auxIP.
m_x = ( auxIP.
m_x / raster1XRescFact ) +
312 auxIP.
m_y = ( auxIP.
m_y / raster1YRescFact ) +
315 raster1InterestPoints.insert( auxIP );
323 progressPtr->pulse();
324 if( ! progressPtr->isActive() )
return false;
335 std::vector< boost::shared_ptr< FloatsMatrix > > raster2Data;
337 double achievedRescaleFactorX = 0;
338 double achievedRescaleFactorY = 0;
355 achievedRescaleFactorX,
356 achievedRescaleFactorY ) )
361 raster2XRescFact = achievedRescaleFactorX;
362 raster2YRescFact = achievedRescaleFactorY;
366 progressPtr->pulse();
367 if( ! progressPtr->isActive() )
return false;
374 boost::shared_ptr< FloatsMatrix > tempMatrix(
381 raster2Data[ 0 ]->getMaxTmpFileSize(),
382 raster2Data[ 0 ]->getMaxMemPercentUsage() ),
383 "Cannot allocate image matrix" );
391 raster2Data[ 0 ] = tempMatrix;
398 progressPtr->pulse();
399 if( ! progressPtr->isActive() )
return false;
408 raster2InterestPoints ) )
422 progressPtr->pulse();
423 if( ! progressPtr->isActive() )
return false;
432 raster2InterestPoints,
437 "Error generating raster 2 features" );
444 raster2InterestPoints.clear();
447 InterestPointsSetT::iterator itB = auxInterestPoints.begin();
448 const InterestPointsSetT::iterator itE = auxInterestPoints.end();
454 auxIP.
m_x = ( auxIP.
m_x / raster2XRescFact ) +
456 auxIP.
m_y = ( auxIP.
m_y / raster2YRescFact ) +
459 raster2InterestPoints.insert( auxIP );
467 progressPtr->pulse();
468 if( ! progressPtr->isActive() )
return false;
479 raster1InterestPoints,
480 raster2InterestPoints,
481 raster1ToRaster2TransfPtr,
482 raster1ToRaster2TransfDMapError,
483 internalMatchedInterestPoints ),
484 "Error matching features" );
488 progressPtr->pulse();
489 if( ! progressPtr->isActive() )
return false;
494 raster1Features.
clear();
495 raster2Features.
clear();
496 raster1InterestPoints.clear();
497 raster2InterestPoints.clear();
503 MatchedInterestPointsSetT::const_iterator itB = internalMatchedInterestPoints.begin();
504 const MatchedInterestPointsSetT::const_iterator itE = internalMatchedInterestPoints.end();
506 float minFeatureValue1 = FLT_MAX;
507 float maxFeatureValue1 = (-1.0) * FLT_MAX;
508 float minFeatureValue2 = FLT_MAX;
509 float maxFeatureValue2 = (-1.0) * FLT_MAX;
510 float minCorrelationValue = FLT_MAX;
511 float maxCorrelationValue = (-1.0) * FLT_MAX;
513 itB = internalMatchedInterestPoints.begin();
516 if( minFeatureValue1 > itB->m_point1.m_feature1 )
518 if( maxFeatureValue1 < itB->m_point1.m_feature1 )
519 maxFeatureValue1 = itB->m_point1.m_feature1;
521 if( minFeatureValue2 > itB->m_point2.m_feature1 )
522 minFeatureValue2 = itB->m_point2.m_feature1;
523 if( maxFeatureValue2 < itB->m_point2.m_feature1 )
524 maxFeatureValue2 = itB->m_point2.m_feature1;
526 if( minCorrelationValue > itB->m_feature )
527 minCorrelationValue = itB->m_feature;
528 if( maxCorrelationValue < itB->m_feature )
529 maxCorrelationValue = itB->m_feature;
534 float featureValue1Range = maxFeatureValue1 - minFeatureValue1;
535 float featureValue2Range = maxFeatureValue2 - minFeatureValue2;
536 float correlationValueRange = maxCorrelationValue - minCorrelationValue;
538 if( ( featureValue1Range == 0.0 ) || ( featureValue2Range == 0.0 ) ||
539 ( correlationValueRange == 0.0 ) )
541 itB = internalMatchedInterestPoints.begin();
545 auxMatchedPoints = *itB;
547 matchedInterestPoints.insert( auxMatchedPoints );
554 itB = internalMatchedInterestPoints.begin();
558 auxMatchedPoints = *itB;
562 ( auxMatchedPoints.
m_feature - minCorrelationValue )
564 correlationValueRange
582 matchedInterestPoints.insert( auxMatchedPoints );
591 progressPtr->pulse();
592 if( ! progressPtr->isActive() )
return false;
602 unsigned int returnValue = 0;
604 const unsigned int maxRastersArea = (
unsigned int)
627 returnValue = maxRastersArea /
636 const double freeVMem = 0.75 * std::min( totalPhysMem, ( totalVMem - usedVMem ) );
638 const double featureElementsNumber = (double)(
642 double maxFeaturesMemory =
647 (-2.0) * featureElementsNumber
650 ( 4.0 * featureElementsNumber * featureElementsNumber )
652 ( 4.0 * freeVMem / ((
double)
sizeof(
float ) ) )
661 (-2.0) * featureElementsNumber
664 ( 4.0 * featureElementsNumber * featureElementsNumber )
666 ( 4.0 * freeVMem / ((
double)
sizeof(
float ) ) )
671 returnValue = std::min( returnValue, (
unsigned int)maxFeaturesMemory );
677 FloatsMatrix& outputData,
const unsigned int iterationsNumber )
679 if( iterationsNumber == 0 )
681 outputData = inputData;
699 const unsigned int lastLineIndex = nLines - 1;
700 const unsigned int lastColIndex = nCols - 1;
701 unsigned int currLine = 0;
702 unsigned int currCol = 0;
708 if( iterationsNumber > 1 )
712 "Cannot allocate image matrix" );
717 for( currLine = 0 ; currLine < nLines ; ++currLine ) {
718 outputData( currLine, 0 ) = 0.0;
719 outputData( currLine, lastColIndex ) = 0.0;
722 for( currCol = 0 ; currCol < nCols ; ++currCol ) {
723 outputData( 0, currCol ) = 0.0;
724 outputData( lastLineIndex, currCol ) = 0.0;
732 float* outputLinePtr = 0;
733 unsigned int prevLine = 0;
734 unsigned int prevCol = 0;
735 unsigned int nextLine = 0;
736 unsigned int nextCol = 0;
738 for(
unsigned int iteration = 0 ; iteration < iterationsNumber ;
743 inputPtr = &inputData;
745 if( iterationsNumber > 1 )
746 outputPtr = &tempMatrix;
748 outputPtr = &outputData;
750 else if( iteration == iterationsNumber - 1 )
752 inputPtr = outputPtr;
753 outputPtr = &outputData;
758 inputPtr = outputPtr;
764 for( currLine = 1 ; currLine < lastLineIndex ; ++currLine )
766 prevLine = currLine - 1;
767 nextLine = currLine + 1;
769 outputLinePtr = outputPtr->operator[]( currLine );
771 for( currCol = 1 ; currCol < lastColIndex ; ++currCol )
773 prevCol = currCol - 1;
774 nextCol = currCol + 1;
776 outputLinePtr[ currCol ] =
778 internalInputMatrix( prevLine, prevCol )
779 + internalInputMatrix( prevLine, currCol )
780 + internalInputMatrix( prevLine, nextCol )
781 + internalInputMatrix( currLine, prevCol )
782 + internalInputMatrix( currLine, nextCol )
783 + internalInputMatrix( nextLine, prevCol )
784 + internalInputMatrix( nextLine, currCol )
785 + internalInputMatrix( nextLine, nextCol )
796 const unsigned int maxInterestPoints,
801 interestPoints.clear();
803 const unsigned int minRasterWidthAndHeight =
807 if( rasterData.
getLinesNumber() < minRasterWidthAndHeight )
return true;
809 bool returnValue =
true;
810 boost::mutex rastaDataAccessMutex;
811 boost::mutex interestPointsAccessMutex;
812 unsigned int nextRasterLinesBlockToProcess = 0;
813 std::vector< InterestPointsSetT > interestPointsSubSectors(
822 &nextRasterLinesBlockToProcess;
836 boost::thread_group threads;
841 threads.add_thread(
new boost::thread(
856 const unsigned int interestPointsSubSectorsSize = interestPointsSubSectors.size();
858 for(
unsigned int interestPointsSubSectorsIdx = 0 ; interestPointsSubSectorsIdx <
859 interestPointsSubSectorsSize ; ++interestPointsSubSectorsIdx )
861 interestPoints.insert(
862 interestPointsSubSectors[ interestPointsSubSectorsIdx ].begin(),
863 interestPointsSubSectors[ interestPointsSubSectorsIdx ].end() );
886 const unsigned int moravecWindowRadius = moravecWindowWidth / 2;
888 const unsigned int rowsBySubSector = (
unsigned int)std::ceil(
891 ((
double)tiePointsSubSectorsSplitFactor) );
892 const unsigned int colsBySubSector = (
unsigned int)std::ceil(
895 ((
double)tiePointsSubSectorsSplitFactor) );
898 const unsigned int lastBufferLineIdx = moravecWindowWidth - 1;
900 const unsigned int rasterBufferLineSizeBytes =
sizeof(
902 const unsigned int maskRasterBufferLineSizeBytes =
sizeof(
906 const unsigned int maxLinesPerProcessingBlock = (
unsigned int)
908 ((
double)rasterLines)
910 ((
double)processingBlocksNumber)
919 if( ! rasterBufferDataHandler.
reset( moravecWindowWidth, bufferCols,
928 boost::scoped_array< float* > rasterBufferHandler(
new float*[ moravecWindowWidth ] );
929 for(
unsigned int rasterBufferDataHandlerLine = 0 ; rasterBufferDataHandlerLine <
930 moravecWindowWidth ; ++rasterBufferDataHandlerLine )
932 rasterBufferHandler[ rasterBufferDataHandlerLine ] = rasterBufferDataHandler[
933 rasterBufferDataHandlerLine ];
936 float** rasterBufferPtr = rasterBufferHandler.get();
942 boost::scoped_array< unsigned char* > maskRasterBufferHandler(
new unsigned char*[ moravecWindowWidth ] );
944 unsigned char** maskRasterBufferPtr = 0;
948 if( ! maskRasterBufferDataHandler.
reset( moravecWindowWidth, bufferCols,
957 for(
unsigned int maskRasterBufferDataHandlerLine = 0 ; maskRasterBufferDataHandlerLine <
958 moravecWindowWidth ; ++maskRasterBufferDataHandlerLine )
960 maskRasterBufferHandler[ maskRasterBufferDataHandlerLine ] = maskRasterBufferDataHandler[
961 maskRasterBufferDataHandlerLine ];
964 maskRasterBufferPtr = maskRasterBufferHandler.get();
970 if( ! maximasBufferDataHandler.
reset( moravecWindowWidth, bufferCols,
979 boost::scoped_array< float* > maximasBufferHandler(
new float*[ moravecWindowWidth ] );
980 float** maximasBufferPtr = maximasBufferHandler.get();
981 unsigned int bufferCol = 0;
982 for(
unsigned int maximasBufferDataHandlerLine = 0 ; maximasBufferDataHandlerLine <
983 moravecWindowWidth ; ++maximasBufferDataHandlerLine )
985 maximasBufferHandler[ maximasBufferDataHandlerLine ] = maximasBufferDataHandler[
986 maximasBufferDataHandlerLine ];
987 for( bufferCol = 0 ; bufferCol < bufferCols ; ++bufferCol )
989 maximasBufferPtr[ maximasBufferDataHandlerLine ][ bufferCol ] = 0;
995 std::vector< InterestPointsSetT > interestPointsSubSectors(
999 for(
unsigned int rasterLinesBlockIdx = 0; rasterLinesBlockIdx <
1000 processingBlocksNumber ; ++rasterLinesBlockIdx )
1012 const unsigned int rasterLinesStart =
1018 (
int)(rasterLinesBlockIdx * maxLinesPerProcessingBlock )
1020 (
int)( 2 * moravecWindowRadius )
1023 const unsigned int rasterLinesEndBound =
1029 ( rasterLinesBlockIdx + 1 )
1031 maxLinesPerProcessingBlock
1034 ( 2 * moravecWindowRadius )
1037 const unsigned int varianceCalcStartRasterLineStart = rasterLinesStart +
1038 moravecWindowRadius;
1039 const unsigned int maximasLocationStartRasterLineStart = rasterLinesStart +
1040 4 * moravecWindowRadius;
1041 unsigned int windowStartBufCol = 0;
1042 const unsigned int windowEndBufColsBound = bufferCols -
1044 unsigned int windowStartBufOffset = 0;
1045 unsigned int windowStartBufXOffset = 0;
1046 unsigned int windowStartBufYOffset = 0;
1051 float diffValue = 0;
1052 bool isLocalMaxima =
false;
1054 float neighborMaximaDif = 0;
1055 unsigned int interestPointsSubSectorsIdx = 0;
1057 for(
unsigned int rasterLine = rasterLinesStart; rasterLine < rasterLinesEndBound ;
1065 memcpy( rasterBufferPtr[ lastBufferLineIdx ],
1067 rasterBufferLineSizeBytes );
1072 roolUpBuffer( maskRasterBufferPtr, moravecWindowWidth );
1073 memcpy( maskRasterBufferPtr[ lastBufferLineIdx ],
1075 maskRasterBufferLineSizeBytes );
1082 if( rasterLine >= varianceCalcStartRasterLineStart )
1086 for( windowStartBufCol = 0 ; windowStartBufCol < windowEndBufColsBound ;
1087 ++windowStartBufCol )
1089 const float& windowCenterPixelValue = rasterBufferPtr[
1090 moravecWindowRadius ][ windowStartBufCol +
1091 moravecWindowRadius ];
1097 for( windowStartBufOffset = 0 ; windowStartBufOffset <
1098 moravecWindowWidth ; ++windowStartBufOffset )
1100 diffValue = windowCenterPixelValue - rasterBufferPtr[
1101 moravecWindowRadius ][ windowStartBufCol +
1102 windowStartBufOffset ];
1103 horVar += ( diffValue * diffValue );
1105 diffValue = windowCenterPixelValue - rasterBufferPtr[
1106 windowStartBufOffset ][ windowStartBufCol +
1107 moravecWindowRadius ];
1108 verVar += ( diffValue * diffValue );
1110 diffValue = windowCenterPixelValue - rasterBufferPtr[
1111 windowStartBufOffset ][ windowStartBufCol +
1112 windowStartBufOffset ];
1113 diagVar += ( diffValue * diffValue );
1115 diffValue = windowCenterPixelValue - rasterBufferPtr[
1116 moravecWindowWidth - 1 - windowStartBufOffset ][ windowStartBufCol +
1117 windowStartBufOffset ];
1118 adiagVar += ( diffValue * diffValue );
1121 maximasBufferPtr[ lastBufferLineIdx ][ windowStartBufCol +
1122 moravecWindowRadius ] = std::min( horVar, std::min(
1123 verVar, std::min( diagVar, adiagVar ) ) );
1129 if( rasterLine >= maximasLocationStartRasterLineStart )
1131 for( windowStartBufCol = 0 ; windowStartBufCol < windowEndBufColsBound ;
1132 ++windowStartBufCol )
1134 isLocalMaxima =
true;
1135 const float& windowCenterPixelValue = maximasBufferPtr[
1136 moravecWindowRadius ][ windowStartBufCol +
1137 moravecWindowRadius ];
1140 for( windowStartBufYOffset = 0 ; windowStartBufYOffset <
1141 moravecWindowWidth ; ++windowStartBufYOffset )
1143 for( windowStartBufXOffset = 0 ; windowStartBufXOffset <
1144 moravecWindowWidth ; ++windowStartBufXOffset )
1146 neighborMaximaDif = windowCenterPixelValue - maximasBufferPtr[
1147 windowStartBufYOffset ][ windowStartBufCol +
1148 windowStartBufXOffset ];
1150 if( neighborMaximaDif < 0.0 )
1152 isLocalMaxima =
false;
1153 windowStartBufYOffset = moravecWindowWidth;
1157 auxInterestPoint.
m_feature1 += std::abs( neighborMaximaDif );
1163 auxInterestPoint.
m_x = windowStartBufCol +
1164 moravecWindowRadius;
1165 auxInterestPoint.
m_y = rasterLine - ( 2 * moravecWindowRadius );
1166 assert( auxInterestPoint.
m_x <
1168 assert( auxInterestPoint.
m_y <
1172 ( maskRasterBufferPtr == 0 )
1174 ( maskRasterBufferPtr[ 0 ][ auxInterestPoint.
m_x ] )
1177 interestPointsSubSectorsIdx = ( ( auxInterestPoint.
m_y /
1178 rowsBySubSector ) * tiePointsSubSectorsSplitFactor ) +
1179 ( auxInterestPoint.
m_x / colsBySubSector );
1180 assert( interestPointsSubSectorsIdx < interestPointsSubSectors.size() );
1182 interestPointsSubSectors[ interestPointsSubSectorsIdx ].insert(
1185 if( interestPointsSubSectors[ interestPointsSubSectorsIdx ].size() >
1186 maxInterestPointsBySubSector )
1188 interestPointsSubSectors[ interestPointsSubSectorsIdx ].erase(
1189 interestPointsSubSectors[ interestPointsSubSectorsIdx ].begin() );
1209 for(
unsigned int interestPointsSubSectorsIdx = 0 ; interestPointsSubSectorsIdx <
1210 interestPointsSubSectors.size() ; ++interestPointsSubSectorsIdx )
1213 interestPointsSubSectorsIdx ).
insert(
1214 interestPointsSubSectors[ interestPointsSubSectorsIdx ].begin(),
1215 interestPointsSubSectors[ interestPointsSubSectorsIdx ].end() );
1218 interestPointsSubSectorsIdx ).size() >
1219 maxInterestPointsBySubSector )
1222 interestPointsSubSectorsIdx ).erase(
1224 interestPointsSubSectorsIdx ).begin() );
1234 const unsigned int correlationWindowWidth,
1239 validInteresPoints.clear();
1244 const unsigned int rotated90CorralationWindowRadius = (
unsigned int)
1251 ( (
double)correlationWindowWidth )
1253 ( (
double)correlationWindowWidth )
1266 const unsigned int rasterDataLines = rasterData.
getLinesNumber();
1267 const unsigned int firstValidInterestPointX =
1268 rotated90CorralationWindowRadius + 1;
1269 const unsigned int lastValidInterestPointX = rasterDataCols
1270 - rotated90CorralationWindowRadius - 2;
1271 const unsigned int firstValidInterestPointY =
1272 rotated90CorralationWindowRadius + 1;
1273 const unsigned int lastValidInterestPointY = rasterDataLines
1274 - rotated90CorralationWindowRadius - 2;
1277 InterestPointsSetT::const_iterator iPointsIt = interestPoints.begin();
1278 const InterestPointsSetT::const_iterator iPointsItEnd = interestPoints.end();
1280 while( iPointsIt != iPointsItEnd )
1282 if( ( iPointsIt->m_x >= firstValidInterestPointX ) &&
1283 ( iPointsIt->m_x <= lastValidInterestPointX ) &&
1284 ( iPointsIt->m_y >= firstValidInterestPointY ) &&
1285 ( iPointsIt->m_y <= lastValidInterestPointY ) )
1287 validInteresPoints.insert( *iPointsIt );
1297 const unsigned int featureElemsNmb = correlationWindowWidth *
1298 correlationWindowWidth;
1299 const unsigned int featureSizeBytes =
sizeof( float ) *
1303 validInteresPoints.size(), featureElemsNmb ),
1304 "Cannot allocate features matrix" );
1308 boost::scoped_array< float > auxFeatureBufferHandler(
1309 new float[ featureElemsNmb ] );
1310 float* auxFeatureBufferPtr = auxFeatureBufferHandler.get();
1314 unsigned int curr_window_x_start = 0;
1315 unsigned int curr_window_y_start = 0;
1316 unsigned int curr_window_x_center = 0;
1317 unsigned int curr_window_y_center = 0;
1318 unsigned int curr_window_x_end = 0;
1319 unsigned int curr_window_y_end = 0;
1320 const unsigned int wind_radius = correlationWindowWidth / 2;
1321 const float wind_radius_double = (float)wind_radius;
1322 unsigned int curr_x = 0;
1323 unsigned int curr_y = 0;
1324 float curr_x_minus_radius = 0;
1325 float curr_y_minus_radius = 0;
1326 unsigned int curr_offset = 0;
1327 float int_x_dir = 0;
1328 float int_y_dir = 0;
1330 float rotated_curr_x = 0;
1331 float rotated_curr_y = 0;
1334 unsigned int rotated_curr_x_img = 0;
1335 unsigned int rotated_curr_y_img = 0;
1336 float featureElementsNormalizeFactor = 0.0;
1337 unsigned int featureElementIdx = 0;
1338 float* featurePtr = 0;
1339 float featureElementValue = 0.0;
1340 float featureElementMaxValue = 0.0;
1341 float featureElementMinValue = 0.0;
1343 InterestPointsSetT::const_iterator viPointsIt = validInteresPoints.begin();
1344 const InterestPointsSetT::const_iterator viPointsItEnd = validInteresPoints.end();
1345 unsigned int validInteresPointsIndex = 0 ;
1347 while( viPointsIt != viPointsItEnd )
1351 curr_window_x_center = viPointsIt->m_x;
1352 assert( curr_window_x_center >= rotated90CorralationWindowRadius );
1354 rotated90CorralationWindowRadius ) );
1355 curr_window_y_center = viPointsIt->m_y;
1356 assert( curr_window_y_center >= rotated90CorralationWindowRadius );
1357 assert( curr_window_y_center < ( rasterData.
getLinesNumber() - 1 -
1358 rotated90CorralationWindowRadius ) );
1360 curr_window_x_start = curr_window_x_center - wind_radius;
1361 curr_window_y_start = curr_window_y_center - wind_radius;
1362 curr_window_x_end = curr_window_x_start + correlationWindowWidth - 1;
1363 curr_window_y_end = curr_window_y_start + correlationWindowWidth - 1;
1369 for( curr_y = curr_window_y_start ; curr_y <= curr_window_y_end ;
1372 for( curr_offset = 0 ; curr_offset < wind_radius ;
1375 int_x_dir += rasterData( curr_y, curr_window_x_end - curr_offset ) -
1376 rasterData( curr_y, curr_window_x_start + curr_offset );
1384 for( curr_x = curr_window_x_start ; curr_x <= curr_window_x_end ;
1387 for( curr_offset = 0 ; curr_offset < wind_radius ; ++curr_offset )
1389 int_y_dir += rasterData( curr_window_y_start + curr_offset, curr_x ) -
1390 rasterData( curr_window_y_end - curr_offset, curr_x );
1396 int_norm = std::sqrt( ( int_x_dir * int_x_dir ) +
1397 ( int_y_dir * int_y_dir ) );
1399 if( int_norm != 0.0 ) {
1400 rot_cos = int_x_dir / int_norm;
1401 rot_sin = int_y_dir / int_norm;
1408 assert( ( rot_cos >= -1.0 ) && ( rot_cos <= 1.0 ) );
1409 assert( ( rot_sin >= -1.0 ) && ( rot_sin <= 1.0 ) );
1424 memset( auxFeatureBufferPtr, 0, featureSizeBytes );
1425 featureElementIdx = 0;
1426 featureElementMaxValue = -1.0 * FLT_MAX;
1427 featureElementMinValue = FLT_MAX;
1429 for( curr_y = 0 ; curr_y < correlationWindowWidth ; ++curr_y )
1431 for( curr_x = 0 ; curr_x < correlationWindowWidth ; ++curr_x )
1435 curr_x_minus_radius = ((float)curr_x) -
1437 curr_y_minus_radius = ((float)curr_y) -
1443 ( rot_cos * curr_x_minus_radius ) +
1444 ( rot_sin * curr_y_minus_radius );
1447 ( rot_cos * curr_y_minus_radius )
1448 - ( rot_sin * curr_x_minus_radius );
1453 rotated_curr_x += wind_radius_double;
1454 rotated_curr_y += wind_radius_double;
1458 rotated_curr_x_img = curr_window_x_start +
1459 (
unsigned int)
ROUND( rotated_curr_x );
1460 rotated_curr_y_img = curr_window_y_start +
1461 (
unsigned int)
ROUND( rotated_curr_y );
1463 featureElementValue = rasterData( rotated_curr_y_img,
1464 rotated_curr_x_img );
1466 auxFeatureBufferPtr[ featureElementIdx++ ] = featureElementValue;
1468 if( featureElementMaxValue < featureElementValue )
1469 featureElementMaxValue = featureElementValue;
1471 if( featureElementMinValue > featureElementValue )
1472 featureElementMinValue = featureElementValue;
1479 if( featureElementMaxValue == featureElementMinValue )
1480 featureElementsNormalizeFactor = 0.0;
1482 featureElementsNormalizeFactor = 1.0f / ( featureElementMaxValue -
1483 featureElementMinValue );
1485 featurePtr = features[ validInteresPointsIndex ];
1487 for( featureElementIdx = 0 ; featureElementIdx < featureElemsNmb ;
1488 ++featureElementIdx )
1490 featurePtr[ featureElementIdx ] = (
1491 ( auxFeatureBufferPtr[ featureElementIdx ] - featureElementMinValue )
1492 * featureElementsNormalizeFactor );
1493 assert( featurePtr[ featureElementIdx ] >= 0.0 );
1494 assert( featurePtr[ featureElementIdx ] <= 1.0 );
1497 ++validInteresPointsIndex;
1510 const double raster1ToRaster2TransfDMapError,
1513 matchedPoints.clear();
1515 const unsigned int interestPointsSet1Size = interestPointsSet1.size();
1516 if( interestPointsSet1Size == 0 )
return true;
1518 const unsigned int interestPointsSet2Size = interestPointsSet2.size();
1519 if( interestPointsSet2Size == 0 )
return true;
1522 assert( featuresSet1.
getLinesNumber() == interestPointsSet1Size );
1523 assert( featuresSet2.
getLinesNumber() == interestPointsSet2Size );
1527 InterestPointsSetT::const_iterator it1 = interestPointsSet1.begin();
1528 boost::scoped_array< InterestPointT > internalInterestPointsSet1(
1530 for(
unsigned int idx1 = 0 ; idx1 < interestPointsSet1Size ; ++idx1 )
1532 internalInterestPointsSet1[ idx1 ] = *it1;
1537 InterestPointsSetT::const_iterator it2 = interestPointsSet2.begin();
1538 boost::scoped_array< InterestPointT > internalInterestPointsSet2(
1540 for(
unsigned int idx2 = 0 ; idx2 < interestPointsSet2Size ; ++idx2 )
1542 internalInterestPointsSet2[ idx2 ] = *it2;
1552 "Error crearting the correlation matrix" );
1554 unsigned int col = 0;
1555 unsigned int line = 0;
1558 for( line = 0 ; line < interestPointsSet1Size ; ++line )
1560 linePtr = corrMatrix[ line ];
1562 for( col = 0 ; col < interestPointsSet2Size ; ++col )
1568 boost::mutex syncMutex;
1569 unsigned int nextFeatureIdx1ToProcess = 0;
1591 boost::thread_group threads;
1593 for(
unsigned int threadIdx = 0 ; threadIdx < procsNumber ;
1596 threads.add_thread(
new boost::thread(
1610 std::vector< float > eachLineMaxABSValues( interestPointsSet1Size,
1612 std::vector< unsigned int > eachLineMaxABSIndexes( interestPointsSet1Size,
1613 interestPointsSet2Size );
1614 std::vector< float > eachColMaxABSValues( interestPointsSet2Size,
1616 std::vector< unsigned int > eachColMaxABSIndexes( interestPointsSet2Size,
1617 interestPointsSet1Size );
1620 for( line = 0 ; line < interestPointsSet1Size ; ++line )
1622 linePtr = corrMatrix[ line ];
1624 for( col = 0 ; col < interestPointsSet2Size ; ++col )
1626 absValue = std::abs( linePtr[ col ] );
1630 if( absValue > eachLineMaxABSValues[ line ] )
1632 eachLineMaxABSValues[ line ] = absValue;
1633 eachLineMaxABSIndexes[ line ] = col;
1636 if( absValue > eachColMaxABSValues[ col ] )
1638 eachColMaxABSValues[ col ] = absValue;
1639 eachColMaxABSIndexes[ col ] = line;
1649 for( line = 0 ; line < interestPointsSet1Size ; ++line )
1651 col = eachLineMaxABSIndexes[ line ];
1653 if( ( col < interestPointsSet2Size ) &&
1654 ( eachColMaxABSIndexes[ col ] == line ) )
1656 auxMatchedPoints.
m_point1 = internalInterestPointsSet1[ line ];
1657 auxMatchedPoints.
m_point2 = internalInterestPointsSet2[ col ];
1658 auxMatchedPoints.
m_feature = std::abs( corrMatrix( line, col ) );
1660 matchedPoints.insert( auxMatchedPoints );
1681 unsigned int feat2Idx = 0;
1682 float const* feat1Ptr = 0;
1683 float const* feat2Ptr = 0;
1684 float* corrMatrixLinePtr = 0;
1685 unsigned int featCol = 0;
1689 float ccorrelation = 0;
1694 std::auto_ptr< te::gm::GeometricTransformation > raster1ToRaster2TransfPtr;
1704 const unsigned int featuresSet1Size =
1706 const unsigned int featuresSet2Size =
1711 std::vector< unsigned int > selectedFeaturesSet2Indexes;
1712 unsigned int selectedFeaturesSet2IndexesSize = 0;
1716 for(
unsigned int feat2Idx = 0 ; feat2Idx < featuresSet2Size ; ++feat2Idx )
1718 interestPointsSet2RTree.
insert(
1729 selectedFeaturesSet2Indexes.resize( featuresSet2Size );
1730 selectedFeaturesSet2IndexesSize = featuresSet2Size;
1731 for(
unsigned int feat2Idx = 0 ; feat2Idx < featuresSet2Size ; ++feat2Idx )
1733 selectedFeaturesSet2Indexes[ feat2Idx ] = feat2Idx;
1739 for(
unsigned int feat1Idx = 0 ; feat1Idx < featuresSet1Size ; ++feat1Idx )
1751 raster1ToRaster2TransfPtr->directMap(
1755 auxEnvelope.
m_lly );
1760 auxEnvelope.
m_llx -= interestPointsSet2RTreeSearchRadius;
1761 auxEnvelope.
m_lly -= interestPointsSet2RTreeSearchRadius;
1762 auxEnvelope.
m_urx += interestPointsSet2RTreeSearchRadius;
1763 auxEnvelope.
m_ury += interestPointsSet2RTreeSearchRadius;
1765 selectedFeaturesSet2Indexes.clear();
1766 interestPointsSet2RTree.
search( auxEnvelope,
1767 selectedFeaturesSet2Indexes );
1769 selectedFeaturesSet2IndexesSize = selectedFeaturesSet2Indexes.size();
1772 corrMatrixLinePtr = paramsPtr->
m_corrMatrixPtr->operator[]( feat1Idx );
1776 for(
unsigned int selectedFSIIdx = 0 ; selectedFSIIdx <
1777 selectedFeaturesSet2IndexesSize ; ++selectedFSIIdx )
1779 feat2Idx = selectedFeaturesSet2Indexes[ selectedFSIIdx ];
1785 for( featCol = 0 ; featCol < featureElementsNmb ; ++featCol )
1787 sumAA += feat1Ptr[ featCol ] * feat1Ptr[ featCol ];
1788 sumBB += feat2Ptr[ featCol ] * feat2Ptr[ featCol ];
1791 cc_norm = std::sqrt( sumAA * sumBB );
1793 if( cc_norm == 0.0 )
1795 corrMatrixLinePtr[ feat2Idx ] = 0;
1800 for( featCol = 0 ; featCol < featureElementsNmb ; ++featCol )
1802 ccorrelation += ( feat1Ptr[ featCol ] * feat2Ptr[ featCol ] ) /
1806 corrMatrixLinePtr[ feat2Idx ] = ccorrelation;
te::rp::TiePointsLocatorInputParameters m_inputParameters
Input parameters.
void clear()
Clear all allocated resources and go back to the initial default parameters.
std::vector< InterestPointsSetT > * m_interestPointsSubSectorsPtr
A pointer to a valid interest points container (one element by subsector)..
bool executeMatchingByCorrelation(const FloatsMatrix &featuresSet1, const FloatsMatrix &featuresSet2, const InterestPointsSetT &interestPointsSet1, const InterestPointsSetT &interestPointsSet2, te::gm::GeometricTransformation const *const raster1ToRaster2TransfPtr, const double raster1ToRaster2TransfDMapError, MatchedInterestPointsSetT &matchedPoints) const
Match each feature using correlation.
TiePointsLocatorMoravecStrategy()
void reset()
Clear all internal allocated resources and go back to the initial not-initialized state...
A class that represents an R-tree.
unsigned int m_y
Point Y coord.
FloatsMatrix const * m_featuresSet2Ptr
static bool applyMeanFilter(const FloatsMatrix &inputData, FloatsMatrix &outputData, const unsigned int iterationsNumber)
Mean Filter.
unsigned int * m_nextFeatureIdx1ToProcessPtr
bool initialize(const te::rp::TiePointsLocatorInputParameters &inputParameters)
Initialize the strategy.
FloatsMatrix const * m_rasterDataPtr
The loaded raster data.
This class can be used to inform the progress of a task.
double m_urx
Upper right corner x-coordinate.
unsigned int m_moravecWindowWidth
Thread return value pointer.
float m_feature
Matched interest feature.
~TiePointsLocatorMoravecStrategy()
FloatsMatrix const * m_featuresSet1Ptr
te::gm::GeometricTransformation const * m_raster1ToRaster2TransfPtr
A pointer to a transformation direct mapping raster 1 indexed coords into raster 2 indexed coords...
MemoryPolicy getMemPolicy() const
Returns the current memory policy.
unsigned int m_x
Point X coord.
#define TERP_TRUE_OR_RETURN_FALSE(value, message)
Checks if value is true. For false values a warning message will be logged and a return of context wi...
Tie-Pointsr locator Moravec strategy.
static bool generateCorrelationFeatures(const InterestPointsSetT &interestPoints, const unsigned int correlationWindowWidth, const FloatsMatrix &rasterData, FloatsMatrix &features, InterestPointsSetT &validInteresPoints)
Generate correlation features ( normalized - unit vector ) matrix for the given interes points...
InterestPointT const * m_interestPointsSet1Ptr
TECOMMONEXPORT unsigned int GetPhysProcNumber()
Returns the number of physical processors.
double m_llx
Lower left corner x-coordinate.
unsigned int m_processingBlocksNumber
The raster data will be splitted into this number of blocks for processing.
An Envelope defines a 2D rectangular region.
TECOMMONEXPORT unsigned long int GetUsedVirtualMemory()
Returns the amount of used virtual memory (bytes) for the current process (physical + swapped)...
InterestPointT const * m_interestPointsSet2Ptr
static void locateMoravecInterestPointsThreadEntry(MoravecLocatorThreadParams *paramsPtr)
Movavec locator thread entry.
unsigned int getColumnsNumber() const
The number of current matrix columns.
mydialect insert("+", new te::da::BinaryOpEncoder("+"))
std::multiset< InterestPointT > InterestPointsSetT
bool m_isInitialized
true if this instance is initialized.
int search(const te::gm::Envelope &mbr, std::vector< DATATYPE > &report) const
Range search query.
unsigned int m_maxInterestPointsBySubSector
The maximum number of interest points by sub-sector.
The parameters passed to the matchCorrelationEuclideanThreadEntry method.
static bool loadRasterData(te::rst::Raster const *rasterPtr, const std::vector< unsigned int > &rasterBands, te::rst::Raster const *maskRasterPtr, const unsigned int maskRasterBand, const unsigned int rasterTargetAreaLineStart, const unsigned int rasterTargetAreaColStart, const unsigned int rasterTargetAreaWidth, const unsigned int rasterTargetAreaHeight, const double desiredRescaleFactorX, const double desiredRescaleFactorY, const te::rst::Interpolator::Method rasterInterpMethod, const unsigned char maxMemPercentUsage, std::vector< boost::shared_ptr< FloatsMatrix > > &loadedRasterData, UCharsMatrix &loadedMaskRasterData, double &achievedRescaleFactorX, double &achievedRescaleFactorY)
Load rasters data (normalized between 0 and 1).
void reset()
Reset (clear) the active instance data.
double m_lly
Lower left corner y-coordinate.
bool getMatchedInterestPoints(te::gm::GeometricTransformation const *const raster1ToRaster2TransfPtr, const double raster1ToRaster2TransfDMapError, MatchedInterestPointsSetT &matchedInterestPoints)
Try to find matched interest points.
unsigned int getAutoMaxTiePointsNumber() const
Returns a automatically calculated optimum maximum amount tie-points following the current parameters...
bool locateMoravecInterestPoints(const unsigned int maxInterestPoints, const FloatsMatrix &rasterData, UCharsMatrix const *maskRasterDataPtr, InterestPointsSetT &interestPoints) const
Moravec interest points locator.
static void roolUpBuffer(BufferElementT **bufferPtr, const unsigned int &bufferLinesNumber)
RoolUp a buffer of lines.
The parameters passed to the moravecLocatorThreadEntry method.
void insert(const te::gm::Envelope &mbr, const DATATYPE &data)
It inserts an item into the tree.
double m_ury
Upper right corner y-coordinate.
unsigned int m_tiePointsSubSectorsSplitFactor
The number of sectors along each direction.
A generic template matrix.
TemplateElementType ElementTypeT
Public matrix element type definition.
InterestPointT m_point2
Interest point 2.
boost::mutex * m_interestPointsAccessMutexPtr
A pointer to a valid mutex to control the output interest points container access.
boost::mutex * m_syncMutexPtr
double m_searchOptTreeSearchRadius
Optimization tree search radius (pixels).
static void executeMatchingByCorrelationThreadEntry(ExecuteMatchingByCorrelationThreadEntryParams *paramsPtr)
Correlation/Euclidean match thread entry.
#define ROUND(x)
Minimum of two values.
UCharsMatrix const * m_maskRasterDataPtr
The loaded mask raster data pointer (or zero if no mask is avaliable).
TECOMMONEXPORT unsigned long int GetTotalPhysicalMemory()
Returns the amount of total physical memory (bytes).
unsigned int * m_nextRasterLinesBlockToProcessValuePtr
A pointer to a valid counter to control the blocks processing sequence.
unsigned int getLinesNumber() const
The number of current matrix lines.
std::multiset< MatchedInterestPointsT > MatchedInterestPointsSetT
#define TERP_TRUE_OR_THROW(value, message)
Checks if value is true and throws an exception if not.
InterestPointT m_point1
Interest point 1.
FloatsMatrix * m_corrMatrixPtr
TECOMMONEXPORT unsigned long int GetTotalVirtualMemory()
Returns the amount of total virtual memory (bytes) that can be claimed by the current process (physic...
float m_feature1
Interest point feature 1 value.
boost::mutex * m_rastaDataAccessMutexPtr
A pointer to a valid mutex to controle raster data access.