Skip to content
Snippets Groups Projects
Commit a373d872 authored by Gaëlle USSEGLIO's avatar Gaëlle USSEGLIO
Browse files

ENH : Override MultiThreading into SARGroupedByOrthoImageFilter to process...

ENH : Override MultiThreading into SARGroupedByOrthoImageFilter to process input region (not output) in parallel
parent 8311bbce
No related branches found
No related tags found
1 merge request!28Interfero dem
......@@ -76,7 +76,7 @@ namespace otb
AddParameter(ParameterType_InputImage, "topographicphase", "Input Topographic Phase (estimation with DEM projection)");
SetParameterDescription("topographicphase", "Input Topographic Phase (estimation with DEM projection).");
MandatoryOff("topographicphase");
AddParameter(ParameterType_InputImage, "incartmeanmaster", "Input Cartesian Mean Master image");
SetParameterDescription("incartmeanmaster", "Input Cartesian Mean Master image.");
......@@ -150,14 +150,11 @@ namespace otb
FloatVectorImageType::Pointer CartMeanPtr;
CartMeanPtr = GetParameterFloatVectorImage("incartmeanmaster");
filterGroupedBy->SetCartesianMeanInput(CartMeanPtr);
filterGroupedBy->SetCompensatedComplexInput(filterCompensatedComplex->GetOutput());
// Main Output
SetParameterOutputImage("out", filterGroupedBy->GetOutput());
}
// Vector for filters
std::vector<itk::ProcessObject::Pointer> m_Ref;
......
......@@ -24,6 +24,8 @@
#include "itkImageToImageFilter.h"
#include "itkSmartPointer.h"
#include "itkPoint.h"
#include "itkSimpleFastMutexLock.h"
#include "itkImageRegionSplitter.h"
#include "itkImageScanlineConstIterator.h"
#include "itkImageScanlineIterator.h"
......@@ -128,6 +130,10 @@ public:
using Point3DType = itk::Point<double,3>;
typedef double ValueType;
typedef itk::SimpleFastMutexLock MutexType;
typedef typename itk::ImageRegionSplitter<2> SplitterType; // Dimension = 2 for inputs
typedef typename SplitterType::Pointer SplitterPointer;
// Typedef for iterators
typedef itk::ImageScanlineConstIterator< ImageVectorType > InputIterator;
typedef itk::ImageScanlineIterator< ImageOutType > OutputIterator;
......@@ -150,6 +156,18 @@ public:
void SetSARImageKeyWorList(ImageKeywordlist sarImageKWL);
void SetDEMImagePtr(ImageDEMPointer demPtr);
// Getter for arrays called inot callback function
void GetArrays(double * &tmpArray_real, double * &tmpArray_imag, double * &tmpArray_mod,
double * &tmpArray_count)
{
// Assign Ptrs (=> needs reference)
tmpArray_real = m_tmpArray_real;
tmpArray_imag = m_tmpArray_imag;
tmpArray_mod = m_tmpArray_mod;
tmpArray_count = m_tmpArray_count;
}
protected:
// Constructor
SARGroupedByOrthoImageFilter();
......@@ -187,14 +205,33 @@ protected:
* Therefore, this implementation provides a ThreadedGenerateData() routine
* which is called for each processing thread. The main output image data is
* allocated automatically by the superclass prior to calling
* ThreadedGenerateData(). ThreadedGenerateData can only write to the
* portion of the output image specified by the parameter
* "outputRegionForThread"
*
* ThreadedGenerateData(). The implementation of multithreading is overrided here with the used of
* MultiThreader directly */
using Superclass::ThreadedGenerateData;
/** inputRegion and outputRegion represent the input/output region for the current streaming.
Temporary arrays are shared between all threads. */
void ThreadedGenerateData(const ImageVectorRegionType& inputRegion,
const ImageOutRegionType& outputRegion,
double * tmpArray_real, double * tmpArray_imag,
double * tmpArray_mod, double * tmpArray_count,
itk::ThreadIdType threadId);
/** 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() */
virtual void ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread,
itk::ThreadIdType threadId) ITK_OVERRIDE;
void GenerateData() override;
private:
......@@ -217,6 +254,13 @@ protected:
int m_nbLinesSAR;
int m_nbColSAR;
MutexType * m_Mutex;
// Temporary arrays
double * m_tmpArray_real;
double * m_tmpArray_imag;
double * m_tmpArray_mod;
double * m_tmpArray_count;
};
} // End namespace otb
......
......@@ -77,7 +77,8 @@ namespace otb
m_geoidEmg96 = new ossimGeoidEgm96(ossimFilename(ConfigurationManager::GetGeoidFile()));
}
}
m_Mutex = new MutexType();
}
/**
......@@ -91,6 +92,10 @@ namespace otb
delete m_geoidEmg96;
m_geoidEmg96 = 0;
}
delete m_Mutex;
m_Mutex = 0;
}
/**
......@@ -196,10 +201,11 @@ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
/////////////////////// Main Output : interferogram (into Ground geometry) ///////////////////////
// Vector Image :
// At Least 3 Components :
// At Least 4 Components :
// _ amplitude
// _ phase
// _ coherency
// _ count
outputPtr->SetNumberOfComponentsPerPixel(4);
// The output is defined with the DEM Image
......@@ -260,12 +266,9 @@ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
int nbColSAR = m_nbColSAR;
int nbLinesSAR = m_nbLinesSAR;
//std::cout << "nbColSAR and nbLineSAR : " << nbColSAR << " , " << nbLinesSAR << std::endl;
// Margin to apply on direct localisation for the four sides.
int margin = 100;
int margin = 1000;
// Indice for the SAR bloc (determined by the Pipeline)
ImageOutIndexType id[4] ;
......@@ -280,7 +283,6 @@ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
Point2DType pixelSAR_Into_MNT_LatLon(0);
Point2DType pixelSAR(0);
// ImageType::IndexType pixelSAR_Into_MNT;
itk::ContinuousIndex<double,2> pixelSAR_Into_MNT;
Point3DType demGeoPoint;
......@@ -297,9 +299,6 @@ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
firstC = nbColSAR;
lastC = 0;
intoInputImage = true;
//std::cout << "At the beggindg of OuputToInput intoInputImage = " << intoInputImage << " And firstC, firstL, lastC, lastL : " << firstC << ", " << firstL << ", " << lastC << ", " << lastL
// << std::endl;
// For each side of the output region
......@@ -318,8 +317,8 @@ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
gptPt.lat = demLatLonPoint[1];
h = m_geoidEmg96->offsetFromEllipsoid(gptPt);
}
//demGeoPoint[2] = h;
demGeoPoint[2] = 0;
demGeoPoint[2] = h;
//demGeoPoint[2] = 0;
// Call the method of sarSensorModelAdapter
m_SarSensorModelAdapter->WorldToLineSampleYZ(demGeoPoint, col_row, y_z);
......@@ -397,15 +396,8 @@ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
inputRequestedRegion.SetIndex(inputRequestedRegionIndex);
inputRequestedRegion.SetSize(inputRequestedRegionSize);
//std::cout << "Pouet : " << inputRequestedRegion << std::endl;
}
// std::cout << "At the end of OuputToInput : " << inputRequestedRegion << " with intoInputImage = " << intoInputImage << " And firstC, firstL, lastC, lastL : " << firstC << ", " << firstL << ", " << lastC << ", " << lastL
// << std::endl;
//std::cout << "At the end of OuputToInput : " << inputRequestedRegion.GetSize() << std::endl;
return inputRequestedRegion;
}
......@@ -430,8 +422,6 @@ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
ImageVectorPointer cartesianMeanPtr = const_cast< ImageVectorType * >( this->GetCartesianMeanInput() );
ImageVectorPointer compensatedComplexPtr = const_cast< ImageVectorType * >( this->GetCompensatedComplexInput() );
//std::cout << "inputRequestedRegion : " << inputRequestedRegion << std::endl;
//if (isIntoInputImages)
// {
cartesianMeanPtr->SetRequestedRegion(inputRequestedRegion);
......@@ -441,224 +431,259 @@ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
// {
// itkExceptionMacro(<<"DEM and SAR geometry do not match.");
// }
}
/**
* Method ThreadedGenerateDataWithoutTopographicPhase
/**
* Method GenerateData
*/
template<class TImageVector, class TImageDEM, class TImageOut>
void
SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
::ThreadedGenerateData(const ImageOutRegionType & outputRegionForThread,
itk::ThreadIdType threadId)
::GenerateData()
{
// Allocate outputs
this->AllocateOutputs();
int nbThreads = this->GetNumberOfThreads();
// Get Output Requested region
ImageOutRegionType outputRegion = this->GetOutput()->GetRequestedRegion();
// Compute corresponding input region
bool isIntoInputImage = true;
ImageVectorRegionType inputRegionForThread = OutputRegionToInputRegion(outputRegionForThread, isIntoInputImage);
ImageVectorRegionType inputRegion = OutputRegionToInputRegion(outputRegion, isIntoInputImage);
//std::cout << "Pouet : " << inputRegionForThread << " With " << isIntoInputImage << std::endl;
if (isIntoInputImage)
{
// Allocation for temporay arrays
int regionSize = outputRegion.GetSize()[0] * outputRegion.GetSize()[1];
m_tmpArray_real = new double[regionSize];
m_tmpArray_imag = new double[regionSize];
m_tmpArray_mod = new double[regionSize];
m_tmpArray_count = new double[regionSize];
for (unsigned int i = 0; i < regionSize; i++)
{
m_tmpArray_real[i] = static_cast<double>(0);
m_tmpArray_imag[i] = static_cast<double>(0);
m_tmpArray_mod[i] = static_cast<double>(0);
m_tmpArray_count[i] = static_cast<double>(0);
}
// Support progress methods/callbacks
itk::ProgressReporter progress( this, threadId, outputRegionForThread.GetNumberOfPixels() );
// Iterator on output (Ground geometry)
OutputIterator OutIt(this->GetOutput(), outputRegionForThread);
OutIt.GoToBegin();
// Temporary arrays to store outputRegion pixels
int regionSize = outputRegionForThread.GetSize()[0] * outputRegionForThread.GetSize()[1];
ValueType * tmpArray_real = new ValueType[regionSize];
ValueType * tmpArray_imag = new ValueType[regionSize];
ValueType * tmpArray_mod = new ValueType[regionSize];
ValueType * tmpArray_count = new ValueType[regionSize];
// Set up the multithreaded processing
ThreadStruct str;
str.Filter = this;
Point2DType demLonLatPoint;
itk::ContinuousIndex<double,2> pixel_Into_DEM_index;
ImageOutIndexType index;
// Setting up multithreader
this->GetMultiThreader()->SetNumberOfThreads(this->GetNumberOfThreads());
this->GetMultiThreader()->SetSingleMethod(this->ThreaderCallback, &str);
// Init Temporary array (to zero)
for (unsigned int i = 0; i < regionSize; i++)
{
tmpArray_real[i] = static_cast<ValueType>(0);
tmpArray_imag[i] = static_cast<ValueType>(0);
tmpArray_mod[i] = static_cast<ValueType>(0);
tmpArray_count[i] = static_cast<ValueType>(0);
}
// multithread the execution
this->GetMultiThreader()->SingleMethodExecute();
// Loop on output
ImageOutPixelType pixelOutput;
pixelOutput.Reserve(4);
int index_intoArray = 0;
// Iterator on output (Ground geometry)
OutputIterator OutIt(this->GetOutput(), outputRegion);
OutIt.GoToBegin();
// Scan input and estimate ortho geometry if output Region corresponds to input SAR image
if (isIntoInputImage)
{
// Iterators on inputs (SAR geometry)
InputIterator CartMeanIt(this->GetCartesianMeanInput(), inputRegionForThread);
InputIterator CompComplexIt(this->GetCompensatedComplexInput(), inputRegionForThread);
CartMeanIt.GoToBegin();
CompComplexIt.GoToBegin();
// Scan the input Region
// for each line
while ( !CartMeanIt.IsAtEnd() && !CompComplexIt.IsAtEnd())
// For each line
while ( !OutIt.IsAtEnd())
{
CartMeanIt.GoToBeginOfLine();
CompComplexIt.GoToBeginOfLine();
// For each column
while (!CartMeanIt.IsAtEndOfLine() && !CompComplexIt.IsAtEndOfLine())
OutIt.GoToBeginOfLine();
// For each colunm
while (!OutIt.IsAtEndOfLine())
{
// Get elt from cartesian mean image (X, Y and Z) and create an ECEF point with its
ossimEcefPoint cartPoint(CartMeanIt.Get()[0], CartMeanIt.Get()[1], CartMeanIt.Get()[2]);
if (m_tmpArray_count[index_intoArray] > 0)
{
double mod_Acc = sqrt(m_tmpArray_real[index_intoArray]*
m_tmpArray_real[index_intoArray] +
m_tmpArray_imag[index_intoArray]*
m_tmpArray_imag[index_intoArray]);
// Amplitude
pixelOutput[0] = m_Gain *
sqrt((m_tmpArray_mod[index_intoArray]/(m_tmpArray_count[index_intoArray])));
// Conversion into World point
ossimGpt worldPoint(cartPoint);
// Phase
pixelOutput[1] = std::atan2(m_tmpArray_imag[index_intoArray],
m_tmpArray_real[index_intoArray]);
demLonLatPoint[0] = worldPoint.lon;
demLonLatPoint[1] = worldPoint.lat;
// Get Nearest index into DEM/Ground Geometry
// Transform into continuous index
m_DEMPtr->TransformPhysicalPointToContinuousIndex(demLonLatPoint,
pixel_Into_DEM_index);
// Integer index
index [0] = static_cast<int>(pixel_Into_DEM_index[0]);
index [1] = static_cast<int>(pixel_Into_DEM_index[1]);
// Check if index into outputRegionForThread
if (index[0] >= outputRegionForThread.GetIndex()[0] &&
index[0] < (outputRegionForThread.GetIndex()[0] +
outputRegionForThread.GetSize()[0]) &&
index[1] >= outputRegionForThread.GetIndex()[1] &&
index[1] < (outputRegionForThread.GetIndex()[1] +
outputRegionForThread.GetSize()[1]))
// Mod 2*Pi
pixelOutput[1] = pixelOutput[1]-(2*M_PI)*floor(pixelOutput[1]/(2*M_PI));
// Coherency
pixelOutput[2] = mod_Acc / m_tmpArray_mod[index_intoArray] ;
// IsData
pixelOutput[3] = m_tmpArray_count[index_intoArray];
}
else
{
//std::cout << "index : " << index << std::endl;
// Accumulations for CompensatedComplex for current ground index
int index_intoArray = static_cast<int>((index[0]-outputRegionForThread.GetIndex()[0]) +
(index[1] - outputRegionForThread.GetIndex()[1]) *
outputRegionForThread.GetSize()[0]);
if (index_intoArray > regionSize)
{
std::cout << "WTF !! " << std::endl;
std::cout << "index[0] = " << index[0] << " And " << "index[1] = " << index[1] <<
std::endl;
std::cout << "outputRegionForThread.GetIndex()[0] = " <<
outputRegionForThread.GetIndex()[0] << " outputRegionForThread.GetIndex()[1] = " <<
outputRegionForThread.GetIndex()[1] << std::endl;
std::cout << "outputRegionForThread.GetSize()[0] = " <<
outputRegionForThread.GetSize()[0] << " outputRegionForThread.GetSize()[1] = " <<
outputRegionForThread.GetSize()[1] << std::endl;
std::cout << "index_intoArray = " << index_intoArray << " And index_intoArray " <<
regionSize << std::endl;
}
if (index[0] == 1261 && index[1] == 634)
{
std::cout << "Rien ????????????" << std::endl;
}
if (index[0] == 1557 && index[1] == 2038)
{
std::cout << "de la donnee !!!!!!" << std::endl;
}
tmpArray_real[index_intoArray] += CompComplexIt.Get()[0];
tmpArray_imag[index_intoArray] += CompComplexIt.Get()[1];
tmpArray_mod[index_intoArray] += CompComplexIt.Get()[2];
++tmpArray_count[index_intoArray];
pixelOutput[0] = 0;
pixelOutput[1] = 0;
pixelOutput[2] = 0;
pixelOutput[3] = 0;
}
++ CartMeanIt;
++ CompComplexIt;
// Assigne Main output (Groun geometry)
OutIt.Set(pixelOutput);
// Next colunm
++OutIt;
++index_intoArray;
}
// Next line
CartMeanIt.NextLine();
CompComplexIt.NextLine();
OutIt.NextLine();
}
} // End scan inputs
// Free Memory
delete m_tmpArray_real;
m_tmpArray_real = 0;
delete m_tmpArray_imag;
m_tmpArray_imag = 0;
delete m_tmpArray_mod;
m_tmpArray_mod = 0;
delete m_tmpArray_count;
m_tmpArray_count = 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 TImageVector, class TImageDEM, class TImageOut>
ITK_THREAD_RETURN_TYPE
SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
::ThreaderCallback(void *arg)
{
ThreadStruct *str;
int threadId, threadCount;
threadId = ((itk::MultiThreader::ThreadInfoStruct *) (arg))->ThreadID;
threadCount = ((itk::MultiThreader::ThreadInfoStruct *) (arg))->NumberOfThreads;
str = (ThreadStruct *) (((itk::MultiThreader::ThreadInfoStruct *) (arg))->UserData);
// Get filter elts (temporary arrays and vector of input regions)
double * tmpArray_real;
double * tmpArray_imag;
double * tmpArray_mod;
double * tmpArray_count;
str->Filter->GetArrays(tmpArray_real, tmpArray_imag, tmpArray_mod, tmpArray_count);
// Get inputRegion for current thread (Split the input region to get smaller ones : one per thread)
SplitterPointer splitter = SplitterType::New();
ImageVectorRegionType inputRegionForThread = splitter->GetSplit(threadId, threadCount,
str->Filter->GetInput()->GetRequestedRegion());
// Loop on output
ImageOutPixelType pixelOutput;
pixelOutput.Reserve(4);
int index_intoArray = 0;
// For each line
while ( !OutIt.IsAtEnd())
str->Filter->ThreadedGenerateData(inputRegionForThread,
str->Filter->GetOutput()->GetRequestedRegion(),
tmpArray_real, tmpArray_imag,
tmpArray_mod, tmpArray_count,
threadId);
return ITK_THREAD_RETURN_VALUE;
}
/**
* Method ThreadedGenerateData
*/
template<class TImageVector, class TImageDEM, class TImageOut>
void
SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
::ThreadedGenerateData(const ImageVectorRegionType& inputRegion,
const ImageOutRegionType& outputRegion,
double * tmpArray_real, double * tmpArray_imag,
double * tmpArray_mod, double * tmpArray_count,
itk::ThreadIdType threadId)
{
Point2DType demLonLatPoint;
itk::ContinuousIndex<double,2> pixel_Into_DEM_index;
ImageOutIndexType index;
int regionSize = outputRegion.GetSize()[0] * outputRegion.GetSize()[1];
// Iterators on inputs (SAR geometry)
InputIterator CartMeanIt(this->GetCartesianMeanInput(), inputRegion);
InputIterator CompComplexIt(this->GetCompensatedComplexInput(), inputRegion);
CartMeanIt.GoToBegin();
CompComplexIt.GoToBegin();
// Scan the input Region
// for each line
while ( !CartMeanIt.IsAtEnd() && !CompComplexIt.IsAtEnd())
{
OutIt.GoToBeginOfLine();
CartMeanIt.GoToBeginOfLine();
CompComplexIt.GoToBeginOfLine();
// For each colunm
while (!OutIt.IsAtEndOfLine())
// For each column
while (!CartMeanIt.IsAtEndOfLine() && !CompComplexIt.IsAtEndOfLine())
{
if (tmpArray_count[index_intoArray] > 0)
{
double mod_Acc = sqrt(tmpArray_real[index_intoArray]*tmpArray_real[index_intoArray] +
tmpArray_imag[index_intoArray]*tmpArray_imag[index_intoArray]);
// Amplitude
pixelOutput[0] = m_Gain *
sqrt((tmpArray_mod[index_intoArray]/(tmpArray_count[index_intoArray])));
// Get elt from cartesian mean image (X, Y and Z) and create an ECEF point with its
ossimEcefPoint cartPoint(CartMeanIt.Get()[0], CartMeanIt.Get()[1], CartMeanIt.Get()[2]);
// Phase
pixelOutput[1] = std::atan2(tmpArray_imag[index_intoArray],
tmpArray_real[index_intoArray]);
// Conversion into World point
ossimGpt worldPoint(cartPoint);
// Mod 2*Pi
pixelOutput[1] = pixelOutput[1]-(2*M_PI)*floor(pixelOutput[1]/(2*M_PI));
// Coherency
pixelOutput[2] = mod_Acc / tmpArray_mod[index_intoArray] ;
// IsData
pixelOutput[3] = tmpArray_count[index_intoArray];
}
else
demLonLatPoint[0] = worldPoint.lon;
demLonLatPoint[1] = worldPoint.lat;
// Get Nearest index into DEM/Ground Geometry
// Transform into continuous index
m_DEMPtr->TransformPhysicalPointToContinuousIndex(demLonLatPoint,
pixel_Into_DEM_index);
// Integer index
index [0] = static_cast<int>(pixel_Into_DEM_index[0]);
index [1] = static_cast<int>(pixel_Into_DEM_index[1]);
// Check if index into outputRegion
if (index[0] >= outputRegion.GetIndex()[0] &&
index[0] < (outputRegion.GetIndex()[0] +
outputRegion.GetSize()[0]) &&
index[1] >= outputRegion.GetIndex()[1] &&
index[1] < (outputRegion.GetIndex()[1] +
outputRegion.GetSize()[1]))
{
pixelOutput[0] = 0;
pixelOutput[1] = 0;
pixelOutput[2] = 0;
pixelOutput[3] = 0;
// Accumulations for CompensatedComplex for current ground index
int index_intoArray = static_cast<int>((index[0]-outputRegion.GetIndex()[0]) +
(index[1] - outputRegion.GetIndex()[1]) *
outputRegion.GetSize()[0]);
// Mutex to protect class arguments (m_*Arrays)
m_Mutex->Lock();
tmpArray_real[index_intoArray] += CompComplexIt.Get()[0];
tmpArray_imag[index_intoArray] += CompComplexIt.Get()[1];
tmpArray_mod[index_intoArray] += CompComplexIt.Get()[2];
++tmpArray_count[index_intoArray];
m_Mutex->Unlock();
}
// Assigne Main output (Groun geometry)
OutIt.Set(pixelOutput);
progress.CompletedPixel();
// Next colunm
++OutIt;
++index_intoArray;
++ CartMeanIt;
++ CompComplexIt;
}
// Next line
OutIt.NextLine();
}
CartMeanIt.NextLine();
CompComplexIt.NextLine();
} // End scan inputs
// Free Memory
delete tmpArray_real;
tmpArray_real = 0;
delete tmpArray_imag;
tmpArray_imag = 0;
delete tmpArray_mod;
tmpArray_mod = 0;
delete tmpArray_count;
tmpArray_count = 0;
}
} /*namespace otb*/
} /*namespace otb*/
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment