From 124c0a42f39e630971c3a99dbcdb919a20651519 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ga=C3=ABlle=20USSEGLIO?= <gaelle.usseglio@cnes.fr>
Date: Fri, 6 Nov 2020 11:01:29 +0000
Subject: [PATCH] ENH : Adapt rough shifts estimation to decrease memory
 consumption

---
 app/otbSARCorrelationGrid.cxx                 | 86 ++++++++++++++-----
 app/otbSARCorrelationRough.cxx                | 81 +++++++++++------
 ...bSARTemporalCorrelationGridImageFilter.txx | 42 +++++----
 3 files changed, 144 insertions(+), 65 deletions(-)

diff --git a/app/otbSARCorrelationGrid.cxx b/app/otbSARCorrelationGrid.cxx
index 27bf362..b8888f3 100644
--- a/app/otbSARCorrelationGrid.cxx
+++ b/app/otbSARCorrelationGrid.cxx
@@ -24,6 +24,8 @@
 #include "itkMacro.h"
 #include "itkFFTNormalizedCorrelationImageFilter.h"
 
+#include "otbMultiToMonoChannelExtractROI.h"
+
 #include "otbSARStreamingMaximumMinimumImageFilter.h"
 #include "otbSARTemporalCorrelationGridImageFilter.h"
 
@@ -49,7 +51,8 @@ public:
   typedef itk::FFTNormalizedCorrelationImageFilter<FloatImageType, FloatImageType>          CorFilterType;
   typedef otb::SARStreamingMaximumMinimumImageFilter< FloatImageType>		   MinMaxFilterType;
   typedef otb::SARTemporalCorrelationGridImageFilter<FloatImageType, FloatVectorImageType>  CorGridFilterType;
-  
+  typedef otb::MultiToMonoChannelExtractROI<FloatVectorImageType::InternalPixelType,
+					    FloatImageType::InternalPixelType> ExtractROIFilterType;
 
 private:
   void DoInit() override
@@ -137,40 +140,67 @@ void DoExecute() override
   otbAppLogINFO(<<"Grid Step for azimut : "<<grid_step_azi);
  
   // Get master and slave image
-  FloatImageType::Pointer MasterPtr = GetParameterFloatImage("inmaster");
-  FloatImageType::Pointer SlavePtr = GetParameterFloatImage("inslave");
+  FloatVectorImageType::Pointer MasterPtr = GetParameterFloatVectorImage("inmaster");
+  FloatVectorImageType::Pointer SlavePtr = GetParameterFloatVectorImage("inslave");
    
