From 004a461ec61c87fa1e7e9eeaeb2d6c3856a573ff Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ga=C3=ABlle=20USSEGLIO?= <gaelle.usseglio@cnes.fr>
Date: Tue, 2 Jun 2020 14:50:01 +0000
Subject: [PATCH] ENH : Fix streaming into TilesAnalysisWith Step filter with
 buffer to store previous output region

---
 app/otbSARPhaseFiltering.cxx                  |  11 +-
 include/otbTilesAnalysisWithStepImageFilter.h |  15 ++
 .../otbTilesAnalysisWithStepImageFilter.txx   | 229 ++++++++++++++----
 3 files changed, 208 insertions(+), 47 deletions(-)

diff --git a/app/otbSARPhaseFiltering.cxx b/app/otbSARPhaseFiltering.cxx
index f7a863c..9307df6 100644
--- a/app/otbSARPhaseFiltering.cxx
+++ b/app/otbSARPhaseFiltering.cxx
@@ -226,7 +226,8 @@ namespace otb
 	functor->SetSizeTiles(sizeTiles);
 	functor->Setalpha(alpha);
 
-	// Instanciate the Coregistration filter
+	// Instanciate the Tiles Analysis (with Step) filter
+	// Limitation : Works only with a streaming by strips
 	TilesFilterType::Pointer filterPhase = TilesFilterType::New();
 	m_Ref.push_back(filterPhase.GetPointer());
 	filterPhase->SetSizeTiles(sizeTiles);
@@ -293,10 +294,10 @@ namespace otb
 	m_Ref.push_back(filterGroupedBy.GetPointer());
 	
 	// Only for phase and coherence
-	std::vector<std::string> bands;
-	bands.push_back("pha");
-	bands.push_back("coh");
-	filterGroupedBy->bandSelection(bands);
+	// std::vector<std::string> bands;
+	// bands.push_back("pha");
+	// bands.push_back("coh");
+	// filterGroupedBy->bandSelection(bands);
 
 	filterGroupedBy->SetInput(filterPhase->GetOutput());
 	filterGroupedBy->SetMLRan(factor_ran);
diff --git a/include/otbTilesAnalysisWithStepImageFilter.h b/include/otbTilesAnalysisWithStepImageFilter.h
index 1d62afe..bd2f65a 100644
--- a/include/otbTilesAnalysisWithStepImageFilter.h
+++ b/include/otbTilesAnalysisWithStepImageFilter.h
@@ -46,6 +46,8 @@ namespace otb
  * _ assign output with accumulations
  * 
  * This filter can be used for phase filtering.
+ *
+ * Limitations : This filter works only with a streaming by strips
  * \ingroup DiapOTBModule
  */
 
@@ -156,6 +158,8 @@ protected:
 
   void BeforeThreadedGenerateData();
 
+  void AfterThreadedGenerateData();
+
  /** Static function used as a "callback" by the MultiThreader.  The threading
   * library will call this routine for each thread, which will delegate the
   * control to ThreadedGenerateData(). */
@@ -183,6 +187,17 @@ protected:
 
   MutexType * m_Mutex;
 
+  // Store previous/effective index and size of output requested region and expect new ones
+  // Used into GenerateInputRequestedRegion to define output requested region
+  ImageOutIndexType m_expectedIndex;
+  ImageOutSizeType m_expectedSize;
+  ImageOutIndexType m_saveIndex;
+  ImageOutSizeType m_saveSize;
+
+  // Store previous output region with it size : to initialize correctly current output region
+  ImageOutSizeType m_oldSize;
+  ImageOutPixelType * m_saveOutputBuffer;
+
 };
 
 } // End namespace otb
diff --git a/include/otbTilesAnalysisWithStepImageFilter.txx b/include/otbTilesAnalysisWithStepImageFilter.txx
index be9861b..208da01 100644
--- a/include/otbTilesAnalysisWithStepImageFilter.txx
+++ b/include/otbTilesAnalysisWithStepImageFilter.txx
@@ -46,6 +46,8 @@ namespace otb
     this->m_Margin = 0;
 
     m_Mutex = new MutexType(); 
