From bf33edd9ed5ab8e6546bf0f07d26ea904074a0eb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ga=C3=ABlle=20USSEGLIO?= <gaelle.usseglio@cnes.fr>
Date: Wed, 16 Oct 2019 14:33:32 +0000
Subject: [PATCH] ENH: Cleaning into SARInterferogram Filter/Application
---
app/otbSARInterferogram.cxx | 37 -
app/otbSARRobustInterferogram.cxx | 47 --
include/otbSARInterferogramImageFilter.h | 42 +-
include/otbSARInterferogramImageFilter.txx | 856 +++++++--------------
python_src/diapOTB.py | 6 -
python_src/diapOTB_S1IW.py | 16 +-
6 files changed, 264 insertions(+), 740 deletions(-)
diff --git a/app/otbSARInterferogram.cxx b/app/otbSARInterferogram.cxx
index 4a4d4d4..3bcc988 100644
--- a/app/otbSARInterferogram.cxx
+++ b/app/otbSARInterferogram.cxx
@@ -70,10 +70,6 @@ namespace otb
AddParameter(ParameterType_InputImage, "insarmaster", "Input SAR Master image");
SetParameterDescription("insarmaster", "Input SAR Master image.");
- AddParameter(ParameterType_InputImage, "indem", "Input DEM (for metadata only)");
- SetParameterDescription("indem", "Input DEM (for metadata only).");
- MandatoryOff("indem");
-
AddParameter(ParameterType_InputImage, "topographicphase", "Input Topographic Phase (estimation with DEM projection)");
SetParameterDescription("topographicphase", "Input Topographic Phase (estimation with DEM projection).");
MandatoryOff("topographicphase");
@@ -108,17 +104,9 @@ namespace otb
SetMinimumParameterFloatValue("gain", 0);
MandatoryOff("gain");
- AddParameter(ParameterType_Bool, "ortho", "Build interferogram into ortho geometry");
- SetParameterDescription("ortho", "If true, then build interferogram into ortho geometry.");
-
AddParameter(ParameterType_OutputImage, "out", "Interferogram");
SetParameterDescription("out","Output Vector Image : Interferogram.");
- AddParameter(ParameterType_OutputImage, "outopt", "Interferogram into ortho geometry");
- SetParameterDescription("outopt","Output Vector Image : Interferogram into ortho geometry.");
- MandatoryOff("outopt");
-
-
AddRAMParameter();
SetDocExampleParameterValue("insarslave","s1b-s4-slc-vv-20160929t014610-20160929t014634-002277-003d71-002.tiff");
@@ -151,16 +139,6 @@ namespace otb
otbAppLogINFO(<<"Averaging Margin on range : "<<margin_Ran);
otbAppLogINFO(<<"Gain Factor : "<<factor_gain);
- // Check if ortho is enabled
- if (IsParameterEnabled("ortho"))
- {
- if (!GetParameterByKey("indem")->HasValue())
- {
- otbAppLogFATAL(<<"ortho activate. indem must be assigned");
- }
- }
-
-
///////////////////////////////////// Interferogram Filter ///////////////////////////////////////////////
// Instanciate the Interferogram filter
InterferogramFilterType::Pointer filterInterferogram = InterferogramFilterType::New();
@@ -170,12 +148,6 @@ namespace otb
filterInterferogram->SetMarginRan(margin_Ran);
filterInterferogram->SetMarginAzi(margin_Azi);
filterInterferogram->SetGain(factor_gain);
- filterInterferogram->SetUseDEMGeoAsOutput(false);
-
- if (GetParameterByKey("indem")->HasValue())
- {
- filterInterferogram->SetDEMPtr(GetParameterFloatImage("indem"));
- }
// Execute the main Pipeline
ComplexFloatImageType::Pointer SARMasterPtr;
@@ -192,20 +164,11 @@ namespace otb
if (GetParameterByKey("topographicphase")->HasValue())
{
filterInterferogram->SetTopographicPhaseInput(GetParameterImage("topographicphase"));
-
- if (IsParameterEnabled("ortho"))
- {
- filterInterferogram->SetUseDEMGeoAsOutput(true);
- }
}
// Main Output
SetParameterOutputImage("out", filterInterferogram->GetOutput());
- if (IsParameterEnabled("ortho"))
- {
- SetParameterOutputImage("outopt", filterInterferogram->GetOutputIntoDEMGeo());
- }
}
// Vector for filters
std::vector<itk::ProcessObject::Pointer> m_Ref;
diff --git a/app/otbSARRobustInterferogram.cxx b/app/otbSARRobustInterferogram.cxx
index 01f93c8..781d7a3 100644
--- a/app/otbSARRobustInterferogram.cxx
+++ b/app/otbSARRobustInterferogram.cxx
@@ -75,10 +75,6 @@ namespace otb
AddParameter(ParameterType_InputImage, "insarmaster", "Input SAR Master image");
SetParameterDescription("insarmaster", "Input SAR Master image.");
- AddParameter(ParameterType_InputImage, "indem", "Input DEM (for metadata only)");
- SetParameterDescription("indem", "Input DEM (for metadata only).");
- MandatoryOff("indem");
-
AddParameter(ParameterType_InputImage, "ingrid", "Input deformation grid");
SetParameterDescription("ingrid", "Input deformation grid.");
@@ -122,17 +118,10 @@ namespace otb
SetDefaultParameterFloat("gain", 0.1);
SetMinimumParameterFloatValue("gain", 0);
MandatoryOff("gain");
-
- AddParameter(ParameterType_Bool, "ortho", "Build interferogram into ortho geometry");
- SetParameterDescription("ortho", "If true, then build interferogram into ortho geometry.");
AddParameter(ParameterType_OutputImage, "out", "Output interferogram (ML geometry)");
SetParameterDescription("out","Output interferogram (ML geometry).");
- AddParameter(ParameterType_OutputImage, "outopt", "Interferogram into ortho geometry");
- SetParameterDescription("outopt","Output Vector Image : Interferogram into ortho geometry.");
- MandatoryOff("outopt");
-
AddRAMParameter();
SetDocExampleParameterValue("insarslave","s1b-s4-slc-vv-20160929t014610-20160929t014634-002277-003d71-002.tiff");
@@ -168,21 +157,6 @@ namespace otb
otbAppLogINFO(<<"Averaging Margin on azimut :"<<margin_Azi);
otbAppLogINFO(<<"Averaging Margin on range : "<<margin_Ran);
otbAppLogINFO(<<"Gain : "<<gain);
-
- // Check if ortho is enabled
- if (IsParameterEnabled("ortho"))
- {
- if (!GetParameterByKey("outopt")->HasValue())
- {
- otbAppLogFATAL(<<"ortho activate. outopt must be assigned");
- }
-
- if (!GetParameterByKey("indem")->HasValue())
- {
- otbAppLogFATAL(<<"ortho activate. indem must be assigned");
- }
- }
-
// Retrive inputs
FloatVectorImageType::Pointer GridPtr = GetParameterFloatVectorImage("ingrid");
@@ -202,11 +176,6 @@ namespace otb
GetInternalApplication("TopographicPhaseApp")->SetParameterInt("gridstepazimut", gridStepAzi);
GetInternalApplication("TopographicPhaseApp")->SetParameterInt("mlran", mlRan);
GetInternalApplication("TopographicPhaseApp")->SetParameterInt("mlazi", mlAzi);
-
- if (IsParameterEnabled("ortho"))
- {
- GetInternalApplication("TopographicPhaseApp")->EnableParameter("cartcopy");
- }
ExecuteInternal("TopographicPhaseApp");
@@ -218,34 +187,18 @@ namespace otb
GetInternalApplication("InterferogramApp")->SetParameterInputImage("topographicphase",
GetInternalApplication("TopographicPhaseApp")->GetParameterOutputImage("out"));
- if (GetParameterByKey("indem")->HasValue())
- {
- FloatImageType::Pointer inputDEM = GetParameterFloatImage("indem");
- GetInternalApplication("InterferogramApp")->SetParameterInputImage("indem",
- inputDEM);
- }
GetInternalApplication("InterferogramApp")->SetParameterInt("mlran", mlRan);
GetInternalApplication("InterferogramApp")->SetParameterInt("mlazi", mlAzi);
GetInternalApplication("InterferogramApp")->SetParameterInt("marginran", margin_Ran);
GetInternalApplication("InterferogramApp")->SetParameterInt("marginazi", margin_Azi);
GetInternalApplication("InterferogramApp")->SetParameterFloat("gain", gain);
- if (IsParameterEnabled("ortho"))
- {
- GetInternalApplication("InterferogramApp")->EnableParameter("ortho");
- }
ExecuteInternal("InterferogramApp");
// Main output
SetParameterOutputImage("out", GetInternalApplication("InterferogramApp")->GetParameterOutputImage("out"));
- // Optionnal output
- if (IsParameterEnabled("ortho"))
- {
- SetParameterOutputImage("outopt", GetInternalApplication("InterferogramApp")->GetParameterOutputImage("outopt"));
- }
-
}
// Vector for filters
std::vector<itk::ProcessObject::Pointer> m_Ref;
diff --git a/include/otbSARInterferogramImageFilter.h b/include/otbSARInterferogramImageFilter.h
index a1b84cd..befd8cc 100644
--- a/include/otbSARInterferogramImageFilter.h
+++ b/include/otbSARInterferogramImageFilter.h
@@ -145,16 +145,13 @@ public:
itkSetMacro(MarginRan, unsigned int);
itkSetMacro(MarginAzi, unsigned int);
itkSetMacro(Gain, float);
- itkSetMacro(UseDEMGeoAsOutput, bool);
- itkSetMacro(ApproxDiapason, bool);
- void SetDEMPtr(ImageDEMPointer DEMPtr);
+
// Getter
itkGetMacro(MLRan, unsigned int);
itkGetMacro(MLAzi, unsigned int);
itkGetMacro(MarginRan, unsigned int);
itkGetMacro(MarginAzi, unsigned int);
itkGetMacro(Gain, float);
- itkGetMacro(UseDEMGeoAsOutput, bool);
// Setter/Getter for inputs (SAR images and potentially topographic phase)
/** Connect one of the operands for interferometry : Master */
@@ -171,13 +168,6 @@ public:
const ImageSARType * GetSlaveInput() const;
const ImagePhaseType * GetTopographicPhaseInput() const;
- // Getter for outputs (Interferogram into ML geometry and potentially into DEM geometry)
- /** Returns the const image */
- const TImageOut * GetOutput() const;
- /** Returns the image */
- TImageOut * GetOutput();
- /** Returns the optionnal output image (into DEM/Ortho geometry) */
- ImageOutType * GetOutputIntoDEMGeo();
protected:
// Constructor
@@ -219,16 +209,16 @@ protected:
* SARInterferogramImageFilter can produce an optionnal image if UseDEMGeoAsOutput = true. The requested
* region for this optionnal output is set to the largest possible region.
*/
- void EnlargeOutputRequestedRegion( itk::DataObject *output ) ITK_OVERRIDE;
+ //void EnlargeOutputRequestedRegion( itk::DataObject *output ) ITK_OVERRIDE;
/**
* SARInterferogramImageFilter can produce an optionnal image according to the functor. An allocation of
* this optionnal output is made into BeforeThreadedGenerateData and free memory into
* AfterThreadedGenerateData.
*/
- void BeforeThreadedGenerateData() ITK_OVERRIDE;
+ //void BeforeThreadedGenerateData() ITK_OVERRIDE;
- void AfterThreadedGenerateData() ITK_OVERRIDE;
+ //void AfterThreadedGenerateData() ITK_OVERRIDE;
/**
* SARInterferogramImageFilterr can be implemented as a multithreaded filter.
@@ -263,12 +253,6 @@ protected:
// gain for amplitude estimation
float m_Gain;
- // DEM Pointer (metadata only)
- ImageDEMPointer m_DEMPtr;
-
- // Use DEM Geometry as optionnal output
- bool m_UseDEMGeoAsOutput;
-
// SAR Dimensions
int m_nbLinesSAR;
int m_nbColSAR;
@@ -277,28 +261,10 @@ protected:
unsigned int m_MarginRan;
unsigned int m_MarginAzi;
- // Counter to synchronize our outputs with calculation
- int m_OutputCounter;
-
- // Mutex for optionnalImage (geo ortho)
- MutexType * m_Mutex;
-
- // Arrays to store ortho accumlations and counters
- double * m_OrthoRealAcc;
- double * m_OrthoImagAcc;
- double * m_OrthoModAcc;
- double * m_OrthoCounter;
- // One arrays by thread to guarantee Thread-Safety
- float ** m_OrthoRealAccByThread;
- float ** m_OrthoImagAccByThread;
- float ** m_OrthoModAccByThread;
- float ** m_OrthoCounterByThread;
-
// Dimensions of DEM
int m_nbLinesDEM;
int m_nbColDEM;
- bool m_ApproxDiapason;
};
} // End namespace otb
diff --git a/include/otbSARInterferogramImageFilter.txx b/include/otbSARInterferogramImageFilter.txx
index 9ae0734..4159ac5 100644
--- a/include/otbSARInterferogramImageFilter.txx
+++ b/include/otbSARInterferogramImageFilter.txx
@@ -35,9 +35,6 @@
#include <algorithm>
#include <omp.h>
-#ifndef M_1_PI
-#define M_1_PI 0.31830988618379067154
-#endif
namespace otb
{
@@ -46,23 +43,12 @@ namespace otb
*/
template <class TImageSAR, class TImageDEM, class TImagePhase, class TImageOut>
SARInterferogramImageFilter< TImageSAR, TImageDEM, TImagePhase, TImageOut >::SARInterferogramImageFilter()
- : m_MLRan(1), m_MLAzi(1), m_Gain(1), m_UseDEMGeoAsOutput(false), m_nbLinesSAR(0), m_nbColSAR(0),
- m_MarginRan(1), m_MarginAzi(1), m_nbLinesDEM(0), m_nbColDEM(0), m_ApproxDiapason(false)
+ : m_MLRan(1), m_MLAzi(1), m_Gain(1), m_nbLinesSAR(0), m_nbColSAR(0),
+ m_MarginRan(1), m_MarginAzi(1), m_nbLinesDEM(0), m_nbColDEM(0)
{
// Inputs required and/or needed
this->SetNumberOfRequiredInputs(2);
this->SetNumberOfIndexedInputs(3);
-
- // Outputs required and/or needed (one required and one optionnal)
- this->SetNumberOfRequiredOutputs(1);
- this->SetNumberOfIndexedOutputs(2);
-
- this->SetNthOutput(0, ImageOutType::New()); // ML Geo
- this->SetNthOutput(1, ImageOutType::New()); // DEM Geo
-
- m_OutputCounter = 0;
-
- m_Mutex = new MutexType();
}
/**
@@ -87,17 +73,6 @@ namespace otb
os << indent << "Gain for amplitude estimation : " << m_Gain << std::endl;
}
- /**
- * SetDEMPtr
- */
- template<class TImageSAR, class TImageDEM, class TImagePhase, class TImageOut>
- void
- SARInterferogramImageFilter< TImageSAR, TImageDEM, TImagePhase, TImageOut >
- ::SetDEMPtr(ImageDEMPointer demPtr)
- {
- m_DEMPtr = demPtr;
- }
-
/**
* Set Master Image
*/
@@ -192,47 +167,6 @@ namespace otb
return static_cast<const ImagePhaseType *>(this->itk::ProcessObject::GetInput(2));
}
-
- /**
- * Get Main Output : Interferogram into ML geometry
- */
- template<class TImageSAR, class TImageDEM, class TImagePhase, class TImageOut>
- TImageOut *
- SARInterferogramImageFilter< TImageSAR, TImageDEM, TImagePhase, TImageOut >
- ::GetOutput()
- {
- if (this->GetNumberOfOutputs() < 1)
- {
- return ITK_NULLPTR;
- }
- return static_cast<ImageOutType *>(this->itk::ProcessObject::GetOutput(0));
- }
-
- template<class TImageSAR, class TImageDEM, class TImagePhase, class TImageOut>
- const TImageOut *
- SARInterferogramImageFilter< TImageSAR, TImageDEM, TImagePhase, TImageOut >
- ::GetOutput() const
- {
- if (this->GetNumberOfOutputs() < 1)
- {
- return 0;
- }
- return static_cast<ImageOutType *>(this->itk::ProcessObject::GetOutput(0));
- }
-
- template<class TImageSAR, class TImageDEM, class TImagePhase, class TImageOut>
- typename SARInterferogramImageFilter< TImageSAR, TImageDEM, TImagePhase, TImageOut >::ImageOutType *
- SARInterferogramImageFilter< TImageSAR, TImageDEM, TImagePhase, TImageOut >
- ::GetOutputIntoDEMGeo()
- {
- if (this->GetNumberOfOutputs() < 2)
- {
- return ITK_NULLPTR;
- }
- return static_cast<ImageOutType *>(this->itk::ProcessObject::GetOutput(1));
- }
-
-
/**
* Method GenerateOutputInformaton()
**/
@@ -327,63 +261,6 @@ namespace otb
// itkExceptionMacro(<<"Margin for averaging can't be superior to ML Factors.");
// }
}
-
-
- //////////////////////////// Optionnal Output : Interferogram into DEM Geo ///////////////////////////
- if (m_UseDEMGeoAsOutput)
- {
- // Get Pointer
- ImageOutPointer outputOptionnalPtr = this->GetOutputIntoDEMGeo();
-
- // Vector Image :
- // At Least 4 Components :
- // _ Amplitude
- // _ Phase
- // _ Coherence
- // _ IsData
- outputOptionnalPtr->SetNumberOfComponentsPerPixel(4);
-
- // This output is defined with the DEM geometry
- outputOptionnalPtr->SetLargestPossibleRegion(m_DEMPtr->GetLargestPossibleRegion());
- outputOptionnalPtr->SetOrigin(m_DEMPtr->GetOrigin());
- outputOptionnalPtr->SetSignedSpacing(m_DEMPtr->GetSignedSpacing());
-
- // Projection Ref
- outputOptionnalPtr->SetProjectionRef(m_DEMPtr->GetProjectionRef());
-
- // Interferogram Ortho requires Topographic phase in input
- if (this->GetTopographicPhaseInput() == nullptr ||
- this->GetTopographicPhaseInput()->GetNumberOfComponentsPerPixel() < 3)
- {
- itkExceptionMacro(<<"Interferogram into ortho geometry requires Topographic phase in input with cartesian coordonates.");
- }
-
- if (m_OutputCounter == 0)
- {
- if (m_DEMPtr == nullptr)
- {
- itkExceptionMacro(<<"To estimate interferogram into ortho geometry, a DEM must be in input.");
- }
-
- // Get DEM's size
- m_nbLinesDEM = m_DEMPtr->GetLargestPossibleRegion().GetSize()[1];
- m_nbColDEM = m_DEMPtr->GetLargestPossibleRegion().GetSize()[0];
-
- // Allocation of arrays for accumulation and counter into ortho geo
- m_OrthoRealAcc = new double[m_nbLinesDEM * m_nbColDEM];
- m_OrthoImagAcc = new double[m_nbLinesDEM * m_nbColDEM];
- m_OrthoModAcc = new double[m_nbLinesDEM * m_nbColDEM];
- m_OrthoCounter = new double[m_nbLinesDEM * m_nbColDEM];
-
- // Initialize to zero
- memset(m_OrthoRealAcc, 0, m_nbLinesDEM * m_nbColDEM * sizeof(double));
- memset(m_OrthoImagAcc, 0, m_nbLinesDEM * m_nbColDEM * sizeof(double));
- memset(m_OrthoModAcc, 0, m_nbLinesDEM * m_nbColDEM * sizeof(double));
- memset(m_OrthoCounter, 0, m_nbLinesDEM * m_nbColDEM * sizeof(double));
- }
- }
-
- ++m_OutputCounter;
}
/**
@@ -571,104 +448,25 @@ namespace otb
// Get Output requested region
ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
- if (m_OutputCounter == 1)
- { ///////////// Find the region into SAR input images (same geometry for SAR after coregistration ////////////
- ImageSARRegionType inputSARRequestedRegion = OutputRegionToInputRegion(outputRequestedRegion, true);
+ ///////////// Find the region into SAR input images (same geometry for SAR after coregistration ////////////
+ ImageSARRegionType inputSARRequestedRegion = OutputRegionToInputRegion(outputRequestedRegion, true);
- ImageSARPointer masterPtr = const_cast< ImageSARType * >( this->GetMasterInput() );
- ImageSARPointer slavePtr = const_cast< ImageSARType * >( this->GetSlaveInput() );
+ ImageSARPointer masterPtr = const_cast< ImageSARType * >( this->GetMasterInput() );
+ ImageSARPointer slavePtr = const_cast< ImageSARType * >( this->GetSlaveInput() );
- masterPtr->SetRequestedRegion(inputSARRequestedRegion);
- slavePtr->SetRequestedRegion(inputSARRequestedRegion);
+ masterPtr->SetRequestedRegion(inputSARRequestedRegion);
+ slavePtr->SetRequestedRegion(inputSARRequestedRegion);
- ///////////// Find the region into topographic phase image (if needed) /////////////
- if (this->GetTopographicPhaseInput() != nullptr)
- {
- bool sameGeoAsMainOutput;
- ImagePhaseRegionType phaseRequestedRegion = OutputRegionToInputPhaseRegion(outputRequestedRegion,
- sameGeoAsMainOutput, true);
- ImagePhasePointer phasePtr = const_cast< ImagePhaseType * >( this->GetTopographicPhaseInput() );
- phasePtr->SetRequestedRegion(phaseRequestedRegion);
- }
- }
- else
+ ///////////// Find the region into topographic phase image (if needed) /////////////
+ if (this->GetTopographicPhaseInput() != nullptr)
{
- // No region requested
- ImageSARRegionType inputSARRequestedRegion;
- ImageSARIndexType index;
- index.Fill(0);
- ImageSARSizeType size;
- size.Fill(0);
- inputSARRequestedRegion.SetIndex(index);
- inputSARRequestedRegion.SetSize(size);
-
- ImageSARPointer masterPtr = const_cast< ImageSARType * >( this->GetMasterInput() );
- ImageSARPointer slavePtr = const_cast< ImageSARType * >( this->GetSlaveInput() );
-
- masterPtr->SetRequestedRegion(inputSARRequestedRegion);
- slavePtr->SetRequestedRegion(inputSARRequestedRegion);
-
-
- ImagePhaseRegionType phaseRequestedRegion;
- ImagePhaseIndexType indexPhase;
- indexPhase.Fill(0);
- ImagePhaseSizeType sizePhase;
- sizePhase.Fill(0);
- phaseRequestedRegion.SetIndex(indexPhase);
- phaseRequestedRegion.SetSize(sizePhase);
-
+ bool sameGeoAsMainOutput;
+ ImagePhaseRegionType phaseRequestedRegion = OutputRegionToInputPhaseRegion(outputRequestedRegion,
+ sameGeoAsMainOutput, true);
ImagePhasePointer phasePtr = const_cast< ImagePhaseType * >( this->GetTopographicPhaseInput() );
phasePtr->SetRequestedRegion(phaseRequestedRegion);
}
}
-
-
- /**
- * Method EnlargeOutputRequestedRegion
- */
- template<class TImageSAR, class TImageDEM, class TImagePhase, class TImageOut>
- void
- SARInterferogramImageFilter< TImageSAR, TImageDEM, TImagePhase, TImageOut >
- ::EnlargeOutputRequestedRegion( itk::DataObject *itkNotUsed(output) )
- {
- // This filter requires all of the optionnal output image (Nth Output = 1).
- if (m_UseDEMGeoAsOutput)
- {
- if ( this->itk::ProcessObject::GetOutput(1) )
- {
- this->itk::ProcessObject::GetOutput(1)->SetRequestedRegionToLargestPossibleRegion();
- }
- }
- }
-
-
- /**
- * Method BeforeThreadedGenerateData
- */
- template<class TImageSAR, class TImageDEM, class TImagePhase, class TImageOut>
- void
- SARInterferogramImageFilter< TImageSAR, TImageDEM, TImagePhase, TImageOut >
- ::BeforeThreadedGenerateData()
- {
- // Allocation for the optionnal image (if needed)
- if (m_UseDEMGeoAsOutput)
- {
- // Allocate optionnal output
- ImageOutPointer outputOptionnalPtr = this->GetOutputIntoDEMGeo();
- outputOptionnalPtr->SetBufferedRegion(outputOptionnalPtr->GetRequestedRegion() );
- outputOptionnalPtr->Allocate();
- }
- }
-
- /**
- * Method AfterThreadedGenerateData
- */
- template<class TImageSAR, class TImageDEM, class TImagePhase, class TImageOut>
- void
- SARInterferogramImageFilter< TImageSAR, TImageDEM, TImagePhase, TImageOut >
- ::AfterThreadedGenerateData()
- {
- }
/**
@@ -951,470 +749,332 @@ namespace otb
ImageOutPixelType outPixel;
outPixel.Reserve(4);
- // Two possible outputs => Two kinds of streaming and calculation :
- // m_OutputCounter = 0 => Main output and main calculation
- // m_OutputCounter = 1 => Optionnal output (ortho geo) and copy of accumulation into ortho geo to ortho
- // image
- if (m_OutputCounter == 1)
- {
- // Compute corresponding input region for SAR images
- ImageSARRegionType inputRegionForThread = OutputRegionToInputRegion(outputRegionForThread, true);
+ // Compute corresponding input region for SAR images
+ ImageSARRegionType inputRegionForThread = OutputRegionToInputRegion(outputRegionForThread, true);
- // Compute corresponding input region for Topographic phase image
- bool sameGeoAsMainOutput; // If true => Geo TopoPhase == main output geo else Geo TopoPhase = SAR Geo
- ImagePhaseRegionType inputPhaseRegionForThread = OutputRegionToInputPhaseRegion(outputRegionForThread,
- sameGeoAsMainOutput, true);
- // Define/declare an iterator that will walk the input regions for this
- // thread.
- InputSARIterator inMasterIt(this->GetMasterInput(), inputRegionForThread);
- InputSARIterator inSlaveIt(this->GetSlaveInput(), inputRegionForThread);
-
- InputPhaseIterator inTopoPhaseIt(this->GetTopographicPhaseInput(), inputPhaseRegionForThread);
-
- // Define/declare an iterator that will walk the output region for this
- // thread.
- OutputIterator outIt(this->GetOutput(), outputRegionForThread);
+ // Compute corresponding input region for Topographic phase image
+ bool sameGeoAsMainOutput; // If true => Geo TopoPhase == main output geo else Geo TopoPhase = SAR Geo
+ ImagePhaseRegionType inputPhaseRegionForThread = OutputRegionToInputPhaseRegion(outputRegionForThread,
+ sameGeoAsMainOutput, true);
+ // Define/declare an iterator that will walk the input regions for this
+ // thread.
+ InputSARIterator inMasterIt(this->GetMasterInput(), inputRegionForThread);
+ InputSARIterator inSlaveIt(this->GetSlaveInput(), inputRegionForThread);
+
+ InputPhaseIterator inTopoPhaseIt(this->GetTopographicPhaseInput(), inputPhaseRegionForThread);
+
+ // Define/declare an iterator that will walk the output region for this
+ // thread.
+ OutputIterator outIt(this->GetOutput(), outputRegionForThread);
- inMasterIt.GoToBegin();
- inSlaveIt.GoToBegin();
- inTopoPhaseIt.GoToBegin();
- outIt.GoToBegin();
-
- // Get the number of column of the region
- const ImageOutSizeType & outputRegionForThreadSize = outputRegionForThread.GetSize();
- unsigned int nbCol = static_cast<unsigned int>(outputRegionForThreadSize[0]);
+ inMasterIt.GoToBegin();
+ inSlaveIt.GoToBegin();
+ inTopoPhaseIt.GoToBegin();
+ outIt.GoToBegin();
+
+ // Get the number of column of the region
+ const ImageOutSizeType & outputRegionForThreadSize = outputRegionForThread.GetSize();
+ unsigned int nbCol = static_cast<unsigned int>(outputRegionForThreadSize[0]);
- // Allocate Arrays for accumulations (size = nbCol)
- double *complexMulConjTab_Re = new double[nbCol];
- double *complexMulConjTab_Im = new double[nbCol];
- double *complexMulConjTab_Mod = new double[nbCol];
- int *isDataCounter = new int[nbCol];
-
- Point2DType groundPoint;
- itk::ContinuousIndex<double, 2> groundContinousIndex;
-
- //// Initialize the indLastL and indFirstL (for spending our averaging with a margin) ////
- int indFirstL = 0;
- int indLastL = m_MLAzi+ m_MarginAzi;
- int marginAzi = (int) m_MarginAzi;
+ // Allocate Arrays for accumulations (size = nbCol)
+ double *complexMulConjTab_Re = new double[nbCol];
+ double *complexMulConjTab_Im = new double[nbCol];
+ double *complexMulConjTab_Mod = new double[nbCol];
+ int *isDataCounter = new int[nbCol];
+
+ Point2DType groundPoint;
+ itk::ContinuousIndex<double, 2> groundContinousIndex;
+
+ //// Initialize the indLastL and indFirstL (for spending our averaging with a margin) ////
+ int indFirstL = 0;
+ int indLastL = m_MLAzi+ m_MarginAzi;
+ int marginAzi = (int) m_MarginAzi;
- if (inMasterIt.GetIndex()[1] >= marginAzi && marginAzi != 0)
- {
- indFirstL = -marginAzi;
- }
+ if (inMasterIt.GetIndex()[1] >= marginAzi && marginAzi != 0)
+ {
+ indFirstL = -marginAzi;
+ }
- if ((inMasterIt.GetIndex()[1] + (indLastL - indFirstL)) > (m_nbLinesSAR-1))
- {
- indLastL = (m_nbLinesSAR-1) + indFirstL - inMasterIt.GetIndex()[1];
- }
+ if ((inMasterIt.GetIndex()[1] + (indLastL - indFirstL)) > (m_nbLinesSAR-1))
+ {
+ indLastL = (m_nbLinesSAR-1) + indFirstL - inMasterIt.GetIndex()[1];
+ }
- bool backTopreviousLines = false;
- bool backTopreviousColunms = false;
+ bool backTopreviousLines = false;
+ bool backTopreviousColunms = false;
- //////// For each Line /////////
- while ( !inMasterIt.IsAtEnd() && !inSlaveIt.IsAtEnd() && !inTopoPhaseIt.IsAtEnd() && !outIt.IsAtEnd())
+ //////// For each Line /////////
+ while ( !inMasterIt.IsAtEnd() && !inSlaveIt.IsAtEnd() && !inTopoPhaseIt.IsAtEnd() && !outIt.IsAtEnd())
+ {
+ // reinitialize
+ for (unsigned int j = 0; j < nbCol; j++)
{
- // reinitialize
- for (unsigned int j = 0; j < nbCol; j++)
- {
- complexMulConjTab_Re[j] = 0;
- complexMulConjTab_Im[j] = 0;
- complexMulConjTab_Mod[j] = 0;
- isDataCounter[j] = 0;
- }
+ complexMulConjTab_Re[j] = 0;
+ complexMulConjTab_Im[j] = 0;
+ complexMulConjTab_Mod[j] = 0;
+ isDataCounter[j] = 0;
+ }
- if (backTopreviousLines)
+ if (backTopreviousLines)
+ {
+ //// Previous index in azimut (for averaging) ////
+ if (inMasterIt.GetIndex()[1] >= marginAzi && marginAzi != 0)
{
- //// Previous index in azimut (for averaging) ////
- if (inMasterIt.GetIndex()[1] >= marginAzi && marginAzi != 0)
- {
- indFirstL = -marginAzi;
- }
+ indFirstL = -marginAzi;
+ }
- ImageSARIndexType indSARL;
- indSARL[0] = inMasterIt.GetIndex()[0];
- indSARL[1] = inMasterIt.GetIndex()[1] - (2*marginAzi);
+ ImageSARIndexType indSARL;
+ indSARL[0] = inMasterIt.GetIndex()[0];
+ indSARL[1] = inMasterIt.GetIndex()[1] - (2*marginAzi);
- ImagePhaseIndexType indPhaL;
- indPhaL[0] = inTopoPhaseIt.GetIndex()[0];
- indPhaL[1] = inTopoPhaseIt.GetIndex()[1] - (2*marginAzi);
+ ImagePhaseIndexType indPhaL;
+ indPhaL[0] = inTopoPhaseIt.GetIndex()[0];
+ indPhaL[1] = inTopoPhaseIt.GetIndex()[1] - (2*marginAzi);
- if ((indSARL[1] + indLastL - indFirstL) > (m_nbLinesSAR-1))
- {
- indLastL = (m_nbLinesSAR - 1) - indSARL[1] + indFirstL;
- }
+ if ((indSARL[1] + indLastL - indFirstL) > (m_nbLinesSAR-1))
+ {
+ indLastL = (m_nbLinesSAR - 1) - indSARL[1] + indFirstL;
+ }
- // Check previous index
- if (indSARL[1] >= 0 && indPhaL[1] >= 0 && marginAzi != 0)
+ // Check previous index
+ if (indSARL[1] >= 0 && indPhaL[1] >= 0 && marginAzi != 0)
+ {
+ // Previous index inputs
+ inMasterIt.SetIndex(indSARL);
+ inSlaveIt.SetIndex(indSARL);
+ // If topographic phase has SAR geo
+ if (!sameGeoAsMainOutput)
{
- // Previous index inputs
- inMasterIt.SetIndex(indSARL);
- inSlaveIt.SetIndex(indSARL);
- // If topographic phase has SAR geo
- if (!sameGeoAsMainOutput)
- {
- inTopoPhaseIt.SetIndex(indPhaL);
- }
-
- indFirstL = -marginAzi;
+ inTopoPhaseIt.SetIndex(indPhaL);
}
+
+ indFirstL = -marginAzi;
}
+ }
- backTopreviousLines = true;
+ backTopreviousLines = true;
- // GenerateInputRequestedRegion function requires an input region with
- // In_nbLine = m_Averaging * Out_nbLine (Averaing in function of ML Factors and Margins)
- for (int i = indFirstL ; i < indLastL; i++)
- {
- inMasterIt.GoToBeginOfLine();
- inSlaveIt.GoToBeginOfLine();
- inTopoPhaseIt.GoToBeginOfLine();
- outIt.GoToBeginOfLine();
+ // GenerateInputRequestedRegion function requires an input region with
+ // In_nbLine = m_Averaging * Out_nbLine (Averaing in function of ML Factors and Margins)
+ for (int i = indFirstL ; i < indLastL; i++)
+ {
+ inMasterIt.GoToBeginOfLine();
+ inSlaveIt.GoToBeginOfLine();
+ inTopoPhaseIt.GoToBeginOfLine();
+ outIt.GoToBeginOfLine();
- int colCounter = 0;
- backTopreviousColunms = false;
+ int colCounter = 0;
+ backTopreviousColunms = false;
- //// Initialize the indLastC and indFirstC (for spending average with a margin) ////
- int indFirstC = 0;
- int indLastC = m_MLRan + m_MarginRan;
- int marginRan = (int) m_MarginRan;
+ //// Initialize the indLastC and indFirstC (for spending average with a margin) ////
+ int indFirstC = 0;
+ int indLastC = m_MLRan + m_MarginRan;
+ int marginRan = (int) m_MarginRan;
- if (inMasterIt.GetIndex()[0] >= marginRan && marginRan != 0)
- {
- indFirstC = -marginRan;
- }
+ if (inMasterIt.GetIndex()[0] >= marginRan && marginRan != 0)
+ {
+ indFirstC = -marginRan;
+ }
- if ((inMasterIt.GetIndex()[0] + (indLastC - indFirstC)) > (m_nbColSAR-1))
- {
- indLastC = (m_nbColSAR-1) + indFirstC - inMasterIt.GetIndex()[0];
- }
+ if ((inMasterIt.GetIndex()[0] + (indLastC - indFirstC)) > (m_nbColSAR-1))
+ {
+ indLastC = (m_nbColSAR-1) + indFirstC - inMasterIt.GetIndex()[0];
+ }
- ///////// For each column /////////
- while (!inMasterIt.IsAtEndOfLine() && !inSlaveIt.IsAtEndOfLine() && !inTopoPhaseIt.IsAtEndOfLine()
- && !outIt.IsAtEndOfLine())
- {
+ ///////// For each column /////////
+ while (!inMasterIt.IsAtEndOfLine() && !inSlaveIt.IsAtEndOfLine() && !inTopoPhaseIt.IsAtEndOfLine()
+ && !outIt.IsAtEndOfLine())
+ {
- if (backTopreviousColunms)
+ if (backTopreviousColunms)
+ {
+ //// Previous index in range (for averaging) ////
+ if (inMasterIt.GetIndex()[0] >= marginRan && marginRan != 0)
{
- //// Previous index in range (for averaging) ////
- if (inMasterIt.GetIndex()[0] >= marginRan && marginRan != 0)
- {
- indFirstC = -marginRan;
- }
+ indFirstC = -marginRan;
+ }
- ImageSARIndexType indSARC;
- indSARC[0] = inMasterIt.GetIndex()[0] - (2*marginRan);
- indSARC[1] = inMasterIt.GetIndex()[1];
+ ImageSARIndexType indSARC;
+ indSARC[0] = inMasterIt.GetIndex()[0] - (2*marginRan);
+ indSARC[1] = inMasterIt.GetIndex()[1];
- ImagePhaseIndexType indPhaC;
- indPhaC[0] = inTopoPhaseIt.GetIndex()[0] - (2*marginRan);
- indPhaC[1] = inTopoPhaseIt.GetIndex()[1];
+ ImagePhaseIndexType indPhaC;
+ indPhaC[0] = inTopoPhaseIt.GetIndex()[0] - (2*marginRan);
+ indPhaC[1] = inTopoPhaseIt.GetIndex()[1];
- if ((indSARC[0] + indLastC - indFirstC) > (m_nbColSAR-1))
- {
- indLastC = (m_nbColSAR-1) - indSARC[0] + indFirstC;
- }
+ if ((indSARC[0] + indLastC - indFirstC) > (m_nbColSAR-1))
+ {
+ indLastC = (m_nbColSAR-1) - indSARC[0] + indFirstC;
+ }
- // Check previous index
- if (indSARC[0] >= 0 && indPhaC[0] >= 0 && marginRan != 0)
+ // Check previous index
+ if (indSARC[0] >= 0 && indPhaC[0] >= 0 && marginRan != 0)
+ {
+ // Previous index inputs
+ inMasterIt.SetIndex(indSARC);
+ inSlaveIt.SetIndex(indSARC);
+ // If topographic phase has SAR geo
+ if (!sameGeoAsMainOutput)
{
- // Previous index inputs
- inMasterIt.SetIndex(indSARC);
- inSlaveIt.SetIndex(indSARC);
- // If topographic phase has SAR geo
- if (!sameGeoAsMainOutput)
- {
- inTopoPhaseIt.SetIndex(indPhaC);
- }
-
- indFirstC = -marginRan;
+ inTopoPhaseIt.SetIndex(indPhaC);
}
-
+
+ indFirstC = -marginRan;
}
+
+ }
- backTopreviousColunms = true;
+ backTopreviousColunms = true;
- for (int k = indFirstC ; k < indLastC; k++)
+ for (int k = indFirstC ; k < indLastC; k++)
+ {
+ if (!inMasterIt.IsAtEndOfLine() && !inSlaveIt.IsAtEndOfLine())
{
- if (!inMasterIt.IsAtEndOfLine() && !inSlaveIt.IsAtEndOfLine())
- {
- // Check isData (grom topographic phase)
- if (inTopoPhaseIt.Get()[1] != 0)
- {
- // Complex Raw interferogram (master * conj(slave))
- double real_RawInterfero = (inMasterIt.Get().real()*inSlaveIt.Get().real() +
- inMasterIt.Get().imag()*inSlaveIt.Get().imag());
-
- double imag_RawInterfero = (inMasterIt.Get().imag()*inSlaveIt.Get().real() -
- inMasterIt.Get().real()*inSlaveIt.Get().imag());
-
- ///////////// Topographoc phase as complex number ////////////////
- double topoPhase = inTopoPhaseIt.Get()[0];
+ // Check isData (grom topographic phase)
+ if (inTopoPhaseIt.Get()[1] != 0)
+ {
+ // Complex Raw interferogram (master * conj(slave))
+ double real_RawInterfero = (inMasterIt.Get().real()*inSlaveIt.Get().real() +
+ inMasterIt.Get().imag()*inSlaveIt.Get().imag());
+
+ double imag_RawInterfero = (inMasterIt.Get().imag()*inSlaveIt.Get().real() -
+ inMasterIt.Get().real()*inSlaveIt.Get().imag());
+
+ ///////////// Topographoc phase as complex number ////////////////
+ double topoPhase = inTopoPhaseIt.Get()[0];
- if (m_ApproxDiapason)
- {
- int topoPhase_int = static_cast<int>(std::round(inTopoPhaseIt.Get()[0]))%
- 256;
- topoPhase = 2*M_PI*topoPhase_int / 256;
- }
-
- double complexTopoPhase_Re = std::cos(topoPhase);
- double complexTopoPhase_Im = std::sin(topoPhase);
+ double complexTopoPhase_Re = std::cos(topoPhase);
+ double complexTopoPhase_Im = std::sin(topoPhase);
- // Multiply the conj(complexTopoPhase) with complex raw interferogram
- double real_interfero = (real_RawInterfero * complexTopoPhase_Re +
- imag_RawInterfero * complexTopoPhase_Im);
- double imag_interfero = (imag_RawInterfero * complexTopoPhase_Re -
- real_RawInterfero * complexTopoPhase_Im);
+ // Multiply the conj(complexTopoPhase) with complex raw interferogram
+ double real_interfero = (real_RawInterfero * complexTopoPhase_Re +
+ imag_RawInterfero * complexTopoPhase_Im);
+ double imag_interfero = (imag_RawInterfero * complexTopoPhase_Re -
+ real_RawInterfero * complexTopoPhase_Im);
- ///////////// Accumulations ///////////////
- complexMulConjTab_Re[colCounter] += real_interfero;
+ ///////////// Accumulations ///////////////
+ complexMulConjTab_Re[colCounter] += real_interfero;
- complexMulConjTab_Im[colCounter] += imag_interfero;
+ complexMulConjTab_Im[colCounter] += imag_interfero;
- complexMulConjTab_Mod[colCounter] += sqrt((real_interfero*real_interfero) +
- (imag_interfero*imag_interfero));
-
- isDataCounter[colCounter] += 1;
-
- // Accumulation for ortho geo (only if k and i betwwen 0 and ml Factors)
- if (m_UseDEMGeoAsOutput && (i >= 0 && i < (int) m_MLAzi) && (k >= 0 && k < (int) m_MLRan))
- {
- // Determine the index on ortho image
- ossimEcefPoint ecefPt;
- ecefPt.x() = inTopoPhaseIt.Get()[2];
- ecefPt.y() = inTopoPhaseIt.Get()[3];
- ecefPt.z() = inTopoPhaseIt.Get()[4];
- ossimGpt gptPt(ecefPt);
- groundPoint[0] = gptPt.lon;
- groundPoint[1] = gptPt.lat;
-
- m_DEMPtr->TransformPhysicalPointToContinuousIndex(groundPoint,
- groundContinousIndex);
-
- int groundIndexC = static_cast<int>(std::round(groundContinousIndex[0]));
- int groundIndexL = static_cast<int>(std::round(groundContinousIndex[1]));
-
- int groundIndex = groundIndexL*m_nbColDEM + groundIndexC;
-
- // Accumulate into OrthoAccByThread at threadId and at groundIndex
- /*m_OrthoRealAccByThread[threadId][groundIndex] += real_interfero;
- m_OrthoImagAccByThread[threadId][groundIndex] += imag_interfero;
- m_OrthoModAccByThread[threadId][groundIndex] += sqrt((real_interfero*
- real_interfero) +
- (imag_interfero*
- imag_interfero));
-
- // Increment counter
- ++m_OrthoCounterByThread[threadId];*/
-
- m_Mutex->Lock();
-
- m_OrthoRealAcc[groundIndex] += real_interfero;
- m_OrthoImagAcc[groundIndex] += imag_interfero;
- m_OrthoModAcc[groundIndex] += sqrt((real_interfero*real_interfero) +
- (imag_interfero*imag_interfero));
-
- ++m_OrthoCounter[groundIndex];
-
- m_Mutex->Unlock();
-
-
- }
- }
-
- // Next colunm inputs
- ++inMasterIt;
- ++inSlaveIt;
- // If topographic phase has SAR geo
- if (!sameGeoAsMainOutput)
- {
- ++inTopoPhaseIt;
- }
+ complexMulConjTab_Mod[colCounter] += sqrt((real_interfero*real_interfero) +
+ (imag_interfero*imag_interfero));
+
+ isDataCounter[colCounter] += 1;
}
- } // End for k (indFirstC to indLastC)
+
+ // Next colunm inputs
+ ++inMasterIt;
+ ++inSlaveIt;
+ // If topographic phase has SAR geo
+ if (!sameGeoAsMainOutput)
+ {
+ ++inTopoPhaseIt;
+ }
+ }
+ } // End for k (indFirstC to indLastC)
// Estiamte and Assigne outIt
//if (i == (m_MLAzi-1))
- if (i == (indLastL-1))
+ if (i == (indLastL-1))
+ {
+ // Check if Data
+ if (isDataCounter[colCounter] >= 1)
{
- // Check if Data
- if (isDataCounter[colCounter] >= 1)
- {
- ///////////// Estimations of amplitude, phase and coherency ///////////////
- double mod_Acc = sqrt(complexMulConjTab_Re[colCounter]*complexMulConjTab_Re[colCounter]
- + complexMulConjTab_Im[colCounter]*
- complexMulConjTab_Im[colCounter]);
+ ///////////// Estimations of amplitude, phase and coherency ///////////////
+ double mod_Acc = sqrt(complexMulConjTab_Re[colCounter]*complexMulConjTab_Re[colCounter]
+ + complexMulConjTab_Im[colCounter]*
+ complexMulConjTab_Im[colCounter]);
- int countPixel = (m_MLRan + 2*m_MarginRan) * (m_MLAzi + 2*m_MarginAzi);
+ int countPixel = (m_MLRan + 2*m_MarginRan) * (m_MLAzi + 2*m_MarginAzi);
- // Amplitude
- outPixel[0] = m_Gain * sqrt((complexMulConjTab_Mod[colCounter]/(countPixel)));
+ // Amplitude
+ outPixel[0] = m_Gain * sqrt((complexMulConjTab_Mod[colCounter]/(countPixel)));
- // Phase
- outPixel[1] = std::atan2(complexMulConjTab_Im[colCounter],
- complexMulConjTab_Re[colCounter]);
-
- // Mod 2*Pi
- outPixel[1] = outPixel[1]-(2*M_PI)*floor(outPixel[1]/(2*M_PI));
+ // Phase
+ outPixel[1] = std::atan2(complexMulConjTab_Im[colCounter],
+ complexMulConjTab_Re[colCounter]);
- if (m_ApproxDiapason)
- {
- outPixel[1] = std::atan2(complexMulConjTab_Im[colCounter],
- complexMulConjTab_Re[colCounter])* 128 * M_1_PI;
+ // Mod 2*Pi
+ outPixel[1] = outPixel[1]-(2*M_PI)*floor(outPixel[1]/(2*M_PI));
- outPixel[1] = outPixel[1]-(256)*floor(outPixel[1]/(256));
-
- outPixel[1] = std::round(outPixel[1]);
- }
-
- // Coherency
- if (complexMulConjTab_Mod[colCounter] != 0)
- {
- outPixel[2] = mod_Acc / complexMulConjTab_Mod[colCounter];
- }
- else
- {
- outPixel[2] = 0;
- }
-
- // isData set to 1
- outPixel[3] = 1;
-
- if (outPixel[0] == 0 && outPixel[1] == 0 && outPixel[2] == 0)
- {
- outPixel[3] = 0;
- }
+ // Coherency
+ if (complexMulConjTab_Mod[colCounter] != 0)
+ {
+ outPixel[2] = mod_Acc / complexMulConjTab_Mod[colCounter];
}
else
{
- outPixel[0] = 0;
- outPixel[1] = 0;
outPixel[2] = 0;
- outPixel[3] = 0;
}
-
- outIt.Set(outPixel);
- progress.CompletedPixel();
- }
- // Next colunm output
- ++outIt;
- // If topographic phase has Main output geo
- if (sameGeoAsMainOutput)
+ // isData set to 1
+ outPixel[3] = 1;
+
+ if (outPixel[0] == 0 && outPixel[1] == 0 && outPixel[2] == 0)
+ {
+ outPixel[3] = 0;
+ }
+ }
+ else
{
- ++inTopoPhaseIt;
+ outPixel[0] = 0;
+ outPixel[1] = 0;
+ outPixel[2] = 0;
+ outPixel[3] = 0;
}
- colCounter++;
+
+ outIt.Set(outPixel);
+ progress.CompletedPixel();
}
-
- // Next line intputs
- inMasterIt.NextLine();
- inSlaveIt.NextLine();
- // If topographic phase has SAR geo
- if (!sameGeoAsMainOutput)
+
+ // Next colunm output
+ ++outIt;
+ // If topographic phase has Main output geo
+ if (sameGeoAsMainOutput)
{
- inTopoPhaseIt.NextLine();
+ ++inTopoPhaseIt;
}
-
- } // End for i (firstL to lastL)
-
- // Next line output
- outIt.NextLine();
- // If topographic phase has Main output geo
- if (sameGeoAsMainOutput)
+ colCounter++;
+ }
+
+ // Next line intputs
+ inMasterIt.NextLine();
+ inSlaveIt.NextLine();
+ // If topographic phase has SAR geo
+ if (!sameGeoAsMainOutput)
{
inTopoPhaseIt.NextLine();
}
-
+
+ } // End for i (firstL to lastL)
+
+ // Next line output
+ outIt.NextLine();
+ // If topographic phase has Main output geo
+ if (sameGeoAsMainOutput)
+ {
+ inTopoPhaseIt.NextLine();
}
+
+ }
delete [] complexMulConjTab_Re;
delete [] complexMulConjTab_Im;
delete [] complexMulConjTab_Mod;
delete [] isDataCounter;
- }
-
- //////////// Optionnal output ////////////////
- // Copy for optionnal output (ortho image)
- if (m_OutputCounter == 2 && m_UseDEMGeoAsOutput)
- {
- // Streming is on main output or orhto geo is very different => to guarantee Thread-safety only
- // threadId = 0 makes the copy (the whole copy)
- if (threadId == 0)
- {
- int sizeDEM = m_nbColDEM * m_nbLinesDEM;
- ImageOutIndexType outDEMIndex;
-
- for (int indOrtho = 0; indOrtho < sizeDEM; indOrtho++)
- {
- // Index
- outDEMIndex[0] = indOrtho % m_nbColDEM;
- outDEMIndex[1] = indOrtho / m_nbColDEM; // ninteger division
-
- // Amplitude, Phase and coherency
- if (m_OrthoCounter[indOrtho] != 0)
- {
- double mod_Acc = sqrt(m_OrthoRealAcc[indOrtho]*m_OrthoRealAcc[indOrtho]
- + m_OrthoImagAcc[indOrtho]* m_OrthoImagAcc[indOrtho]);
-
-
- // Amplitude
- outPixel[0] = m_Gain * sqrt(m_OrthoModAcc[indOrtho]/m_OrthoCounter[indOrtho]);
-
- // Phase
- outPixel[1] = std::atan2(m_OrthoImagAcc[indOrtho],
- m_OrthoRealAcc[indOrtho]);
-
- // Mod 2*Pi
- outPixel[1] = outPixel[1]-(2*M_PI)*floor(outPixel[1]/(2*M_PI));
-
- if (m_ApproxDiapason)
- {
- outPixel[1] = std::atan2(m_OrthoImagAcc[indOrtho],
- m_OrthoRealAcc[indOrtho])* 128 * M_1_PI;
-
- //outPixel[1] = outPixel[1]-(256)*floor(outPixel[1]/(256));
-
- outPixel[1] = std::round(outPixel[1]);
- }
-
- // Coherency
- if (m_OrthoModAcc[indOrtho] != 0)
- {
- outPixel[2] = mod_Acc / m_OrthoModAcc[indOrtho];
- }
- else
- {
- outPixel[2] = 0;
- }
-
- // isData set to 1
- outPixel[3] = 1;
-
- }
- else
- {
- outPixel[0] = 0;
- outPixel[1] = 0;
- outPixel[2] = 0;
- outPixel[3] = 0;
- }
+
- this->GetOutputIntoDEMGeo()->SetPixel(outDEMIndex, outPixel);
- }
- }
- }
}
-} /*namespace otb*/
+ } /*namespace otb*/
#endif
diff --git a/python_src/diapOTB.py b/python_src/diapOTB.py
index e59c71d..9255952 100644
--- a/python_src/diapOTB.py
+++ b/python_src/diapOTB.py
@@ -140,7 +140,6 @@ if __name__ == "__main__":
ml_geoGrid_range = ml_range
ml_geoGrid_azimut = ml_azimut
gain_interfero = dict_DInSAR['parameter']['Interferogram_gain']
- activateOrthoInterferogram = dict_DInSAR['parameter']['Interferogram_ortho']
if (geoGrid_threshold < 0) or (geoGrid_threshold > 1) :
print("Wrong Threshold for fine deformation grid")
@@ -413,7 +412,6 @@ if __name__ == "__main__":
######## SARRobustInterferogram Application (interf step) #######
print("\n SARRobustInterferogram Application \n")
interferogram = "interferogram.tif"
- interferogram_ortho = "interferogram_ortho.tif"
appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram")
appInterferogram.SetParameterString("insarmaster", master_Image)
appInterferogram.SetParameterString("incoregistratedslave", os.path.join(output_dir, slave_Image_CoRe))
@@ -427,8 +425,4 @@ if __name__ == "__main__":
appInterferogram.SetParameterFloat("gain", gain_interfero)
appInterferogram.SetParameterString("ram", "4000")
appInterferogram.SetParameterString("out", os.path.join(output_dir, interferogram))
- if activateOrthoInterferogram :
- appInterferogram.SetParameterString("ortho", "true")
- appInterferogram.SetParameterString("indem", dem)
- appInterferogram.SetParameterString("outopt", os.path.join(output_dir, interferogram_ortho))
appInterferogram.ExecuteAndWriteOutput()
diff --git a/python_src/diapOTB_S1IW.py b/python_src/diapOTB_S1IW.py
index 352de09..0dcf6ed 100644
--- a/python_src/diapOTB_S1IW.py
+++ b/python_src/diapOTB_S1IW.py
@@ -164,7 +164,6 @@ if __name__ == "__main__":
ml_geoGrid_range = ml_range
ml_geoGrid_azimut = ml_azimut
gain_interfero = dict_DInSAR['parameter']['Interferogram_gain']
- activateOrthoInterferogram = dict_DInSAR['parameter']['Interferogram_ortho']
# esd loop
esd_AutoMode = False # automatic mode to apply a threshold inside the esd loop
esd_NbIter = 0
@@ -591,8 +590,7 @@ if __name__ == "__main__":
######## SARRobustInterferogram Application (interf step) #######
print("\n SARRobustInterferogram Application \n")
interferogram_path = os.path.join(burst_dir, "interferogram" + "_burst" + str(burstId) + ".tif")
- interferogram_ortho_path = os.path.join(burst_dir, "interferogram_ortho" + "_burst" + str(burstId) + ".tif")
-
+
if esd_NbIter > 0 :
esd_dir = os.path.join(burst_dir, "esd")
@@ -600,8 +598,7 @@ if __name__ == "__main__":
os.makedirs(esd_dir)
interferogram_path = os.path.join(esd_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif")
- interferogram_ortho_path = os.path.join(esd_dir, "interferogram_ortho" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif")
-
+
appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram")
appInterferogram.SetParameterString("insarmaster", os.path.join(burst_dir, burstM))
appInterferogram.SetParameterString("incoregistratedslave", os.path.join(burst_dir, burstCoReRerampS))
@@ -617,10 +614,6 @@ if __name__ == "__main__":
appInterferogram.SetParameterFloat("gain", gain_interfero)
appInterferogram.SetParameterString("ram", "4000")
appInterferogram.SetParameterString("out", interferogram_path)
- if activateOrthoInterferogram :
- appInterferogram.SetParameterString("ortho", "true")
- appInterferogram.SetParameterString("indem", dem)
- appInterferogram.SetParameterString("outopt", interferogram_ortho_path)
appInterferogram.ExecuteAndWriteOutput()
list_of_Interferograms.append(interferogram_path)
@@ -749,7 +742,6 @@ if __name__ == "__main__":
######## SARRobustInterferogram Application (interf step) #######
print("\n SARRobustInterferogram Application \n")
interferogram_path = os.path.join(esd_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(iter_esd) + ".tif")
- interferogram_ortho_path = os.path.join(esd_dir, "interferogram_ortho" + "_burst" + str(burstId) + "_iter" + str(iter_esd) + ".tif")
appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram")
appInterferogram.SetParameterString("insarmaster", os.path.join(burst_dir, burstM))
@@ -768,10 +760,6 @@ if __name__ == "__main__":
appInterferogram.SetParameterFloat("gain", gain_interfero)
appInterferogram.SetParameterString("ram", "4000")
appInterferogram.SetParameterString("out", interferogram_path)
- if activateOrthoInterferogram :
- appInterferogram.SetParameterString("ortho", "true")
- appInterferogram.SetParameterString("indem", dem)
- appInterferogram.SetParameterString("outopt", interferogram_ortho_path)
appInterferogram.ExecuteAndWriteOutput()
# Store the new interferogram paths
--
GitLab