-  // Check input size (Master size) to protect memory. 
-  // If size > 15000 raise an exception and indicate more appropriate ml factors
-  if (MasterPtr->GetLargestPossibleRegion().GetSize()[0] > 15000 || 
-      MasterPtr->GetLargestPossibleRegion().GetSize()[1] > 15000)
-    {
-      // integer divisions
-      int intermediate_mlran = MasterPtr->GetLargestPossibleRegion().GetSize()[0]/15000 + 1;
-      int intermediate_mlazi = MasterPtr->GetLargestPossibleRegion().GetSize()[1]/15000 + 1;
+  // Initialize values 
+  float shiftML_range = 0;
+  float shiftML_azimut = 0;
+  int startx = 0;
+  int starty = 0;
+  int sizex = std::min(MasterPtr->GetLargestPossibleRegion().GetSize()[0],
+		       SlavePtr->GetLargestPossibleRegion().GetSize()[0]);
+  int sizey = std::min(MasterPtr->GetLargestPossibleRegion().GetSize()[1],
+		       SlavePtr->GetLargestPossibleRegion().GetSize()[1]);
 
-      otbAppLogCRITICAL(<<"ML Factors are not appropriate for these estimations. Please use other factors such as : "<< intermediate_mlran << " x " << intermediate_mlazi );
+  
+  // Extract ROI from Master and Slave Image to estimate global shifts without a too large amount
+  // of RAM => size = 3000*3000 at the center
+  if (sizex > 3000 && 
+      sizey > 3000)
+    {
+      startx = MasterPtr->GetLargestPossibleRegion().GetSize()[0]/2 - 1500;
+      starty = MasterPtr->GetLargestPossibleRegion().GetSize()[1]/2 - 1500;
 
-      itkExceptionMacro(<<"ML Factors are not appropriate for these estimations");
+      sizex = 3000;
+      sizey = 3000;
     }
 
+  ExtractROIFilterType::Pointer extractROIFilter_master = ExtractROIFilterType::New();
+  m_Ref.push_back(extractROIFilter_master.GetPointer());
+  extractROIFilter_master->SetInput(MasterPtr);
+  extractROIFilter_master->SetChannel(1);
+  extractROIFilter_master->SetStartX(startx);
+  extractROIFilter_master->SetStartY(starty);
+  extractROIFilter_master->SetSizeX(sizex);
+  extractROIFilter_master->SetSizeY(sizey);
+  
+  ExtractROIFilterType::Pointer extractROIFilter_slave = ExtractROIFilterType::New();
+  m_Ref.push_back(extractROIFilter_slave.GetPointer());
+  extractROIFilter_slave->SetInput(SlavePtr);
+  extractROIFilter_slave->SetChannel(1);
+  extractROIFilter_slave->SetStartX(startx);
+  extractROIFilter_slave->SetStartY(starty);
+  extractROIFilter_slave->SetSizeX(sizex);
+  extractROIFilter_slave->SetSizeY(sizey);
+
   // Correlation Filter
   CorFilterType::Pointer correlationFilter = CorFilterType::New();
   m_Ref.push_back(correlationFilter.GetPointer());
-  correlationFilter->SetFixedImage(MasterPtr);
-  correlationFilter->SetMovingImage(SlavePtr);
+  correlationFilter->SetFixedImage(extractROIFilter_master->GetOutput());
+  correlationFilter->SetMovingImage(extractROIFilter_slave->GetOutput());
   
   MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
   minMaxFilter->SetInput(correlationFilter->GetOutput());
   // Adapt streaming with ram parameter (default 256 MB)
   minMaxFilter->GetStreamer()->SetAutomaticStrippedStreaming(GetParameterInt("ram"));
   minMaxFilter->Update();
-
-  float shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
-			  minMaxFilter->GetIndexOfMax()[0]);
-  float shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
-			   -minMaxFilter->GetIndexOfMax()[1]);
-
+  
+  shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
+		   minMaxFilter->GetIndexOfMax()[0]);
+  shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
+		    -minMaxFilter->GetIndexOfMax()[1]);
+  
   // Correlation Grid Filter 
   CorGridFilterType::Pointer filterCorrelationGrid = CorGridFilterType::New();
   m_Ref.push_back(filterCorrelationGrid.GetPointer());
@@ -188,8 +218,18 @@ void DoExecute() override
 
   
   // Define the main pipeline 
-  filterCorrelationGrid->SetMasterInput(MasterPtr);
-  filterCorrelationGrid->SetSlaveInput(SlavePtr);
+  ExtractROIFilterType::Pointer extractROIFilter_masterfull = ExtractROIFilterType::New();
+  m_Ref.push_back(extractROIFilter_masterfull.GetPointer());
+  extractROIFilter_masterfull->SetInput(MasterPtr);
+  extractROIFilter_masterfull->SetChannel(1);
+  ExtractROIFilterType::Pointer extractROIFilter_slavefull = ExtractROIFilterType::New();
+  m_Ref.push_back(extractROIFilter_slavefull.GetPointer());
+  extractROIFilter_slavefull->SetInput(SlavePtr);
+  extractROIFilter_slavefull->SetChannel(1);
+
+  
+  filterCorrelationGrid->SetMasterInput(extractROIFilter_masterfull->GetOutput());
+  filterCorrelationGrid->SetSlaveInput(extractROIFilter_slavefull->GetOutput());
   SetParameterOutputImage("out", filterCorrelationGrid->GetOutput());
 }
    // Vector for filters 
diff --git a/app/otbSARCorrelationRough.cxx b/app/otbSARCorrelationRough.cxx
index 36f6668..13b3914 100644
--- a/app/otbSARCorrelationRough.cxx
+++ b/app/otbSARCorrelationRough.cxx
@@ -24,6 +24,8 @@
 #include "itkMacro.h"
 #include "itkFFTNormalizedCorrelationImageFilter.h"
 