+
+    m_saveOutputBuffer = 0;
   }
 
   /** 
@@ -57,6 +59,12 @@ namespace otb
   {
     delete m_Mutex;
     m_Mutex = 0;
+
+    if (m_saveOutputBuffer)
+      {
+	delete m_saveOutputBuffer;
+	m_saveOutputBuffer = 0;
+      }
   }
  
 
@@ -126,27 +134,49 @@ namespace otb
   {
     ////// Set Output Requestion Region (Contruct Tiles with given size and step for output) ///////
     ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
-    
-    //std::cout << "Before outputRequestedRegion : " << outputRequestedRegion << std::endl;
 
-    // New Index : closest inf multiple of Step
+    // New Index : closest inf multiple of Step and shift to this->m_SizeTiles - m_Step 
     ImageOutIndexType outputIndex;
-    outputIndex[0] = (static_cast<int>(outputRequestedRegion.GetIndex()[0]/this->m_Step))*
-      this->m_Step;
+
     
-    outputIndex[1] = (static_cast<int>(outputRequestedRegion.GetIndex()[1]/this->m_Step))*
-      this->m_Step;
+    if (outputRequestedRegion.GetIndex()[0]/this->m_Step < (this->m_SizeTiles/m_Step))
+      {
+	outputIndex[0] = outputRequestedRegion.GetIndex()[0];
 
-    // New Size : closest sup multiple of SizeTiles with addtional when difference between new and 
-    // old index is too important
+	m_expectedIndex[0] = outputIndex[0]; 
+      }
+    else
+      {
+	outputIndex[0] = (static_cast<int>(outputRequestedRegion.GetIndex()[0]/this->m_Step))*
+	  this->m_Step - (this->m_SizeTiles - m_Step) ;
+
+	m_expectedIndex[0] = m_saveIndex[0] + m_saveSize[0] - (this->m_SizeTiles - m_Step); 
+      }
+
+    if (outputRequestedRegion.GetIndex()[1]/this->m_Step < (this->m_SizeTiles/m_Step))
+      {
+	outputIndex[1] = outputRequestedRegion.GetIndex()[1];
+	m_expectedIndex[1] = outputIndex[1];
+      }
+    else
+      {
+	outputIndex[1] = (static_cast<int>(outputRequestedRegion.GetIndex()[1]/this->m_Step))*
+	 this->m_Step - (this->m_SizeTiles - m_Step);
+
+	m_expectedIndex[1] = m_saveIndex[1] + m_saveSize[1] - (this->m_SizeTiles - m_Step);
+      }
+
+
+    // New Size : closest sup multiple of Step with addtional when difference between new and 
+    // old index and SizeTiles to handles the last tiles
     ImageOutSizeType outputSize;
+
     outputSize[0] = (static_cast<int>((outputRequestedRegion.GetIndex()[0] - outputIndex[0] +
-				       outputRequestedRegion.GetSize()[0])/this->m_SizeTiles) + 
-		                       1)*this->m_SizeTiles + m_Step;
+				       outputRequestedRegion.GetSize()[0])/this->m_Step))*
+                                       this->m_Step + this->m_SizeTiles;
     outputSize[1] = (static_cast<int>((outputRequestedRegion.GetIndex()[1] - outputIndex[1] +
-				       outputRequestedRegion.GetSize()[1])/this->m_SizeTiles) + 
-		                       1)*this->m_SizeTiles + m_Step;
-    
+				       outputRequestedRegion.GetSize()[1])/this->m_Step))*
+		                       this->m_Step + this->m_SizeTiles;
 
     // Region affectation (with Crop to respect image dimension)
     outputRequestedRegion.SetSize(outputSize);
@@ -154,13 +184,133 @@ namespace otb
     outputRequestedRegion.Crop(this->GetOutput()->GetLargestPossibleRegion());
     this->GetOutput()->SetRequestedRegion(outputRequestedRegion);
 
-    //std::cout << "After outputRequestedRegion : " << outputRequestedRegion << std::endl;
 
-    ///////////// With the new output requested region, find the region into input image /////////////
+    // input requested region same than output
     ImageInRegionType inputRequestedRegion = outputRequestedRegion;
     ImageInPointer  inputPtr = const_cast< ImageInType * >( this->GetImageInput() );
     inputPtr->SetRequestedRegion(inputRequestedRegion);
 
+    // Store index and size (for next call of this function)
+    m_saveIndex[0] = outputIndex[0]; 
+    m_saveIndex[1] = outputIndex[1];
+    m_saveSize[0] = outputSize[0]; 
+    m_saveSize[1] = outputSize[1];
+  }
+
+  /**
+   * Method BeforeThreadedGenerateData
+   */
+  template<class TImageIn, class TImageOut, class TFunction>
+  void
+  TilesAnalysisWithStepImageFilter< TImageIn, TImageOut, TFunction >
+ ::BeforeThreadedGenerateData()
+  {
+    // Typedef for Iterators on selected region
+    typedef itk::ImageRegionIterator<ImageOutType> OutItType;
+
+    // Initialize output region 
+    ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
+
+    
+    OutItType OutIt(this->GetOutput(), outputRequestedRegion);
+    OutIt.GoToBegin();
+
+    ImageOutIndexType first_index = OutIt.GetIndex();
+   
+    // Limitations : Works only with a streaing by strips
+    // TODO : See to allow others kind of streamings
+    int counter = (m_oldSize[1] - (this->m_SizeTiles - m_Step + 
+				   m_expectedIndex[1]-m_saveIndex[1])) * m_oldSize[0];
+
+
+    while (!OutIt.IsAtEnd())
+      {
+	ImageOutIndexType index_current = OutIt.GetIndex();
+	
+	// First Streaming
+	if (first_index[0] < (this->m_SizeTiles-m_Step) && 
+	    first_index[1] < (this->m_SizeTiles-m_Step))
+	  {
+	    OutIt.Set(0);
+	  }
+	else
+	  {
+	    // First pixels have already been initialized or estimated => saveBuffer 
+	    // and 0 for others
+	    // TODO : See to allow another kind of streaming (not only strips)
+	    if ((index_current[1]) >= (m_expectedIndex[1] + this->m_SizeTiles - m_Step)) 
+	      {
+		// Init for accumulations
+		OutIt.Set(0);
+	      }
+	    else
+	      {	
+		// Keep values from last streaming
+		OutIt.Set(m_saveOutputBuffer[counter]);
+		
+		counter++; 
+	      }
+
+	  }
+	
+	++OutIt;
+      }
+
+    // Free memory
+    delete m_saveOutputBuffer;
+    m_saveOutputBuffer = 0;
+  }
+
+
+  /**
+   * Method AfterThreadedGenerateData
+   */
+  template<class TImageIn, class TImageOut, class TFunction>
+  void
+  TilesAnalysisWithStepImageFilter< TImageIn, TImageOut, TFunction >
+  ::AfterThreadedGenerateData()
+  {
+    ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
+
+    // Store all output region 
+    ImageOutIndexType saveOutputIndex;
+    saveOutputIndex[0] = outputRequestedRegion.GetIndex()[0];
+    saveOutputIndex[1] = outputRequestedRegion.GetIndex()[1];
+
+    ImageOutSizeType saveOutputSize;
+    saveOutputSize[0] = outputRequestedRegion.GetSize()[0];
+    saveOutputSize[1] = outputRequestedRegion.GetSize()[1];
+
+    ImageOutRegionType saveOutputRegion;
+    saveOutputRegion.SetSize(saveOutputSize);
+    saveOutputRegion.SetIndex(saveOutputIndex);
+    saveOutputRegion.Crop(this->GetOutput()->GetLargestPossibleRegion());
+
+    m_oldSize[0] = saveOutputSize[0];
+    m_oldSize[1] = saveOutputSize[1];
+    
+    // Typedef for Iterators on selected region
+    typedef itk::ImageRegionIterator<ImageOutType> OutItType;
+
+    OutItType OutIt(this->GetOutput(), saveOutputRegion);
+    OutIt.GoToBegin();
+
+
+    m_saveOutputBuffer = new ImageOutPixelType[outputRequestedRegion.GetSize()[0]*
+					       outputRequestedRegion.GetSize()[1]];
+
+
+    // Loop for assignations
+    int counter = 0;
+    while (!OutIt.IsAtEnd())
+      {
+	
+	m_saveOutputBuffer[counter] = OutIt.Get();
+
+	++counter; 
+	++OutIt;
+      }
+    
   }
 
 
