diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 29d85143800acde0b28afb991052e04f5e06a882..151f0f4f5952232c9698010e7e135452e4339a6d 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -103,3 +103,8 @@ OTB_CREATE_APPLICATION(NAME SARMetadataCorrection
SOURCES otbSARMetadataCorrection.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
)
+
+OTB_CREATE_APPLICATION(NAME SAROrthoInterferogram
+ SOURCES otbSAROrthoInterferogram.cxx
+ LINK_LIBRARIES ${${otb-module}_LIBRARIES}
+)
diff --git a/app/otbSAROrthoInterferogram.cxx b/app/otbSAROrthoInterferogram.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8943cf069370caf37f00679a9b08b54908ebaa31
--- /dev/null
+++ b/app/otbSAROrthoInterferogram.cxx
@@ -0,0 +1,173 @@
+/*
+ * Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "otbWrapperApplication.h"
+#include "otbWrapperApplicationFactory.h"
+
+#include "otbSARCompensatedComplexImageFilter.h"
+#include "otbSARGroupedByOrthoImageFilter.h"
+
+#include "otbWrapperOutputImageParameter.h"
+#include "otbWrapperTypes.h"
+
+#include <iostream>
+#include <string>
+#include <fstream>
+
+namespace otb
+{
+
+ namespace Wrapper
+ {
+
+ class SAROrthoInterferogram : public Application
+ {
+ public:
+ typedef SAROrthoInterferogram Self;
+ typedef itk::SmartPointer<Self> Pointer;
+
+ itkNewMacro(Self);
+ itkTypeMacro(SAROrthoInterferogram, otb::Wrapper::Application);
+
+ // Filters
+ typedef otb::SARCompensatedComplexImageFilter<ComplexFloatImageType, FloatVectorImageType, FloatVectorImageType> CompensatedComplexFilterType;
+ typedef otb::SARGroupedByOrthoImageFilter<FloatVectorImageType, FloatImageType, FloatVectorImageType> GroupedByOrthoFilterType;
+
+ private:
+ void DoInit() override
+ {
+ SetName("SAROrthoInterferogram");
+ SetDescription("Interferogram into ground geometry between two SAR images.");
+
+ SetDocLongDescription("This application builds the interferogram into ground geometry between two"
+ " SAR images.");
+
+ //Optional descriptors
+ SetDocLimitations("Only Sentinel 1 (IW and StripMap mode) and Cosmo products are supported for now.");
+ SetDocAuthors("OTB-Team");
+ SetDocSeeAlso(" ");
+ AddDocTag(Tags::SAR);
+ AddDocTag("DiapOTB");
+
+ //Parameter declarations
+ AddParameter(ParameterType_InputImage, "insarslave", "Input SAR Slave image (Coregistrated image)");
+ SetParameterDescription("insarslave", "Input SAR Slave image (Coregistrated image).");
+
+ AddParameter(ParameterType_InputImage, "insarmaster", "Input SAR Master image");
+ SetParameterDescription("insarmaster", "Input SAR Master image.");
+
+ 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.");
+
+ AddParameter(ParameterType_InputImage, "indem", "Input DEM");
+ SetParameterDescription("indem", "GetGround geometry thanks to DEM (only metadata).");
+
+ AddParameter(ParameterType_Float, "gain", "Gain to apply for amplitude estimation");
+ SetParameterDescription("gain","Gain to apply for amplitude estimation");
+ SetDefaultParameterFloat("gain", 0.1);
+ SetMinimumParameterFloatValue("gain", 0);
+ MandatoryOff("gain");
+
+ AddParameter(ParameterType_OutputImage, "out", "Interferogram");
+ SetParameterDescription("out","Output Vector Image : Interferogram.");
+
+ AddRAMParameter();
+
+ SetDocExampleParameterValue("insarslave","s1b-s4-slc-vv-20160929t014610-20160929t014634-002277-003d71-002.tiff");
+ SetDocExampleParameterValue("insarmaster","s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002_SLC.tiff");
+ SetDocExampleParameterValue("incartmeanmaster","CartesianMaster.tiff");
+ SetDocExampleParameterValue("indem","S21E055.hgt");
+ SetDocExampleParameterValue("gain","0.1");
+ SetDocExampleParameterValue("out","s1b-s1a-s4-interferogram.tiff");
+ }
+
+ void DoUpdateParameters() override
+ {
+ // Nothing to do here : all parameters are independent
+ }
+
+ void DoExecute() override
+ {
+ // Get numeric parameters
+ double factor_gain = GetParameterFloat("gain");
+
+ otbAppLogINFO(<<"Gain Factor : "<<factor_gain);
+
+ /////////////////////////////////// Compensated Complex Filter ////////////////////////////////////////
+ // Instanciate the first filter
+ CompensatedComplexFilterType::Pointer filterCompensatedComplex = CompensatedComplexFilterType::New();
+ m_Ref.push_back(filterCompensatedComplex.GetPointer());
+
+ // Execute the Pipeline
+ ComplexFloatImageType::Pointer SARMasterPtr;
+ SARMasterPtr = GetParameterComplexFloatImage("insarmaster");
+
+ ComplexFloatImageType::Pointer SARSlavePtr;
+ SARSlavePtr = GetParameterComplexFloatImage("insarslave");
+
+ // Two Main inputs
+ filterCompensatedComplex->SetMasterInput(SARMasterPtr);
+ filterCompensatedComplex->SetSlaveInput(SARSlavePtr);
+
+ // One optionnal input
+ if (GetParameterByKey("topographicphase")->HasValue())
+ {
+ filterCompensatedComplex->SetTopographicPhaseInput(GetParameterImage("topographicphase"));
+ }
+
+
+ /////////////////////////////// GroupedBy Filter (for ground geometry) //////////////////////////////////
+ // Instanciate the second filter
+ FloatImageType::Pointer inputDEM = GetParameterFloatImage("indem");
+ GroupedByOrthoFilterType::Pointer filterGroupedBy = GroupedByOrthoFilterType::New();
+ m_Ref.push_back(filterGroupedBy.GetPointer());
+ filterGroupedBy->SetSARImageKeyWorList(SARMasterPtr->GetImageKeywordlist());
+ filterGroupedBy->SetDEMImagePtr(inputDEM);
+
+
+ // Execute the Pipeline
+ 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;
+
+ };
+
+ }
+
+}
+
+
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::SAROrthoInterferogram)
diff --git a/include/otbSARCompensatedComplexImageFilter.h b/include/otbSARCompensatedComplexImageFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..38a9490b58540364c9567db6c4f48a0b5912b42c
--- /dev/null
+++ b/include/otbSARCompensatedComplexImageFilter.h
@@ -0,0 +1,207 @@
+/*
+ * Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbSARCompensatedComplexImageFilter_h
+#define otbSARCompensatedComplexImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkSmartPointer.h"
+#include "itkPoint.h"
+
+#include "itkSimpleFastMutexLock.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageScanlineIterator.h"
+
+#include "otbImageKeywordlist.h"
+
+
+namespace otb
+{
+/** \class SARCompensatedComplexImageFilter
+ * \brief Creates an temporary image (to estimate conjugate multiplication) between two images.
+ *
+ * This filter built the conjugate multiplication between a master and a slave image.
+ * A topographic phase can be used to create this temporary image.
+ *
+ * The output is a vector image with real part, imag part and modulus of conjugate multiplication.
+ *
+ * \ingroup DiapOTBModule
+ */
+
+ template <typename TImageSAR, typename TImagePhase, typename TImageOut>
+ class ITK_EXPORT SARCompensatedComplexImageFilter :
+ public itk::ImageToImageFilter<TImageSAR,TImageOut>
+{
+public:
+
+ // Standard class typedefs
+ typedef SARCompensatedComplexImageFilter Self;
+ typedef itk::ImageToImageFilter<TImageSAR,TImageOut> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ // Method for creation through object factory
+ itkNewMacro(Self);
+ // Run-time type information
+ itkTypeMacro(SARCompensatedComplexImageFilter,ImageToImageFilter);
+
+ /** Typedef to image input type (master and slave reech images) : otb::Image (Complex) */
+ typedef TImageSAR ImageSARType;
+ /** Typedef to describe the inout image pointer type. */
+ typedef typename ImageSARType::Pointer ImageSARPointer;
+ typedef typename ImageSARType::ConstPointer ImageSARConstPointer;
+ /** Typedef to describe the inout image region type. */
+ typedef typename ImageSARType::RegionType ImageSARRegionType;
+ /** Typedef to describe the type of pixel and point for inout image. */
+ typedef typename ImageSARType::PixelType ImageSARPixelType;
+ typedef typename ImageSARType::PointType ImageSARPointType;
+ /** Typedef to describe the image index, size types and spacing for inout image. */
+ typedef typename ImageSARType::IndexType ImageSARIndexType;
+ typedef typename ImageSARType::IndexValueType ImageSARIndexValueType;
+ typedef typename ImageSARType::SizeType ImageSARSizeType;
+ typedef typename ImageSARType::SizeValueType ImageSARSizeValueType;
+ typedef typename ImageSARType::SpacingType ImageSARSpacingType;
+ typedef typename ImageSARType::SpacingValueType ImageSARSpacingValueType;
+
+ /** Typedef to image output type : otb::VectorImage */
+ typedef TImageOut ImageOutType;
+ /** Typedef to describe the output image pointer type. */
+ typedef typename ImageOutType::Pointer ImageOutPointer;
+ typedef typename ImageOutType::ConstPointer ImageOutConstPointer;
+ /** Typedef to describe the output image region type. */
+ typedef typename ImageOutType::RegionType ImageOutRegionType;
+ /** Typedef to describe the type of pixel and point for output image. */
+ typedef typename ImageOutType::PixelType ImageOutPixelType;
+ typedef typename ImageOutType::PointType ImageOutPointType;
+ /** Typedef to describe the image index, size types and spacing for output image. */
+ typedef typename ImageOutType::IndexType ImageOutIndexType;
+ typedef typename ImageOutType::IndexValueType ImageOutIndexValueType;
+ typedef typename ImageOutType::SizeType ImageOutSizeType;
+ typedef typename ImageOutType::SizeValueType ImageOutSizeValueType;
+ typedef typename ImageOutType::SpacingType ImageOutSpacingType;
+ typedef typename ImageOutType::SpacingValueType ImageOutSpacingValueType;
+
+
+ // Define Point2DType and Point3DType
+ using Point2DType = itk::Point<double,2>;
+ using Point3DType = itk::Point<double,3>;
+
+ // Typedef for topographic Phase
+ typedef TImagePhase ImagePhaseType;
+ typedef typename ImagePhaseType::Pointer ImagePhasePointer;
+ typedef typename ImagePhaseType::ConstPointer ImagePhaseConstPointer;
+ typedef typename ImagePhaseType::RegionType ImagePhaseRegionType;
+ typedef typename ImagePhaseType::IndexType ImagePhaseIndexType;
+ typedef typename ImagePhaseType::IndexValueType ImagePhaseIndexValueType;
+ typedef typename ImagePhaseType::SizeType ImagePhaseSizeType;
+ typedef typename ImagePhaseType::SizeValueType ImagePhaseSizeValueType;
+
+ // Typedef for iterators
+ typedef itk::ImageScanlineConstIterator< ImageSARType > InputSARIterator;
+ typedef itk::ImageScanlineConstIterator< ImagePhaseType > InputPhaseIterator;
+ typedef itk::ImageScanlineIterator< ImageOutType > OutputIterator;
+
+ typedef itk::SimpleFastMutexLock MutexType;
+
+
+ // Setter/Getter for inputs (SAR images and potentially topographic phase)
+ /** Connect one of the operands for interferometry : Master */
+ void SetMasterInput( const ImageSARType* image);
+
+ /** Connect one of the operands for interferometry : Coregistrated Slave */
+ void SetSlaveInput(const ImageSARType* image);
+
+ /** Connect one of the operands for interferometry : Coregistrated Slave */
+ void SetTopographicPhaseInput(const ImagePhaseType* imagePhaseTopo);
+
+ /** Get the inputs */
+ const ImageSARType* GetMasterInput() const;
+ const ImageSARType * GetSlaveInput() const;
+ const ImagePhaseType * GetTopographicPhaseInput() const;
+
+
+protected:
+ // Constructor
+ SARCompensatedComplexImageFilter();
+
+ // Destructor
+ virtual ~SARCompensatedComplexImageFilter() ITK_OVERRIDE;
+
+ // Print
+ void PrintSelf(std::ostream & os, itk::Indent indent) const ITK_OVERRIDE;
+
+ /** SARCompensatedComplexImageFilter produces an image which are into output geometry. The differences between
+ * output and input images can be the size of images and the dimensions.
+ * As such, SARCompensatedComplexImageFilter needs to provide an implementation for
+ * GenerateOutputInformation() in order to inform the pipeline execution model.
+ */
+ virtual void GenerateOutputInformation() ITK_OVERRIDE;
+
+ /** SARCompensatedComplexImageFilter 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.
+ * As such, TilesAnalysesImageFilter needs to provide an implementation for
+ * GenerateInputRequestedRegion() in order to inform the pipeline execution model.
+ * \sa ProcessObject::GenerateInputRequestedRegion() */
+ virtual void GenerateInputRequestedRegion() ITK_OVERRIDE;
+
+
+ /**
+ * SARCompensatedComplexImageFilterr can be implemented as a multithreaded filter.
+ * 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"
+ *
+ * \sa ImageToImageFilter::ThreadedGenerateData(),
+ * ImageToImageFilter::GenerateData() */
+ virtual void ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread,
+ itk::ThreadIdType threadId) ITK_OVERRIDE;
+
+ virtual void ThreadedGenerateDataWithoutTopographicPhase(const ImageOutRegionType& outputRegionForThread,
+ itk::ThreadIdType threadId);
+ virtual void ThreadedGenerateDataWithTopographicPhase(const ImageOutRegionType& outputRegionForThread,
+ itk::ThreadIdType threadId);
+
+
+ private:
+ SARCompensatedComplexImageFilter(const Self&); // purposely not implemented
+ void operator=(const Self &); // purposely not
+
+ // SAR Dimensions
+ int m_nbLinesSAR;
+ int m_nbColSAR;
+
+
+
+};
+
+} // End namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSARCompensatedComplexImageFilter.txx"
+#endif
+
+
+
+#endif
diff --git a/include/otbSARCompensatedComplexImageFilter.txx b/include/otbSARCompensatedComplexImageFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..8644b0b259151b46edcfc6278f1dfda46ed14674
--- /dev/null
+++ b/include/otbSARCompensatedComplexImageFilter.txx
@@ -0,0 +1,461 @@
+/*
+ * Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbSARCompensatedComplexImageFilter_txx
+#define otbSARCompensatedComplexImageFilter_txx
+
+#include "otbSARCompensatedComplexImageFilter.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageScanlineIterator.h"
+#include "itkProgressReporter.h"
+#include "itkNumericTraitsPointPixel.h"
+#include "itkContinuousIndex.h"
+
+#include "ossim/ossimSarSensorModel.h"
+
+#include <cmath>
+#include <algorithm>
+
+
+namespace otb
+{
+ /**
+ * Constructor with default initialization
+ */
+ template <class TImageSAR, class TImagePhase, class TImageOut>
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >::SARCompensatedComplexImageFilter()
+ : m_nbLinesSAR(0), m_nbColSAR(0)
+ {
+ // Inputs required and/or needed
+ this->SetNumberOfRequiredInputs(2);
+ this->SetNumberOfIndexedInputs(3);
+ }
+
+ /**
+ * Destructor
+ */
+ template <class TImageSAR, class TImagePhase, class TImageOut>
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >::~SARCompensatedComplexImageFilter()
+ {
+ }
+
+ /**
+ * Print
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ void
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::PrintSelf(std::ostream & os, itk::Indent indent) const
+ {
+ Superclass::PrintSelf(os, indent);
+ }
+
+ /**
+ * Set Master Image
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ void
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::SetMasterInput(const ImageSARType* image )
+ {
+ // Process object is not const-correct so the const casting is required.
+ this->SetNthInput(0, const_cast<ImageSARType *>(image));
+ }
+
+ /**
+ * Set Coregistrated Slave Image
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ void
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::SetSlaveInput(const ImageSARType* image )
+ {
+ // Process object is not const-correct so the const casting is required.
+ this->SetNthInput(1, const_cast<ImageSARType *>(image));
+ }
+
+ /**
+ * Set Topographic phase
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ void
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::SetTopographicPhaseInput(const ImagePhaseType* phaseTopoImage)
+ {
+ // Process object is not const-correct so the const casting is required.
+ this->SetNthInput(2, const_cast<ImagePhaseType *>(phaseTopoImage));
+ }
+
+ /**
+ * Get Master Image
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ const typename SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >::ImageSARType *
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::GetMasterInput() const
+ {
+ if (this->GetNumberOfInputs()<1)
+ {
+ return 0;
+ }
+ return static_cast<const ImageSARType *>(this->itk::ProcessObject::GetInput(0));
+ }
+
+ /**
+ * Get Coregistrated Slave Image
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ const typename SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >::ImageSARType *
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::GetSlaveInput() const
+ {
+ if (this->GetNumberOfInputs()<2)
+ {
+ return 0;
+ }
+ return static_cast<const ImageSARType *>(this->itk::ProcessObject::GetInput(1));
+ }
+
+ /**
+ * Get Phase Topographic
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ const typename SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >::ImagePhaseType *
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::GetTopographicPhaseInput() const
+ {
+ if (this->GetNumberOfInputs()<3)
+ {
+ return 0;
+ }
+ return static_cast<const ImagePhaseType *>(this->itk::ProcessObject::GetInput(2));
+ }
+
+ /**
+ * Method GenerateOutputInformaton()
+ **/
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ void
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::GenerateOutputInformation()
+ {
+ // Call superclass implementation
+ Superclass::GenerateOutputInformation();
+
+ // Get pointers to the inputs and output
+ ImageSARConstPointer masterPtr = this->GetMasterInput();
+ ImageSARConstPointer slavePtr = this->GetSlaveInput();
+ ImageOutPointer outputPtr = this->GetOutput();
+
+ // KeyWordList
+ ImageKeywordlist masterKWL = masterPtr->GetImageKeywordlist();
+
+ // Master SAR Dimensions
+ m_nbLinesSAR = this->GetMasterInput()->GetLargestPossibleRegion().GetSize()[1];
+ m_nbColSAR = this->GetMasterInput()->GetLargestPossibleRegion().GetSize()[0];
+
+ /////////////////////// Main Output : conjugate multiplication (into SAR geometry) ///////////////////////
+ // Vector Image :
+ // At Least 3 Components :
+ // _ real part
+ // _ imag part
+ // _ modulus
+ outputPtr->SetNumberOfComponentsPerPixel(3);
+
+ // The output is defined with the Master SAR Image
+ // Origin, Spacing and Size (SAR master geometry)
+ ImageOutPointType outOrigin;
+ outOrigin = masterPtr->GetOrigin();
+ ImageOutSpacingType outSP;
+ outSP = masterPtr->GetSpacing();
+
+ // Define Output Largest Region
+ ImageOutRegionType outputLargestPossibleRegion = masterPtr->GetLargestPossibleRegion();
+
+ outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
+ outputPtr->SetOrigin(outOrigin);
+ outputPtr->SetSpacing(outSP);
+
+ // Add ML factors and bands meaning into keyWordList
+ ImageKeywordlist outputKWL = masterKWL;
+ outputKWL.AddKey("support_data.band.Amplitude", std::to_string(0));
+ outputKWL.AddKey("support_data.band.Phase", std::to_string(1));
+ outputKWL.AddKey("support_data.band.Coherency", std::to_string(2));
+
+ // Set new keyword list to output image
+ outputPtr->SetImageKeywordList(outputKWL);
+
+ ///////// Checks (with input topographic phase keywordlists/metadata) /////////////
+ // Check ML Factors (must be 1 to have the same geometry than SAR images)
+ if (this->GetTopographicPhaseInput() != nullptr)
+ {
+ // Get ImagKeyWordList
+ ImageKeywordlist topoKWL = this->GetTopographicPhaseInput()->GetImageKeywordlist();
+
+ if (topoKWL.HasKey("support_data.ml_ran") && topoKWL.HasKey("support_data.ml_azi"))
+ {
+ // Get Master ML Factors
+ unsigned int topo_MLRan = atoi(topoKWL.GetMetadataByKey("support_data.ml_ran").c_str());
+ unsigned int topo_MLAzi = atoi(topoKWL.GetMetadataByKey("support_data.ml_azi").c_str());
+
+ if ((topo_MLRan != 1) || (topo_MLAzi != 1))
+ {
+ itkExceptionMacro(<<"Error, ML Factor for topographic phase are different than 1.");
+ }
+ }
+ }
+ }
+
+
+
+ /**
+ * Method GenerateInputRequestedRegion
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ void
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::GenerateInputRequestedRegion()
+ {
+ // call the superclass' implementation of this method
+ Superclass::GenerateInputRequestedRegion();
+
+ // Get Output requested region
+ ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
+
+ // Input Region = Output region
+ ImageSARPointer masterPtr = const_cast< ImageSARType * >( this->GetMasterInput() );
+ ImageSARPointer slavePtr = const_cast< ImageSARType * >( this->GetSlaveInput() );
+
+ masterPtr->SetRequestedRegion(outputRequestedRegion);
+ slavePtr->SetRequestedRegion(outputRequestedRegion);
+
+ ///////////// Find the region into topographic phase image (if needed) /////////////
+ if (this->GetTopographicPhaseInput() != nullptr)
+ {
+ ImagePhasePointer phasePtr = const_cast< ImagePhaseType * >( this->GetTopographicPhaseInput() );
+ phasePtr->SetRequestedRegion(outputRequestedRegion);
+ }
+ }
+
+
+ /**
+ * Method ThreadedGenerateData
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ void
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::ThreadedGenerateData(const ImageOutRegionType & outputRegionForThread,
+ itk::ThreadIdType threadId)
+ {
+ if (this->GetTopographicPhaseInput() == nullptr)
+ {
+ // Process for line
+ this->ThreadedGenerateDataWithoutTopographicPhase(outputRegionForThread,threadId);
+ }
+ else
+ {
+ // Process for column
+ this->ThreadedGenerateDataWithTopographicPhase(outputRegionForThread,threadId);
+ }
+ }
+
+ /**
+ * Method ThreadedGenerateDataWithoutTopographicPhase
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ void
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::ThreadedGenerateDataWithoutTopographicPhase(const ImageOutRegionType & outputRegionForThread,
+ itk::ThreadIdType threadId)
+ {
+ // Define/declare an iterator that will walk the input regions for this
+ // thread.
+ InputSARIterator inMasterIt(this->GetMasterInput(), outputRegionForThread);
+ InputSARIterator inSlaveIt(this->GetSlaveInput(), outputRegionForThread);
+
+ // Define/declare an iterator that will walk the output region for this
+ // thread.
+ OutputIterator outIt(this->GetOutput(), outputRegionForThread);
+
+ // Support progress methods/callbacks
+ itk::ProgressReporter progress( this, threadId, outputRegionForThread.GetNumberOfPixels() );
+
+ inMasterIt.GoToBegin();
+ inSlaveIt.GoToBegin();
+ outIt.GoToBegin();
+
+ // Output Pixel (VectorImage Pixel)
+ ImageOutPixelType outPixel;
+ outPixel.Reserve(3);
+
+ // For each Line
+ while ( !inMasterIt.IsAtEnd() && !inSlaveIt.IsAtEnd() && !outIt.IsAtEnd())
+ {
+ inMasterIt.GoToBeginOfLine();
+ inSlaveIt.GoToBeginOfLine();
+ outIt.GoToBeginOfLine();
+
+ // For each column
+ while (!inMasterIt.IsAtEndOfLine() && !inSlaveIt.IsAtEndOfLine() && !outIt.IsAtEndOfLine())
+ {
+ // Conjugate multiplication (master * conj(slave))
+ double real_interfero = (inMasterIt.Get().real()*inSlaveIt.Get().real() +
+ inMasterIt.Get().imag()*inSlaveIt.Get().imag());
+
+ double imag_interfero = (inMasterIt.Get().imag()*inSlaveIt.Get().real() -
+ inMasterIt.Get().real()*inSlaveIt.Get().imag());
+
+ ///////////// Output assignations ///////////////
+ outPixel[0] = real_interfero;
+
+
+ outPixel[1] = imag_interfero;
+
+
+ outPixel[2] = sqrt((real_interfero*real_interfero) +
+ (imag_interfero*imag_interfero));
+
+
+ outIt.Set(outPixel);
+
+ // Next colunm inputs
+ ++inMasterIt;
+ ++inSlaveIt;
+ // Next colunm output
+ ++outIt;
+ }
+
+ // Next line intputs
+ inMasterIt.NextLine();
+ inSlaveIt.NextLine();
+ // Next line output
+ outIt.NextLine();
+ }
+ }
+
+ /**
+ * Method ThreadedGenerateDataWithtTopographicPhase
+ */
+ template<class TImageSAR, class TImagePhase, class TImageOut>
+ void
+ SARCompensatedComplexImageFilter< TImageSAR, TImagePhase, TImageOut >
+ ::ThreadedGenerateDataWithTopographicPhase(const ImageOutRegionType & outputRegionForThread,
+ itk::ThreadIdType threadId)
+ {
+ // Define/declare an iterator that will walk the input regions for this
+ // thread.
+ InputSARIterator inMasterIt(this->GetMasterInput(), outputRegionForThread);
+ InputSARIterator inSlaveIt(this->GetSlaveInput(), outputRegionForThread);
+
+ // Define/declare an iterator that will walk the output region for this
+ // thread.
+ OutputIterator outIt(this->GetOutput(), outputRegionForThread);
+
+ // Topographic phase
+ InputPhaseIterator inTopoPhaseIt(this->GetTopographicPhaseInput(), outputRegionForThread);
+
+ // Support progress methods/callbacks
+ itk::ProgressReporter progress( this, threadId, outputRegionForThread.GetNumberOfPixels() );
+
+ inMasterIt.GoToBegin();
+ inSlaveIt.GoToBegin();
+ inTopoPhaseIt.GoToBegin();
+ outIt.GoToBegin();
+
+ // Output Pixel (VectorImage Pixel)
+ ImageOutPixelType outPixel;
+ outPixel.Reserve(3);
+
+ // For each Line
+ while ( !inMasterIt.IsAtEnd() && !inSlaveIt.IsAtEnd() && !outIt.IsAtEnd())
+ {
+ inMasterIt.GoToBeginOfLine();
+ inSlaveIt.GoToBeginOfLine();
+ inTopoPhaseIt.GoToBeginOfLine();
+ outIt.GoToBeginOfLine();
+
+ // For each column
+ while (!inMasterIt.IsAtEndOfLine() && !inSlaveIt.IsAtEndOfLine() && !outIt.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];
+
+ 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);
+
+ ///////////// Output assignations ///////////////
+ outPixel[0] = real_interfero;
+
+
+ outPixel[1] = imag_interfero;
+
+
+ outPixel[2] = sqrt((real_interfero*real_interfero) +
+ (imag_interfero*imag_interfero));
+ }
+ else
+ {
+ outPixel[0] = 0;
+ outPixel[1] = 0;
+ outPixel[2] = 0;
+ outPixel[3] = 0;
+ }
+
+ outIt.Set(outPixel);
+
+ // Next colunm inputs
+ ++inMasterIt;
+ ++inSlaveIt;
+ ++inTopoPhaseIt;
+ // Next colunm output
+ ++outIt;
+ }
+
+ // Next line intputs
+ inMasterIt.NextLine();
+ inSlaveIt.NextLine();
+ inTopoPhaseIt.NextLine();
+ // Next line output
+ outIt.NextLine();
+ }
+ }
+
+
+ } /*namespace otb*/
+
+#endif
diff --git a/include/otbSARGroupedByOrthoImageFilter.h b/include/otbSARGroupedByOrthoImageFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..37f647a1ad98b9b5a028f9aa7687dfe2f18a7c0f
--- /dev/null
+++ b/include/otbSARGroupedByOrthoImageFilter.h
@@ -0,0 +1,230 @@
+/*
+ * Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbSARGroupedByOrthoImageFilter_h
+#define otbSARGroupedByOrthoImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkSmartPointer.h"
+#include "itkPoint.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageScanlineIterator.h"
+
+#include "otbImageKeywordlist.h"
+#include "otbSarSensorModelAdapter.h"
+
+
+#if defined(__GNUC__) || defined(__clang__)
+# pragma GCC diagnostic push
+# pragma GCC diagnostic ignored "-Woverloaded-virtual"
+#if defined(__GNUC__) && (__GNUC__ > 5)
+#pragma GCC diagnostic ignored "-Wnonnull-compare"
+#endif
+# include "ossim/base/ossimGeoidEgm96.h"
+# pragma GCC diagnostic pop
+#else
+# include "ossim/base/ossimGeoidEgm96.h"
+#endif
+
+
+namespace otb
+{
+/** \class SARGroupedByOrthoImageFilter
+ * \brief Creates an interferogram into ground geometry between two images.
+ *
+ * This filter built the Ortho interferogram between a master and a slave image.
+ *
+ * The output is a vector image with amplitude, phase and coherency.
+ *
+ * \ingroup DiapOTBModule
+ */
+
+ template <typename TImageVector, typename TImageDEM, typename TImageOut>
+ class ITK_EXPORT SARGroupedByOrthoImageFilter :
+ public itk::ImageToImageFilter<TImageVector,TImageOut>
+{
+public:
+
+ // Standard class typedefs
+ typedef SARGroupedByOrthoImageFilter Self;
+ typedef itk::ImageToImageFilter<TImageVector,TImageOut> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ // Method for creation through object factory
+ itkNewMacro(Self);
+ // Run-time type information
+ itkTypeMacro(SARGroupedByOrthoImageFilter,ImageToImageFilter);
+
+ /** Typedef to vector image input type : otb::VectorImage */
+ typedef TImageVector ImageVectorType;
+ /** Typedef to describe the output image pointer type. */
+ typedef typename ImageVectorType::Pointer ImageVectorPointer;
+ typedef typename ImageVectorType::ConstPointer ImageVectorConstPointer;
+ /** Typedef to describe the output image region type. */
+ typedef typename ImageVectorType::RegionType ImageVectorRegionType;
+ /** Typedef to describe the type of pixel and point for output image. */
+ typedef typename ImageVectorType::PixelType ImageVectorPixelType;
+ typedef typename ImageVectorType::PointType ImageVectorPointType;
+ /** Typedef to describe the image index, size types and spacing for output image. */
+ typedef typename ImageVectorType::IndexType ImageVectorIndexType;
+ typedef typename ImageVectorType::IndexValueType ImageVectorIndexValueType;
+ typedef typename ImageVectorType::SizeType ImageVectorSizeType;
+ typedef typename ImageVectorType::SizeValueType ImageVectorSizeValueType;
+ typedef typename ImageVectorType::SpacingType ImageVectorSpacingType;
+ typedef typename ImageVectorType::SpacingValueType ImageVectorSpacingValueType;
+
+ /** Typedef to image output type : otb::VectorImage */
+ typedef TImageOut ImageOutType;
+ /** Typedef to describe the output image pointer type. */
+ typedef typename ImageOutType::Pointer ImageOutPointer;
+ typedef typename ImageOutType::ConstPointer ImageOutConstPointer;
+ /** Typedef to describe the output image region type. */
+ typedef typename ImageOutType::RegionType ImageOutRegionType;
+ /** Typedef to describe the type of pixel and point for output image. */
+ typedef typename ImageOutType::PixelType ImageOutPixelType;
+ typedef typename ImageOutType::PointType ImageOutPointType;
+ /** Typedef to describe the image index, size types and spacing for output image. */
+ typedef typename ImageOutType::IndexType ImageOutIndexType;
+ typedef typename ImageOutType::IndexValueType ImageOutIndexValueType;
+ typedef typename ImageOutType::SizeType ImageOutSizeType;
+ typedef typename ImageOutType::SizeValueType ImageOutSizeValueType;
+ typedef typename ImageOutType::SpacingType ImageOutSpacingType;
+ typedef typename ImageOutType::SpacingValueType ImageOutSpacingValueType;
+
+ /** Typedef to dem type : otb::Image */
+ typedef TImageDEM ImageDEMType;
+ /** Typedef to describe the output image pointer type. */
+ typedef typename ImageDEMType::Pointer ImageDEMPointer;
+ typedef typename ImageDEMType::ConstPointer ImageDEMConstPointer;
+ typedef typename ImageDEMType::SpacingType ImageDEMSpacingType;
+ typedef typename ImageDEMType::SpacingValueType ImageDEMSpacingValueType;
+
+ // DEMHandler
+ typedef otb::DEMHandler DEMHandlerType;
+ typedef typename DEMHandlerType::Pointer DEMHandlerPointerType;
+
+
+ // Define Point2DType and Point3DType
+ using Point2DType = itk::Point<double,2>;
+ using Point3DType = itk::Point<double,3>;
+ typedef double ValueType;
+
+ // Typedef for iterators
+ typedef itk::ImageScanlineConstIterator< ImageVectorType > InputIterator;
+ typedef itk::ImageScanlineIterator< ImageOutType > OutputIterator;
+
+
+
+ // Setter/Getter for inputs (Cartesian mean and GroupedByOrtho)
+ /** Connect one of the operands for interferometry : Master */
+ void SetCartesianMeanInput( const ImageVectorType* image);
+
+ /** Connect one of the operands for interferometry : Coregistrated Slave */
+ void SetCompensatedComplexInput(const ImageVectorType* image);
+
+ /** Get the inputs */
+ const ImageVectorType* GetCartesianMeanInput() const;
+ const ImageVectorType* GetCompensatedComplexInput() const;
+
+
+ // Setter for metadata
+ void SetSARImageKeyWorList(ImageKeywordlist sarImageKWL);
+ void SetDEMImagePtr(ImageDEMPointer demPtr);
+
+protected:
+ // Constructor
+ SARGroupedByOrthoImageFilter();
+
+ // Destructor
+ virtual ~SARGroupedByOrthoImageFilter() ITK_OVERRIDE;
+
+ // Print
+ void PrintSelf(std::ostream & os, itk::Indent indent) const ITK_OVERRIDE;
+
+ /** SARGroupedByOrthoImageFilter produces an image which are into output geometry. The differences between
+ * output and input images can be the size of images and the dimensions.
+ * As such, SARGroupedByOrthoImageFilter needs to provide an implementation for
+ * GenerateOutputInformation() in order to inform the pipeline execution model.
+ */
+ virtual void GenerateOutputInformation() ITK_OVERRIDE;
+
+ /** SARGroupedByOrthoImageFilter 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.
+ * As such, TilesAnalysesImageFilter needs to provide an implementation for
+ * GenerateInputRequestedRegion() in order to inform the pipeline execution model.
+ * \sa ProcessObject::GenerateInputRequestedRegion() */
+ virtual void GenerateInputRequestedRegion() ITK_OVERRIDE;
+
+ /**
+ * OutputRegionToInputRegion returns the input SAR region
+ */
+ ImageVectorRegionType OutputRegionToInputRegion(const ImageOutRegionType& outputRegion,
+ bool & intoInputImage) const;
+
+
+ /**
+ * SARGroupedByOrthoImageFilterr can be implemented as a multithreaded filter.
+ * 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"
+ *
+ * \sa ImageToImageFilter::ThreadedGenerateData(),
+ * ImageToImageFilter::GenerateData() */
+ virtual void ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread,
+ itk::ThreadIdType threadId) ITK_OVERRIDE;
+
+
+ private:
+ SARGroupedByOrthoImageFilter(const Self&); // purposely not implemented
+ void operator=(const Self &); // purposely not
+
+ // Instance of SarSensorModelAdapter
+ SarSensorModelAdapter::Pointer m_SarSensorModelAdapter;
+
+ // SAR Image KeyWorldList
+ ImageKeywordlist m_SarImageKwl;
+
+ ossimGeoidEgm96 * m_geoidEmg96;
+
+ // DEM
+ ImageDEMConstPointer m_DEMPtr;
+
+ float m_Gain;
+
+ int m_nbLinesSAR;
+ int m_nbColSAR;
+
+};
+
+} // End namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSARGroupedByOrthoImageFilter.txx"
+#endif
+
+
+
+#endif
diff --git a/include/otbSARGroupedByOrthoImageFilter.txx b/include/otbSARGroupedByOrthoImageFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..bb6105306270dab41d481853fa55f659c82872fd
--- /dev/null
+++ b/include/otbSARGroupedByOrthoImageFilter.txx
@@ -0,0 +1,664 @@
+/*
+ * Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbSARGroupedByOrthoImageFilter_txx
+#define otbSARGroupedByOrthoImageFilter_txx
+
+#include "otbSARGroupedByOrthoImageFilter.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageScanlineIterator.h"
+#include "itkProgressReporter.h"
+#include "itkNumericTraitsPointPixel.h"
+#include "itkContinuousIndex.h"
+
+#include "otbConfigurationManager.h"
+
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+#include "ossim/base/ossimGpt.h"
+#include "ossim/base/ossimEcefPoint.h"
+#pragma GCC diagnostic pop
+#else
+#include "ossim/base/ossimGpt.h"
+#include "ossim/base/ossimEcefPoint.h"
+#endif
+
+#include <cmath>
+#include <algorithm>
+
+
+namespace otb
+{
+ /**
+ * Constructor with default initialization
+ */
+ template <class TImageVector, class TImageDEM, class TImageOut>
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >::SARGroupedByOrthoImageFilter()
+ : m_SarSensorModelAdapter(ITK_NULLPTR), m_geoidEmg96(NULL), m_DEMPtr(ITK_NULLPTR), m_Gain(1.)
+ {
+ // Inputs required
+ this->SetNumberOfRequiredInputs(2);
+
+ std::cout << "ConfigurationManager::GetGeoidFile() = " << ConfigurationManager::GetGeoidFile() << std::endl;
+
+ DEMHandlerPointerType DEMHandler = DEMHandler::Instance();
+
+ m_geoidEmg96 = 0;
+ const std::string isEmg96 = "egm96";
+
+ if (!(ConfigurationManager::GetGeoidFile().empty()))
+ {
+ // DEMHandler instance to specify the geoid file
+ DEMHandler->OpenGeoidFile(ConfigurationManager::GetGeoidFile());
+
+ // If Geoid by default (emg96) instanciate directly a ossimGeoidEgm96 (increase performance)
+ std::size_t found = ConfigurationManager::GetGeoidFile().find(isEmg96);
+ if (found!=std::string::npos)
+ {
+ m_geoidEmg96 = new ossimGeoidEgm96(ossimFilename(ConfigurationManager::GetGeoidFile()));
+ }
+ }
+
+ }
+
+ /**
+ * Destructor
+ */
+ template <class TImageVector, class TImageDEM, class TImageOut>
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >::~SARGroupedByOrthoImageFilter()
+ {
+ if (m_geoidEmg96)
+ {
+ delete m_geoidEmg96;
+ m_geoidEmg96 = 0;
+ }
+ }
+
+ /**
+ * Print
+ */
+ template<class TImageVector, class TImageDEM, class TImageOut>
+ void
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+ ::PrintSelf(std::ostream & os, itk::Indent indent) const
+ {
+ Superclass::PrintSelf(os, indent);
+ }
+
+ /**
+ * Set Master Image
+ */
+ template<class TImageVector, class TImageDEM, class TImageOut>
+ void
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+ ::SetCartesianMeanInput(const ImageVectorType* image )
+ {
+ // Process object is not const-correct so the const casting is required.
+ this->SetNthInput(0, const_cast<ImageVectorType *>(image));
+ }
+
+ /**
+ * Set Coregistrated Slave Image
+ */
+ template<class TImageVector, class TImageDEM, class TImageOut>
+ void
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+ ::SetCompensatedComplexInput(const ImageVectorType* image )
+ {
+ // Process object is not const-correct so the const casting is required.
+ this->SetNthInput(1, const_cast<ImageVectorType *>(image));
+ }
+
+ /**
+ * Get Master Image
+ */
+ template<class TImageVector, class TImageDEM, class TImageOut>
+ const typename SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >::ImageVectorType *
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+ ::GetCartesianMeanInput() const
+ {
+ if (this->GetNumberOfInputs()<1)
+ {
+ return 0;
+ }
+ return static_cast<const ImageVectorType *>(this->itk::ProcessObject::GetInput(0));
+ }
+
+ /**
+ * Get Coregistrated Slave Image
+ */
+ template<class TImageVector, class TImageDEM, class TImageOut>
+ const typename SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >::ImageVectorType *
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+ ::GetCompensatedComplexInput() const
+ {
+ if (this->GetNumberOfInputs()<2)
+ {
+ return 0;
+ }
+ return static_cast<const ImageVectorType *>(this->itk::ProcessObject::GetInput(1));
+ }
+
+/**
+ * Set Sar Image keyWordList
+ */
+template<class TImageVector, class TImageDEM, class TImageOut>
+void
+SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+::SetSARImageKeyWorList(ImageKeywordlist sarImageKWL)
+{
+ m_SarImageKwl = sarImageKWL;
+}
+
+/**
+ * Set DEM Pointer
+ */
+template<class TImageVector, class TImageDEM, class TImageOut>
+void
+SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+::SetDEMImagePtr(ImageDEMPointer demPtr)
+{
+ m_DEMPtr = demPtr;
+}
+
+ /**
+ * Method GenerateOutputInformaton()
+ **/
+ template<class TImageVector, class TImageDEM, class TImageOut>
+ void
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+ ::GenerateOutputInformation()
+ {
+ // Call superclass implementation
+ Superclass::GenerateOutputInformation();
+
+ // Get pointer to output
+ ImageOutPointer outputPtr = this->GetOutput();
+
+ /////////////////////// Main Output : interferogram (into Ground geometry) ///////////////////////
+ // Vector Image :
+ // At Least 3 Components :
+ // _ amplitude
+ // _ phase
+ // _ coherency
+ outputPtr->SetNumberOfComponentsPerPixel(4);
+
+ // The output is defined with the DEM Image
+ // Origin, Spacing and Size (Ground geometry)
+ ImageOutPointType outOrigin;
+ outOrigin = m_DEMPtr->GetOrigin();
+
+ // Define Output Largest Region
+ ImageOutRegionType outputLargestPossibleRegion = m_DEMPtr->GetLargestPossibleRegion();
+ outputLargestPossibleRegion.SetSize(m_DEMPtr->GetLargestPossibleRegion().GetSize());
+ outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
+ outputPtr->SetOrigin(outOrigin);
+ outputPtr->SetSignedSpacing(m_DEMPtr->GetSignedSpacing());
+
+ // Projection Ref
+ outputPtr->SetProjectionRef(m_DEMPtr->GetProjectionRef());
+
+ // Create and Initilaze the SarSensorModelAdapter
+ m_SarSensorModelAdapter = SarSensorModelAdapter::New();
+ bool loadOk = m_SarSensorModelAdapter->LoadState(m_SarImageKwl);
+
+ if(!loadOk || !m_SarSensorModelAdapter->IsValidSensorModel())
+ {
+ itkExceptionMacro(<<"SAR image does not contain a valid SAR sensor model.");
+ }
+
+ // Add Bands meaning into keyWordList
+ ImageKeywordlist outputKWL;
+ outputKWL.AddKey("support_data.band.Amplitude", std::to_string(0));
+ outputKWL.AddKey("support_data.band.Phase", std::to_string(1));
+ outputKWL.AddKey("support_data.band.Coherency", std::to_string(2));
+
+ // Set new keyword list to output image
+ outputPtr->SetImageKeywordList(outputKWL);
+
+ m_nbLinesSAR = this->GetCompensatedComplexInput()->GetLargestPossibleRegion().GetSize()[1];
+ m_nbColSAR = this->GetCompensatedComplexInput()->GetLargestPossibleRegion().GetSize()[0];
+ }
+
+/**
+* Method OutputREgionToInputRegion
+*/
+template<class TImageVector, class TImageDEM, class TImageOut>
+typename SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >::ImageVectorRegionType
+SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+::OutputRegionToInputRegion(const ImageOutRegionType& outputRegion,
+ bool & intoInputImage) const
+{
+ // Compute the input requested region (size and start index)
+ // Use the image transformations to insure an input requested region
+ ImageOutSizeType outputRequestedRegionSize = outputRegion.GetSize();
+ ImageOutIndexType outputRequestedRegionIndex = outputRegion.GetIndex();
+
+ // Here we are breaking traditional pipeline steps because we need to access the input field data
+ // so as to compute the input image requested region
+ // The direct localisation (SAR to DEM) is performed in order to determine the DEM region (input)
+ // represented by the SAR region (output)
+
+ 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;
+
+ // Indice for the SAR bloc (determined by the Pipeline)
+ ImageOutIndexType id[4] ;
+ id[0][0] = outputRequestedRegionIndex[0];
+ id[0][1] = outputRequestedRegionIndex[1];
+ id[1][0] = outputRequestedRegionIndex[0];
+ id[1][1] = outputRequestedRegionIndex[1] + outputRequestedRegionSize[1];
+ id[2][0] = outputRequestedRegionIndex[0] + outputRequestedRegionSize[0];
+ id[2][1] = outputRequestedRegionIndex[1];
+ id[3][0] = outputRequestedRegionIndex[0] + outputRequestedRegionSize[0];
+ id[3][1] = outputRequestedRegionIndex[1] + outputRequestedRegionSize[1];
+
+ Point2DType pixelSAR_Into_MNT_LatLon(0);
+ Point2DType pixelSAR(0);
+ // ImageType::IndexType pixelSAR_Into_MNT;
+ itk::ContinuousIndex<double,2> pixelSAR_Into_MNT;
+
+ Point3DType demGeoPoint;
+ Point2DType demLatLonPoint;
+ Point2DType col_row(0);
+ Point2DType y_z(0);
+ ossimGpt gptPt;
+
+ // Initialie the first and last DEM indice
+ int firstL, firstC, lastL, lastC;
+
+ firstL = nbLinesSAR;
+ lastL = 0;
+ 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
+ for (int i = 0; i < 4; i++)
+ {
+ // Transform index to Lat/Lon Point
+ m_DEMPtr->TransformIndexToPhysicalPoint(id[i], demLatLonPoint);
+ demGeoPoint[0] = demLatLonPoint[0];
+ demGeoPoint[1] = demLatLonPoint[1];
+
+ // Get elevation from earth geoid thanks to DEMHandler or ossimgeoidEmg96
+ double h = 0;
+ if (m_geoidEmg96)
+ {
+ gptPt.lon = demLatLonPoint[0];
+ gptPt.lat = demLatLonPoint[1];
+ h = m_geoidEmg96->offsetFromEllipsoid(gptPt);
+ }
+ //demGeoPoint[2] = h;
+ demGeoPoint[2] = 0;
+
+ // Call the method of sarSensorModelAdapter
+ m_SarSensorModelAdapter->WorldToLineSampleYZ(demGeoPoint, col_row, y_z);
+
+ // Assigne the limits
+ if (firstL > col_row[1])
+ {
+ firstL = col_row[1];
+ }
+
+ if (lastL < col_row[1])
+ {
+ lastL = col_row[1];
+ }
+
+ if (firstC > col_row[0])
+ {
+ firstC = col_row[0];
+ }
+
+ if (lastC < col_row[0])
+ {
+ lastC = col_row[0];
+ }
+
+ }
+
+ // Apply the marge
+ firstL -= margin;
+ firstC -= margin;
+ lastL += margin;
+ lastC += margin;
+
+ // Check the limits
+ if (firstC < 0)
+ {
+ firstC = 0;
+ }
+ if (firstL < 0)
+ {
+ firstL = 0;
+ }
+ if (lastC > nbColSAR-1)
+ {
+ lastC = nbColSAR-1;
+ }
+ if (lastL > nbLinesSAR-1)
+ {
+ lastL = nbLinesSAR-1;
+ }
+
+
+ // Check if into input Image
+
+ if (firstC > nbColSAR-1 || firstL > nbLinesSAR-1)
+ {
+ intoInputImage = false;
+ }
+ if (lastC <= 0 || lastL <= 0)
+ {
+ intoInputImage = false;
+ }
+
+
+ // Transform sides to region
+ ImageVectorIndexType inputRequestedRegionIndex;
+ ImageVectorRegionType inputRequestedRegion = outputRegion;
+ if (intoInputImage)
+ {
+ inputRequestedRegionIndex[0] = static_cast<ImageVectorIndexValueType>(firstC);
+ inputRequestedRegionIndex[1] = static_cast<ImageVectorIndexValueType>(firstL);
+ ImageVectorSizeType inputRequestedRegionSize;
+ inputRequestedRegionSize[0] = static_cast<ImageVectorIndexValueType>(lastC - firstC);
+ inputRequestedRegionSize[1] = static_cast<ImageVectorIndexValueType>(lastL - firstL);
+
+ 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;
+}
+
+ /**
+ * Method GenerateInputRequestedRegion
+ */
+ template<class TImageVector, class TImageDEM, class TImageOut>
+ void
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+ ::GenerateInputRequestedRegion()
+ {
+ // call the superclass' implementation of this method
+ Superclass::GenerateInputRequestedRegion();
+
+ // Get Output requested region
+ ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
+
+ // Input Region with OutputToInputRegion
+ bool isIntoInputImages;
+ ImageVectorRegionType inputRequestedRegion = OutputRegionToInputRegion(outputRequestedRegion,
+ isIntoInputImages);
+ 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);
+ compensatedComplexPtr->SetRequestedRegion(inputRequestedRegion);
+ // }
+ // else
+ // {
+ // itkExceptionMacro(<<"DEM and SAR geometry do not match.");
+ // }
+ }
+
+
+ /**
+ * Method ThreadedGenerateDataWithoutTopographicPhase
+ */
+ template<class TImageVector, class TImageDEM, class TImageOut>
+ void
+ SARGroupedByOrthoImageFilter< TImageVector, TImageDEM, TImageOut >
+ ::ThreadedGenerateData(const ImageOutRegionType & outputRegionForThread,
+ itk::ThreadIdType threadId)
+ {
+
+ // Compute corresponding input region
+ bool isIntoInputImage = true;
+ ImageVectorRegionType inputRegionForThread = OutputRegionToInputRegion(outputRegionForThread, isIntoInputImage);
+
+
+ //std::cout << "Pouet : " << inputRegionForThread << " With " << isIntoInputImage << std::endl;
+
+ // 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];
+
+ Point2DType demLonLatPoint;
+ itk::ContinuousIndex<double,2> pixel_Into_DEM_index;
+ ImageOutIndexType index;
+
+ // 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);
+ }
+
+
+ // 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())
+ {
+ CartMeanIt.GoToBeginOfLine();
+ CompComplexIt.GoToBeginOfLine();
+
+ // For each column
+ while (!CartMeanIt.IsAtEndOfLine() && !CompComplexIt.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]);
+
+ // Conversion into World point
+ ossimGpt worldPoint(cartPoint);
+
+ 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]))
+ {
+ //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];
+ }
+ ++ CartMeanIt;
+ ++ CompComplexIt;
+
+ }
+
+ // Next line
+ CartMeanIt.NextLine();
+ CompComplexIt.NextLine();
+ }
+
+ } // End scan inputs
+
+
+ // Loop on output
+ ImageOutPixelType pixelOutput;
+ pixelOutput.Reserve(4);
+ int index_intoArray = 0;
+
+ // For each line
+ while ( !OutIt.IsAtEnd())
+ {
+ OutIt.GoToBeginOfLine();
+
+ // For each colunm
+ while (!OutIt.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])));
+
+ // Phase
+ pixelOutput[1] = std::atan2(tmpArray_imag[index_intoArray],
+ tmpArray_real[index_intoArray]);
+
+ // 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
+ {
+ pixelOutput[0] = 0;
+ pixelOutput[1] = 0;
+ pixelOutput[2] = 0;
+ pixelOutput[3] = 0;
+ }
+
+ // Assigne Main output (Groun geometry)
+ OutIt.Set(pixelOutput);
+ progress.CompletedPixel();
+
+ // Next colunm
+ ++OutIt;
+ ++index_intoArray;
+ }
+
+ // Next line
+ OutIt.NextLine();
+ }
+
+
+ // 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*/
+
+#endif