+#include "otbMultiToMonoChannelExtractROI.h"
+
 #include "otbSARStreamingMaximumMinimumImageFilter.h"
 
 #include <iostream>
@@ -47,6 +49,8 @@ public:
   // Filters
   typedef itk::FFTNormalizedCorrelationImageFilter<FloatImageType, FloatImageType> CorFilterType;
   typedef otb::SARStreamingMaximumMinimumImageFilter< FloatImageType>		   MinMaxFilterType;
+  typedef otb::MultiToMonoChannelExtractROI<FloatVectorImageType::InternalPixelType,
+					    FloatImageType::InternalPixelType> ExtractROIFilterType;
 
 private:
   void DoInit() override
@@ -129,46 +133,73 @@ void DoExecute() override
   otbAppLogINFO(<<"RAM : "<< GetParameterInt("ram"));
 
   // Get master and slave image
-  FloatImageType::Pointer MasterPtr = GetParameterFloatImage("inmaster");
-  FloatImageType::Pointer SlavePtr = GetParameterFloatImage("inslave");
+  FloatVectorImageType::Pointer MasterPtr = GetParameterFloatVectorImage("inmaster");
+  FloatVectorImageType::Pointer SlavePtr = GetParameterFloatVectorImage("inslave");
+   
+  // Initialize values 
+  float shiftML_range = 0;
+  float shiftML_azimut = 0;
+  float shiftSLC_range = 0;
+  float shiftSLC_azimut = 0;
+  int startx = 0;
+  int starty = 0;
+  int sizex = std::min(MasterPtr->GetLargestPossibleRegion().GetSize()[0],
+		       SlavePtr->GetLargestPossibleRegion().GetSize()[0]);
+  int sizey = std::min(MasterPtr->GetLargestPossibleRegion().GetSize()[1],
+		       SlavePtr->GetLargestPossibleRegion().GetSize()[1]);
+
   
-  // Check input size (Master size) to protect memory. 
-  // If size > 15000 raise an exception and indicate more appropriate ml factors
-  if (MasterPtr->GetLargestPossibleRegion().GetSize()[0] > 15000 || 
-      MasterPtr->GetLargestPossibleRegion().GetSize()[1] > 15000)
+  // Extract ROI from Master and Slave Image to estimate global shifts without a too large amount
+  // of RAM => size = 3000*3000 at the center
+  if (sizex > 3000 && 
+      sizey > 3000)
     {
-      // integer divisions
-      int intermediate_mlran = MasterPtr->GetLargestPossibleRegion().GetSize()[0]/15000 + 1;
-      int intermediate_mlazi = MasterPtr->GetLargestPossibleRegion().GetSize()[1]/15000 + 1;
+      startx = MasterPtr->GetLargestPossibleRegion().GetSize()[0]/2 - 1500;
+      starty = MasterPtr->GetLargestPossibleRegion().GetSize()[1]/2 - 1500;
 
-      otbAppLogCRITICAL(<<"ML Factors are not appropriate for these estimations. Please use other factors such as : "<< intermediate_mlran << " x " << intermediate_mlazi );
-
-      itkExceptionMacro(<<"ML Factors are not appropriate for these estimations");
+      sizex = 3000;
+      sizey = 3000;
     }
- 
+
+  ExtractROIFilterType::Pointer extractROIFilter_master = ExtractROIFilterType::New();
+  m_Ref.push_back(extractROIFilter_master.GetPointer());
+  extractROIFilter_master->SetInput(MasterPtr);
+  extractROIFilter_master->SetChannel(1);
+  extractROIFilter_master->SetStartX(startx);
+  extractROIFilter_master->SetStartY(starty);
+  extractROIFilter_master->SetSizeX(sizex);
+  extractROIFilter_master->SetSizeY(sizey);
   
+  ExtractROIFilterType::Pointer extractROIFilter_slave = ExtractROIFilterType::New();
+  m_Ref.push_back(extractROIFilter_slave.GetPointer());
+  extractROIFilter_slave->SetInput(SlavePtr);
+  extractROIFilter_slave->SetChannel(1);
+  extractROIFilter_slave->SetStartX(startx);
+  extractROIFilter_slave->SetStartY(starty);
+  extractROIFilter_slave->SetSizeX(sizex);
+  extractROIFilter_slave->SetSizeY(sizey);
+
   // Correlation Filter
   CorFilterType::Pointer correlationFilter = CorFilterType::New();
   m_Ref.push_back(correlationFilter.GetPointer());