@@ -181,9 +331,10 @@ namespace otb
     ImageOutRegionType outputRegion = this->GetOutput()->GetRequestedRegion();
             
     // Define the number of tiles for processing
-    int nbTiles_Lines = outputRegion.GetSize()[1]/m_Step - this->m_SizeTiles/m_Step; 
-    int nbTiles_Col = outputRegion.GetSize()[0]/m_Step - this->m_SizeTiles/m_Step;
+    int nbTiles_Lines = (outputRegion.GetSize()[1]-this->m_SizeTiles)/m_Step; 
+    int nbTiles_Col = (outputRegion.GetSize()[0]-this->m_SizeTiles)/m_Step;
     
+
     // Check if additionnal tiles are required
     if (outputRegion.GetSize()[1] % m_Step != 0)  
       {
@@ -195,9 +346,19 @@ namespace otb
       }
     
     // First Tile = index of output requested region
-    int IndL_FirstTiles = outputRegion.GetIndex()[1];
-    int IndC_FirstTiles = outputRegion.GetIndex()[0];
+    //int IndL_FirstTiles = outputRegion.GetIndex()[1];
+    //int IndC_FirstTiles = outputRegion.GetIndex()[0];
+    
+    // First Tile = index of expected output requested region (real output region might be larger 
+    // than expected) => do not process all tiles because some of its had been estimated.
+    int IndL_FirstTiles = m_expectedIndex[1];
+    int IndC_FirstTiles = m_expectedIndex[0];
      
