From 13ad2a305086ca488e8852238a58e685d4038c43 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ga=C3=ABlle=20USSEGLIO?= <gaelle.usseglio@cnes.fr>
Date: Wed, 13 May 2020 15:55:18 +0000
Subject: [PATCH] WIP: New application for phase filtering (part 2 : Add new
output estimation with a step for window extractions into inherited class
TilesAnalysisWithStep)
---
include/otbSARTilesPhaseFilteringFunctor.h | 11 +-
include/otbTilesAnalysisImageFilter.h | 16 +-
include/otbTilesAnalysisWithStepImageFilter.h | 43 ++-
.../otbTilesAnalysisWithStepImageFilter.txx | 277 +++++++++++++++++-
4 files changed, 335 insertions(+), 12 deletions(-)
diff --git a/include/otbSARTilesPhaseFilteringFunctor.h b/include/otbSARTilesPhaseFilteringFunctor.h
index f353588..873c0fe 100644
--- a/include/otbSARTilesPhaseFilteringFunctor.h
+++ b/include/otbSARTilesPhaseFilteringFunctor.h
@@ -243,6 +243,8 @@ public:
{
w[i] = 2.0*(i+1)/(float)m_SizeTiles;
w[m_SizeTiles-i-1]=w[i];
+ //w[i] = 1;
+ //w[m_SizeTiles-i-1]=w[i];
}
for (unsigned int i=0; i < m_SizeTiles; i++)
@@ -309,6 +311,9 @@ public:
{
peakAmp = bufAmp[k];
}
+
+ FFTInvIn[k][0] = 0;
+ FFTInvIn[k][1] = 0;
}
if (peakAmp != 0)
@@ -342,8 +347,10 @@ public:
if (norm != 0)
{
- TilesResult[k].real( (FFTInvOut[k][0]*wf[i1][i2])/norm );
- TilesResult[k].imag( (FFTInvOut[k][1]*wf[i1][i2])/norm );
+ //TilesResult[k].real( (FFTInvOut[k][0]*wf[i1][i2])/norm );
+ //TilesResult[k].imag( (FFTInvOut[k][1]*wf[i1][i2])/norm );
+ TilesResult[k].real( (FFTInvOut[k][0]*wf[i1][i2]) );
+ TilesResult[k].imag( (FFTInvOut[k][1]*wf[i1][i2]) );
}
else
{
diff --git a/include/otbTilesAnalysisImageFilter.h b/include/otbTilesAnalysisImageFilter.h
index 2505d7b..bf0edef 100644
--- a/include/otbTilesAnalysisImageFilter.h
+++ b/include/otbTilesAnalysisImageFilter.h
@@ -232,6 +232,15 @@ protected:
// Function to apply on Output Tiles
FunctionPointer m_FunctionOnTiles;
+
+ // Nb Tiles for each Pipeline (each call to GenerateData)
+ unsigned int m_NbTiles;
+ unsigned int m_NbTilesCol;
+
+ // Index for the first tile of each Pipeline (each call to GenerateData)
+ unsigned int m_IndLFirstTile;
+ unsigned int m_IndCFirstTile;
+
private:
TilesAnalysisImageFilter(const Self&); // purposely not implemented
@@ -254,13 +263,6 @@ protected:
// Step (SLC or SAR geo) for input deformation grid
GridSizeType m_GridStep;
- // Nb Tiles for each Pipeline (each call to GenerateData)
- unsigned int m_NbTiles;
- unsigned int m_NbTilesCol;
-
- // Index for the first tile of each Pipeline (each call to GenerateData)
- unsigned int m_IndLFirstTile;
- unsigned int m_IndCFirstTile;
};
} // End namespace otb
diff --git a/include/otbTilesAnalysisWithStepImageFilter.h b/include/otbTilesAnalysisWithStepImageFilter.h
index 2af29d5..1d62afe 100644
--- a/include/otbTilesAnalysisWithStepImageFilter.h
+++ b/include/otbTilesAnalysisWithStepImageFilter.h
@@ -25,6 +25,7 @@
#include "itkSmartPointer.h"
#include "itkPoint.h"
+#include "itkSimpleFastMutexLock.h"
#include "itkImageScanlineConstIterator.h"
#include "itkImageScanlineIterator.h"
@@ -36,7 +37,7 @@
namespace otb
{
/** \class TilesAnalysisWithStepImageFilter
- * \brief Creates/Analyses tiles fixed by a given size and a input step.
+ * \brief Creates/Analyses tiles fixed by a given size and an input step.
*
* This filter :
* _ creates output tiles with given size and step
@@ -105,6 +106,8 @@ public:
using Point2DType = itk::Point<double,2>;
using Point3DType = itk::Point<double,3>;
+ typedef itk::SimpleFastMutexLock MutexType;
+
// Typedef for iterators
typedef itk::ImageScanlineConstIterator< ImageInType > InputIterator;
typedef itk::ImageScanlineIterator< ImageOutType > OutputIterator;
@@ -123,7 +126,7 @@ protected:
TilesAnalysisWithStepImageFilter();
// Destructor
- ~TilesAnalysisWithStepImageFilter() = default;
+ ~TilesAnalysisWithStepImageFilter();
// Print
void PrintSelf(std::ostream & os, itk::Indent indent) const override;
@@ -134,7 +137,41 @@ protected:
* GenerateOutputInformation() in order to inform the pipeline execution model.
*/
void GenerateOutputInformation() override;
+
+
+ /** TilesAnalysisWithStepImageFilter needs a input requested region that corresponds to the margin and shifts
+ * into the requested region of our output requested region. The output requested region needs to be modified
+ * to construct as wanted tiles with input size and step.
+ * As such, TilesAnalysesImageFilter needs to provide an implementation for
+ * GenerateInputRequestedRegion() in order to inform the pipeline execution model.
+ * \sa ProcessObject::GenerateInputRequestedRegion() */
+ void GenerateInputRequestedRegion() override;
+ using Superclass::ThreadedGenerateData;
+ /** startIndex and stopIndex represent the indices of tiles to process in thread threadId and index for the
+ first tile */
+ void ThreadedGenerateData(unsigned int startTileId, unsigned int stopTileId,
+ unsigned int firstTileIndexL, unsigned int firstTileIndexC,
+ itk::ThreadIdType threadId);
+
+ void BeforeThreadedGenerateData();
+
+ /** 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(). */
+ static ITK_THREAD_RETURN_TYPE ThreaderCallback(void *arg);
+
+ /** Internal structure used for passing image data into the threading library */
+ struct ThreadStruct
+ {
+ Pointer Filter;
+ };
+
+ /**
+ * \sa ImageToImageFilter::ThreadedGenerateData(),
+ * ImageToImageFilter::GenerateData() */
+ void GenerateData() override;
+
private:
TilesAnalysisWithStepImageFilter(const Self&); // purposely not implemented
@@ -144,6 +181,8 @@ protected:
// Step to scan input image
unsigned int m_Step;
+ MutexType * m_Mutex;
+
};
} // End namespace otb
diff --git a/include/otbTilesAnalysisWithStepImageFilter.txx b/include/otbTilesAnalysisWithStepImageFilter.txx
index 3f3f62a..4a2d867 100644
--- a/include/otbTilesAnalysisWithStepImageFilter.txx
+++ b/include/otbTilesAnalysisWithStepImageFilter.txx
@@ -44,8 +44,20 @@ namespace otb
{
this->m_SizeTiles = 64;
this->m_Margin = 0;
+
+ m_Mutex = new MutexType();
+ }
+
+ /**
+ * Destructor with default initialization
+ */
+ template <class TImageIn, class TImageOut, class TFunction>
+ TilesAnalysisWithStepImageFilter< TImageIn, TImageOut, TFunction >::
+ ~TilesAnalysisWithStepImageFilter()
+ {
+ delete m_Mutex;
+ m_Mutex = 0;
}
-
/**
@@ -92,6 +104,269 @@ namespace otb
}
+
+ /**
+ * Method GenerateInputRequestedRegion
+ */
+ template<class TImageIn, class TImageOut, class TFunction>
+ void
+ TilesAnalysisWithStepImageFilter< TImageIn, TImageOut, TFunction >
+ ::GenerateInputRequestedRegion()
+ {
+ ////// Set Output Requestion Region (Contruct Tiles with given size and step for output) ///////
+ ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
+
+ // New Index : closest inf multiple of SizeTiles
+ ImageOutIndexType outputIndex;
+ outputIndex[0] = (static_cast<int>(outputRequestedRegion.GetIndex()[0]/this->m_SizeTiles))*
+ this->m_SizeTiles;
+ outputIndex[1] = (static_cast<int>(outputRequestedRegion.GetIndex()[1]/this->m_SizeTiles))*
+ this->m_SizeTiles;
+
+ // New Size : closest sup multiple of SizeTiles with addtional when difference between new and old index is
+ // too important
+ ImageOutSizeType outputSize;
+ outputSize[0] = (static_cast<int>((outputRequestedRegion.GetIndex()[0] - outputIndex[0] +
+ outputRequestedRegion.GetSize()[0])/this->m_SizeTiles) +
+ this->m_SizeTiles/m_Step)*this->m_SizeTiles;
+ outputSize[1] = (static_cast<int>((outputRequestedRegion.GetIndex()[1] - outputIndex[1] +
+ outputRequestedRegion.GetSize()[1])/this->m_SizeTiles) +
+ this->m_SizeTiles/m_Step)*this->m_SizeTiles;
+
+
+ std::cout << "outputSize = " << outputSize << std::endl;
+ std::cout << "outputIndex = " << outputIndex << std::endl;
+
+ // Region affectation (with Crop to respect image dimension)
+ outputRequestedRegion.SetSize(outputSize);
+ outputRequestedRegion.SetIndex(outputIndex);
+ outputRequestedRegion.Crop(this->GetOutput()->GetLargestPossibleRegion());
+ this->GetOutput()->SetRequestedRegion(outputRequestedRegion);
+
+ std::cout << "outputRequestedRegion Size : " << outputRequestedRegion << std::endl;
+
+ ///////////// With the new output requested region, find the region into input image /////////////
+ ImageInRegionType inputRequestedRegion = outputRequestedRegion;
+ ImageInPointer inputPtr = const_cast< ImageInType * >( this->GetImageInput() );
+ inputPtr->SetRequestedRegion(inputRequestedRegion);
+
+ }
+
+
+ /**
+ * Method GenerateData
+ */
+ template<class TImageIn, class TImageOut, class TFunction>
+ void
+ TilesAnalysisWithStepImageFilter< TImageIn, TImageOut, TFunction >
+ ::GenerateData()
+ {
+ // Allocate outputs
+ this->AllocateOutputs();
+
+ int nbThreads = this->GetNumberOfThreads();
+
+ // Get Output Requested region
+ 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;
+
+ // Check if additionnal tiles are required
+ if (outputRegion.GetSize()[1] % m_Step != 0)
+ {
+ ++nbTiles_Lines;
+ }
+ if (outputRegion.GetSize()[0] % m_Step != 0)
+ {
+ ++nbTiles_Col;
+ }
+
+ // First Tile = index of output requested region
+ int IndL_FirstTiles = outputRegion.GetIndex()[1];
+ int IndC_FirstTiles = outputRegion.GetIndex()[0];
+
+ ///////////////////////// Pre Processing ///////////////////////
+ // Call Function initialization
+ this->m_FunctionOnTiles->Initialize(nbThreads);
+
+ this->m_NbTiles = nbTiles_Lines*nbTiles_Col;
+ this->m_NbTilesCol = nbTiles_Col;
+ this->m_IndLFirstTile = IndL_FirstTiles;
+ this->m_IndCFirstTile = IndC_FirstTiles;
+
+ std::cout << "Into GenerateData : " << std::endl;
+ std::cout << this->m_NbTiles << std::endl;
+ std::cout << this->m_NbTilesCol << std::endl;
+ std::cout << this->m_IndLFirstTile << std::endl;
+ std::cout << this->m_IndCFirstTile << std::endl;
+
+ this->BeforeThreadedGenerateData();
+
+ ///////////////////////// Processing ///////////////////////
+ // Set up the multithreaded processing
+ ThreadStruct str;
+ str.Filter = this;
+
+ // Setting up multithreader
+ this->GetMultiThreader()->SetNumberOfThreads(this->GetNumberOfThreads());
+ this->GetMultiThreader()->SetSingleMethod(this->ThreaderCallback, &str);
+
+ // multithread the execution
+ this->GetMultiThreader()->SingleMethodExecute();
+
+ ///////////// Post Processing ////////////
+ // 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())
+ {
+ // Accumulations
+ OutIt.Set(0);
+ ++OutIt;
+ }
+ }
+
+
+ template<class TImageIn, class TImageOut, class TFunction>
+ void
+ TilesAnalysisWithStepImageFilter< TImageIn, TImageOut, TFunction >
+ ::ThreadedGenerateData(unsigned int startTileId, unsigned int stopTileId,
+ unsigned int firstTileIndexL, unsigned int firstTileIndexC,
+ itk::ThreadIdType threadId)
+ {
+ ImageInPointer inputPtr = const_cast< ImageInType * >( this->GetImageInput() );
+
+ // Typedef for Iterators on selected region
+ typedef itk::ImageRegionIterator<ImageOutType> OutItType;
+
+ for (unsigned int k = startTileId; k < stopTileId; k++) //<= ??????
+ {
+ // Transform Tiles id into index
+ ImageOutIndexType index_current;
+ int id_L = k/this->m_NbTilesCol; // quotient of Euclidean division (integer division)
+ int id_C = k - id_L*this->m_NbTilesCol; // remainder of Euclidean division
+ index_current[1] = firstTileIndexL + id_L*m_Step;
+ index_current[0] = firstTileIndexC + id_C*m_Step;
+
+ // Construct the current tile
+ ImageOutIndexType currentOutTilesIndex;
+ ImageOutSizeType currentOutTilesSize;
+
+ // Create regions with index and the given size
+ // Output tiles
+ currentOutTilesIndex[0] = index_current[0];
+ currentOutTilesIndex[1] = index_current[1];
+ currentOutTilesSize[0] = this->m_SizeTiles;
+ currentOutTilesSize[1] = this->m_SizeTiles;
+ ImageOutRegionType outputCurrentRegion;
+ outputCurrentRegion.SetIndex(currentOutTilesIndex);
+ outputCurrentRegion.SetSize(currentOutTilesSize);
+
+ outputCurrentRegion.Crop(this->GetOutput()->GetLargestPossibleRegion());
+
+ // Process the current tile
+ ImageOutPixelType * TilesCache = new ImageOutPixelType[outputCurrentRegion.GetSize()[1]*
+ outputCurrentRegion.GetSize()[0]];
+
+ // Call Function Process :
+ // Inputs : Ptr Image In and the output tile
+ // Output : Ptr on result array
+ this->m_FunctionOnTiles->Process(inputPtr, outputCurrentRegion, TilesCache, threadId);
+
+ int compteur = 0;
+
+ // Assignate Output with function result
+ // Iterators on output
+ OutItType OutIt(this->GetOutput(), outputCurrentRegion);
+ OutIt.GoToBegin();
+ m_Mutex->Lock();
+ while (!OutIt.IsAtEnd())
+ {
+ // Accumulations
+ ImageOutPixelType outTmp = OutIt.Get();
+ OutIt.Set(outTmp + TilesCache[compteur]);
+ ++compteur;
+ ++OutIt;
+ }
+ m_Mutex->Unlock();
+
+ delete [] TilesCache;
+ TilesCache = 0;
+
+ }
+ }
+
+
+ // Callback routine used by the threading library. This routine just calls
+ // the ThreadedGenerateData method after setting the correct region for this
+ // thread.
+ template<class TImageIn, class TImageOut, class TFunction>
+ ITK_THREAD_RETURN_TYPE
+ TilesAnalysisWithStepImageFilter< TImageIn, TImageOut, TFunction >
+ ::ThreaderCallback(void *arg)
+ {
+ ThreadStruct *str;
+ int threadId, threadCount;
+ unsigned int total, start, stop, indLFirstTile, indCFirstTile;
+
+ threadId = ((itk::MultiThreader::ThreadInfoStruct *) (arg))->ThreadID;
+ threadCount = ((itk::MultiThreader::ThreadInfoStruct *) (arg))->NumberOfThreads;
+ str = (ThreadStruct *) (((itk::MultiThreader::ThreadInfoStruct *) (arg))->UserData);
+
+ total = str->Filter->GetNbTiles();
+ indLFirstTile = str->Filter->GetIndLFirstTile();
+ indCFirstTile = str->Filter->GetIndCFirstTile();
+
+
+ std::cout << "Into TheaderCallback : " << std::endl;
+ std::cout << total << std::endl;
+ std::cout << indLFirstTile << std::endl;
+ std::cout << indCFirstTile << std::endl;
+
+ if (threadId < static_cast<int>(total))
+ {
+ start =
+ static_cast<unsigned int>(threadId *
+ std::ceil(static_cast<double>(total)/(static_cast<double>(threadCount))));
+ stop =
+ 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
+ // case the vertex at that index will be processed in the next strip.
+ if (start != stop)
+ {
+ str->Filter->ThreadedGenerateData(start, stop, indLFirstTile, indCFirstTile, threadId);
+ }
+ }
+ // else
+ // {
+ // otherwise don't use this thread. Sometimes the threads don't
+ // break up very well and it is just as efficient to leave a
+ // few threads idle.
+ // }
+
+ return ITK_THREAD_RETURN_VALUE;
+ }
--
GitLab