-  correlationFilter->SetFixedImage(MasterPtr);
-  correlationFilter->SetMovingImage(SlavePtr);
-
+  correlationFilter->SetFixedImage(extractROIFilter_master->GetOutput());
+  correlationFilter->SetMovingImage(extractROIFilter_slave->GetOutput());
+  
   MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
   minMaxFilter->SetInput(correlationFilter->GetOutput());
   // Adapt streaming with ram parameter (default 256 MB)
   minMaxFilter->GetStreamer()->SetAutomaticStrippedStreaming(GetParameterInt("ram"));
   minMaxFilter->Update();
 
+  shiftSLC_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
+		    minMaxFilter->GetIndexOfMax()[0]) * static_cast<float>(factorML_ran);
+  shiftSLC_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
+		     -minMaxFilter->GetIndexOfMax()[1]) * static_cast<float>(factorML_azi);
 
-  float shiftSLC_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
-		       minMaxFilter->GetIndexOfMax()[0]) * static_cast<float>(factorML_ran);
-  float shiftSLC_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
-			-minMaxFilter->GetIndexOfMax()[1]) * static_cast<float>(factorML_azi);
-
-  float shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
-			  minMaxFilter->GetIndexOfMax()[0]);
-  float shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
-			   -minMaxFilter->GetIndexOfMax()[1]);
+  shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2.)-
+		   minMaxFilter->GetIndexOfMax()[0]);
+  shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2.)
+		    -minMaxFilter->GetIndexOfMax()[1]);
 
   // Assigne Outputs
   SetParameterFloat("shiftranslc", shiftSLC_range);
diff --git a/include/otbSARTemporalCorrelationGridImageFilter.txx b/include/otbSARTemporalCorrelationGridImageFilter.txx
index 1e67691..a341f17 100644
--- a/include/otbSARTemporalCorrelationGridImageFilter.txx
+++ b/include/otbSARTemporalCorrelationGridImageFilter.txx
@@ -231,12 +231,15 @@ namespace otb
     int firstL_master = 0;
     if (m_FirstLCOffset)
       {
-	int firstC = std::max (1 + m_WinCor_ran, 1 + m_WinCor_ran - m_RoughShift_ran + m_WinAround_ShiftRan);
-	int firstL = std::max (1 + m_WinCor_azi, 1 + m_WinCor_azi - m_RoughShift_azi + m_WinAround_ShiftAzi);
-	int lastC = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[0] - m_RoughShift_ran - 
+	int abs_RoughShift_ran = std::abs(m_RoughShift_ran);
+	int abs_RoughShift_azi = std::abs(m_RoughShift_azi);
+	
+	int firstC = std::max (1 + m_WinCor_ran, 1 + m_WinCor_ran - abs_RoughShift_ran + m_WinAround_ShiftRan);
+	int firstL = std::max (1 + m_WinCor_azi, 1 + m_WinCor_azi - abs_RoughShift_azi + m_WinAround_ShiftAzi);
+	int lastC = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[0] - abs_RoughShift_ran - 
 			     m_WinAround_ShiftRan - m_WinCor_ran, 
 			     masterPtr->GetLargestPossibleRegion().GetSize()[0]- m_WinCor_ran);
-	int lastL = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[1] - m_RoughShift_azi - 
+	int lastL = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[1] - abs_RoughShift_azi - 
 			     m_WinAround_ShiftAzi- m_WinCor_azi, 
 			     masterPtr->GetLargestPossibleRegion().GetSize()[1] - m_WinCor_azi);
     
@@ -328,15 +331,17 @@ namespace otb
     
     int firstC_Grid = 0;
     int firstL_Grid = 0;
-
+    int abs_RoughShift_ran = std::abs(m_RoughShift_ran);
+    int abs_RoughShift_azi = std::abs(m_RoughShift_azi);
+    
      if (m_FirstLCOffset)
        {
-	 int firstC = std::max (1 + m_WinCor_ran, 1 + m_WinCor_ran - m_RoughShift_ran + m_WinAround_ShiftRan);
-	 int firstL = std::max (1 + m_WinCor_azi, 1 + m_WinCor_azi - m_RoughShift_azi + m_WinAround_ShiftAzi);
-	 int lastC = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[0] - m_RoughShift_ran - 
+	 int firstC = std::max (1 + m_WinCor_ran, 1 + m_WinCor_ran - abs_RoughShift_ran + m_WinAround_ShiftRan);
+	 int firstL = std::max (1 + m_WinCor_azi, 1 + m_WinCor_azi - abs_RoughShift_azi + m_WinAround_ShiftAzi);
+	 int lastC = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[0] - abs_RoughShift_ran - 
 			      m_WinAround_ShiftRan - m_WinCor_ran, 
 			      masterPtr->GetLargestPossibleRegion().GetSize()[0]- m_WinCor_ran);
-	 int lastL = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[1] - m_RoughShift_azi - 
+	 int lastL = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[1] - abs_RoughShift_azi - 
 			      m_WinAround_ShiftAzi- m_WinCor_azi, 
 			      masterPtr->GetLargestPossibleRegion().GetSize()[1] - m_WinCor_azi);
 	 
@@ -390,11 +395,11 @@ namespace otb
     ImageInSizeType slaveSize;
 
     // Set the max_shift (margin)
-    firstL = masterRequestedIndex[1] - m_RoughShift_azi - m_WinAround_ShiftAzi - m_WinCor_azi;
-    firstC = masterRequestedIndex[0] - m_RoughShift_ran - m_WinAround_ShiftRan - m_WinCor_ran;
-    lastL = masterRequestedIndex[1] + masterRequestedSize[1] + m_RoughShift_azi + m_WinAround_ShiftAzi + 
+    firstL = masterRequestedIndex[1] - abs_RoughShift_azi - m_WinAround_ShiftAzi - m_WinCor_azi;
+    firstC = masterRequestedIndex[0] - abs_RoughShift_ran - m_WinAround_ShiftRan - m_WinCor_ran;
+    lastL = masterRequestedIndex[1] + masterRequestedSize[1] + abs_RoughShift_azi + m_WinAround_ShiftAzi + 
       m_WinCor_azi;
-    lastC =  masterRequestedIndex[0] + masterRequestedSize[0] + m_RoughShift_ran + m_WinAround_ShiftRan + m_WinCor_ran;
+    lastC =  masterRequestedIndex[0] + masterRequestedSize[0] + abs_RoughShift_ran + m_WinAround_ShiftRan + m_WinCor_ran;
 
     // Check the limits
     if (firstC < 0)
@@ -518,12 +523,15 @@ namespace otb
 
     if (m_FirstLCOffset)
        {
-	 int firstC = std::max (1 + m_WinCor_ran, 1 + m_WinCor_ran - m_RoughShift_ran + m_WinAround_ShiftRan);
-	 int firstL = std::max (1 + m_WinCor_azi, 1 + m_WinCor_azi - m_RoughShift_azi + m_WinAround_ShiftAzi);
-	 int lastC = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[0] - m_RoughShift_ran - 
+	 int abs_RoughShift_ran = std::abs(m_RoughShift_ran);
+	 int abs_RoughShift_azi = std::abs(m_RoughShift_azi);
+	 
+	 int firstC = std::max (1 + m_WinCor_ran, 1 + m_WinCor_ran - abs_RoughShift_ran + m_WinAround_ShiftRan);
+	 int firstL = std::max (1 + m_WinCor_azi, 1 + m_WinCor_azi - abs_RoughShift_azi + m_WinAround_ShiftAzi);
+	 int lastC = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[0] - abs_RoughShift_ran - 
 			      m_WinAround_ShiftRan - m_WinCor_ran, 
 			      masterPtr->GetLargestPossibleRegion().GetSize()[0]- m_WinCor_ran);
-	 int lastL = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[1] - m_RoughShift_azi - 
+	 int lastL = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[1] - abs_RoughShift_azi - 
 			      m_WinAround_ShiftAzi- m_WinCor_azi, 
 			      masterPtr->GetLargestPossibleRegion().GetSize()[1] - m_WinCor_azi);
 	 
-- 
GitLab