+    nbTiles_Lines -=  (m_expectedIndex[1] - outputRegion.GetIndex()[1])/m_Step;
+    nbTiles_Col -=  (m_expectedIndex[0] - outputRegion.GetIndex()[0])/m_Step;
+
+    nbTiles_Lines++;
+
     ///////////////////////// Pre Processing ///////////////////////
     // Call Function initialization
     this->m_FunctionOnTiles->Initialize(nbThreads);
@@ -222,33 +383,13 @@ namespace otb
     this->GetMultiThreader()->SingleMethodExecute();
     
     ///////////// Post Processing ////////////
+    this->AfterThreadedGenerateData();
+    
     // Call Function terminate
     this->m_FunctionOnTiles->Terminate();
   }
 
 
-  template<class TImageIn, class TImageOut, class TFunction>
-  void
-  TilesAnalysisWithStepImageFilter< TImageIn, TImageOut, TFunction >
- ::BeforeThreadedGenerateData()
-  {
-    // Typedef for Iterators on selected region
-    typedef itk::ImageRegionIterator<ImageOutType> OutItType;
-
-    // Initialize output region 
-    ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
-
-    
-    OutItType OutIt(this->GetOutput(), outputRequestedRegion);
-    OutIt.GoToBegin();
-    while (!OutIt.IsAtEnd())
-      {
-	// For accumulations
-	OutIt.Set(0);
-	++OutIt;
-      }
-  }
-
 
  template<class TImageIn, class TImageOut, class TFunction>
   void
@@ -262,6 +403,7 @@ namespace otb
    // Typedef for Iterators on selected region
    typedef itk::ImageRegionIterator<ImageOutType> OutItType;
 
+
    for (unsigned int k = startTileId; k < stopTileId; k++)  //<= ??????
      {
        // Transform Tiles id into index
@@ -299,6 +441,7 @@ namespace otb
        int compteur = 0;
 
        // Assignate Output with function result
+       // Mutex to guarantee Thread-Safety
        // Iterators on output
        OutItType OutIt(this->GetOutput(), outputCurrentRegion);
        OutIt.GoToBegin();
@@ -308,6 +451,7 @@ namespace otb
 	   // Accumulations
 	   ImageOutPixelType outTmp = OutIt.Get();
 	   OutIt.Set(outTmp + TilesCache[compteur]);
+
 	   ++compteur; 
 	   ++OutIt;
 	 }
@@ -350,6 +494,7 @@ namespace otb
 	  static_cast<unsigned int>((threadId+1) * 
 				    std::ceil(static_cast<double>(total)/(static_cast<double>(threadCount))));
 
+
 	if (stop > total) stop = total;
 
 	// For very small graphs it might occur that start = stop. In this
-- 
GitLab