diff --git a/app/otbSAROrthoInterferogram.cxx b/app/otbSAROrthoInterferogram.cxx
index 8943cf069370caf37f00679a9b08b54908ebaa31..9342fce752eb575edceb0f29df5f4a42df5e0ffe 100644
--- a/app/otbSAROrthoInterferogram.cxx
+++ b/app/otbSAROrthoInterferogram.cxx
@@ -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;
diff --git a/include/otbSARGroupedByOrthoImageFilter.h b/include/otbSARGroupedByOrthoImageFilter.h
index 37f647a1ad98b9b5a028f9aa7687dfe2f18a7c0f..bd2ae6a339ab7b8150b2416b16c9ac3043367fa4 100644
--- a/include/otbSARGroupedByOrthoImageFilter.h
+++ b/include/otbSARGroupedByOrthoImageFilter.h
@@ -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
diff --git a/include/otbSARGroupedByOrthoImageFilter.txx b/include/otbSARGroupedByOrthoImageFilter.txx
index bb6105306270dab41d481853fa55f659c82872fd..4dc7ba1ca258a7363118a692aba1726769103959 100644
--- a/include/otbSARGroupedByOrthoImageFilter.txx
+++ b/include/otbSARGroupedByOrthoImageFilter.txx
@@ -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