diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 4de89d4f93954d2fa444d36e90750c5440a97bc6..220370d9b4f894d0242dcaad87fdcc2bcd298917 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -103,3 +103,8 @@ OTB_CREATE_APPLICATION(NAME SARCorrelationRough
SOURCES otbSARCorrelationRough.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
)
+
+OTB_CREATE_APPLICATION(NAME SARCorrelationPrecise
+ SOURCES otbSARCorrelationPrecise.cxx
+ LINK_LIBRARIES ${${otb-module}_LIBRARIES}
+)
diff --git a/app/otbSARCorrelationPrecise.cxx b/app/otbSARCorrelationPrecise.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..15a0a6fb91df5f8499778695f21b7d35acdc24fd
--- /dev/null
+++ b/app/otbSARCorrelationPrecise.cxx
@@ -0,0 +1,185 @@
+/*
+ * 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 "itkFFTNormalizedCorrelationImageFilter.h"
+#include "itkMinimumMaximumImageCalculator.h"
+
+#include "otbSARTemporalCorrelationGridImageFilter.h"
+
+#include <iostream>
+#include <string>
+#include <fstream>
+
+namespace otb
+{
+namespace Wrapper
+{
+
+class SARCorrelationPrecise : public Application
+{
+public:
+ typedef SARCorrelationPrecise Self;
+ typedef itk::SmartPointer<Self> Pointer;
+
+ itkNewMacro(Self);
+ itkTypeMacro(SARCorrelationPrecise, otb::Wrapper::Application);
+
+ // Filters
+ typedef itk::FFTNormalizedCorrelationImageFilter<FloatImageType, FloatImageType> CorFilterType;
+ typedef itk::MinimumMaximumImageCalculator< FloatImageType> MinMaxCalculatorType;
+ typedef otb::SARTemporalCorrelationGridImageFilter<FloatImageType, FloatVectorImageType> CorGridFilterType;
+
+
+private:
+ void DoInit() override
+ {
+ SetName("SARCorrelationPrecise");
+ SetDescription("Computes SAR correlation shift (into temporal domain).");
+
+ SetDocLongDescription("This application computes correlation shifts between two images : "
+ "shift in range and shift in azimut. "
+ "The inputs of this application are MultiLooked images (real images).");
+
+ //Optional descriptors
+ SetDocLimitations("Only Sentinel 1 products and Cosmo are supported for now.");
+ SetDocAuthors("OTB-Team");
+ SetDocSeeAlso(" ");
+ AddDocTag(Tags::SAR);
+
+ //Parameter declarations
+ AddParameter(ParameterType_InputImage, "inmaster", "Input Master image (real image)");
+ SetParameterDescription("inmaster", "Master Image (real image).");
+
+ AddParameter(ParameterType_InputImage, "inslave", "Input Slave image (real image)");
+ SetParameterDescription("inslave", "Slave Image (real image).");
+
+ AddParameter(ParameterType_Int, "mlran", "MultiLook factor on distance");
+ SetParameterDescription("mlran","MultiLook factor on distance.");
+ SetDefaultParameterInt("mlran", 3);
+ SetMinimumParameterIntValue("mlran", 1);
+ MandatoryOff("mlran");
+
+ AddParameter(ParameterType_Int, "mlazi", "MultiLook factor on azimut");
+ SetParameterDescription("mlazi","MultiLook factor on azimut.");
+ SetDefaultParameterInt("mlazi", 3);
+ SetMinimumParameterIntValue("mlazi", 1);
+ MandatoryOff("mlazi");
+
+ AddParameter(ParameterType_Int, "gridsteprange", "Grid step for range dimension (into SLC/SAR geometry)");
+ SetParameterDescription("gridsteprange","Grid step for range dimension (into SLC/SAR geometry).");
+ SetDefaultParameterInt("gridsteprange", 150);
+ SetMinimumParameterIntValue("gridsteprange", 1);
+ MandatoryOff("gridsteprange");
+
+ AddParameter(ParameterType_Int, "gridstepazimut", "Grid step for azimut dimension (into SLC/SAR geometry)");
+ SetParameterDescription("gridstepazimut","Grid step for azimut dimension (into SLC/SAR geometry).");
+ SetDefaultParameterInt("gridstepazimut", 150);
+ SetMinimumParameterIntValue("gridstepazimut", 1);
+ MandatoryOff("gridstepazimut");
+
+ // Parameter Output
+ AddParameter(ParameterType_OutputImage, "out", "Output Correlation grid (Vector Image)");
+ SetParameterDescription("out","Output Correlation Grid Vector Image (Shift_ran, Shift_azi, Correlation_rate).");
+
+ AddRAMParameter();
+
+ SetDocExampleParameterValue("inmaster","s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002_ML.tiff");
+ SetDocExampleParameterValue("inslave","s1b-s4-slc-vv-20160929t014610-20160929t014634-002277-003d71-002_ML.tiff");
+ SetDocExampleParameterValue("out", "out_CorrelationGrid.tiff");
+ }
+
+ void DoUpdateParameters() override
+ {
+// Nothing to do here : all parameters are independent
+ }
+
+void DoExecute() override
+{
+ // Get numeric parameters
+ int factorML_azi = GetParameterInt("mlazi");
+ int factorML_ran = GetParameterInt("mlran");
+ int grid_step_azi = GetParameterInt("gridstepazimut");
+ int grid_step_ran = GetParameterInt("gridsteprange");
+
+ // Log information
+ otbAppLogINFO(<<"ML Factor on azimut :"<<factorML_azi);
+ otbAppLogINFO(<<"ML Factor on range : "<<factorML_ran);
+ otbAppLogINFO(<<"Grid Step for range : "<<grid_step_ran);
+ otbAppLogINFO(<<"Grid Step for azimut : "<<grid_step_azi);
+
+ // Get master and slave image
+ FloatImageType::Pointer MasterPtr = GetParameterFloatImage("inmaster");
+ FloatImageType::Pointer SlavePtr = GetParameterFloatImage("inslave");
+
+
+ // Correlation Filter
+ CorFilterType::Pointer correlationFilter = CorFilterType::New();
+ m_Ref.push_back(correlationFilter.GetPointer());
+ correlationFilter->SetFixedImage(MasterPtr);
+ correlationFilter->SetMovingImage(SlavePtr);
+ correlationFilter->Update();
+
+ // Min/Max calculator
+ MinMaxCalculatorType::Pointer minMaxFilter = MinMaxCalculatorType::New();
+ minMaxFilter->SetImage(correlationFilter->GetOutput());
+ minMaxFilter->ComputeMaximum();
+
+ float shiftSLC_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2)-
+ minMaxFilter->GetIndexOfMaximum()[0]) * static_cast<float>(factorML_ran);
+ float shiftSLC_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2)
+ -minMaxFilter->GetIndexOfMaximum()[1]) * static_cast<float>(factorML_azi);
+
+ float shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2)-
+ minMaxFilter->GetIndexOfMaximum()[0]);
+ float shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2)
+ -minMaxFilter->GetIndexOfMaximum()[1]);
+
+ // Correlation Grid Filter
+ CorGridFilterType::Pointer filterCorrelationGrid = CorGridFilterType::New();
+ m_Ref.push_back(filterCorrelationGrid.GetPointer());
+
+ // Configure CorrelationGrid Filter
+ filterCorrelationGrid->SetMLran(factorML_ran);
+ filterCorrelationGrid->SetMLazi(factorML_azi);
+ filterCorrelationGrid->SetGridStep(grid_step_ran, grid_step_azi);
+ filterCorrelationGrid->SetRoughShift_ran(shiftML_range);
+ filterCorrelationGrid->SetRoughShift_azi(shiftML_azimut);
+
+ //filterCorrelationGrid->SetRoughShift_ran(42);
+ //filterCorrelationGrid->SetRoughShift_azi(692);
+
+ // Define the main pipeline
+ filterCorrelationGrid->SetMasterInput(MasterPtr);
+ filterCorrelationGrid->SetSlaveInput(SlavePtr);
+ SetParameterOutputImage("out", filterCorrelationGrid->GetOutput());
+}
+ // Vector for filters
+ std::vector<itk::ProcessObject::Pointer> m_Ref;
+};
+
+
+}
+
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::SARCorrelationPrecise)
diff --git a/app/otbSARCorrelationRough.cxx b/app/otbSARCorrelationRough.cxx
index 1011102bb967cb7b0d2fa9a715dfc8819bda1358..644fb89ee3046a13344f4b8e425f620afdbd7e37 100644
--- a/app/otbSARCorrelationRough.cxx
+++ b/app/otbSARCorrelationRough.cxx
@@ -82,15 +82,25 @@ private:
MandatoryOff("mlazi");
// Parameter Output
- AddParameter(ParameterType_Float,"shiftran","rough shift in range dimension (output of this application)");
- SetParameterDescription("shiftran", "rough shift in range.");
- SetParameterRole("shiftran", Role_Output);
- SetDefaultParameterFloat("shiftran", 0);
-
- AddParameter(ParameterType_Float,"shiftazi","rough shift in azimut dimension (output of this application)");
- SetParameterDescription("shiftazi", "rough shift in azimut.");
- SetParameterRole("shiftazi", Role_Output);
- SetDefaultParameterFloat("shiftazi", 0);
+ AddParameter(ParameterType_Float,"shiftranslc","rough shift in range dimension into SLC geometry (output of this application)");
+ SetParameterDescription("shiftranslc", "rough shift in range into SLC geometry.");
+ SetParameterRole("shiftranslc", Role_Output);
+ SetDefaultParameterFloat("shiftranslc", 0);
+
+ AddParameter(ParameterType_Float,"shiftazislc","rough shift in azimut dimension into SLC geometry (output of this application)");
+ SetParameterDescription("shiftazislc", "rough shift in azimut into SLC geometry .");
+ SetParameterRole("shiftazislc", Role_Output);
+ SetDefaultParameterFloat("shiftazislc", 0);
+
+ AddParameter(ParameterType_Float,"shiftranml","rough shift in range dimension into ML geometry (output of this application)");
+ SetParameterDescription("shiftranml", "rough shift in range into ML geometry.");
+ SetParameterRole("shiftranml", Role_Output);
+ SetDefaultParameterFloat("shiftranml", 0);
+
+ AddParameter(ParameterType_Float,"shiftaziml","rough shift in azimut dimension into ML geometry (output of this application)");
+ SetParameterDescription("shiftaziml", "rough shift in azimut into ML geometry .");
+ SetParameterRole("shiftaziml", Role_Output);
+ SetDefaultParameterFloat("shiftaziml", 0);
AddRAMParameter();
@@ -130,14 +140,22 @@ void DoExecute() override
minMaxFilter->SetImage(correlationFilter->GetOutput());
minMaxFilter->ComputeMaximum();
- float shift_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2)-
+ float shiftSLC_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2)-
minMaxFilter->GetIndexOfMaximum()[0]) * static_cast<float>(factorML_ran);
- float shift_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2)
+ float shiftSLC_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2)
-minMaxFilter->GetIndexOfMaximum()[1]) * static_cast<float>(factorML_azi);
+ float shiftML_range = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[0]/2)-
+ minMaxFilter->GetIndexOfMaximum()[0]);
+ float shiftML_azimut = ((correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize()[1]/2)
+ -minMaxFilter->GetIndexOfMaximum()[1]);
+
// Assigne Outputs
- SetParameterFloat("shiftran", shift_range);
- SetParameterFloat("shiftazi", shift_azimut);
+ SetParameterFloat("shiftranslc", shiftSLC_range);
+ SetParameterFloat("shiftazislc", shiftSLC_azimut);
+
+ SetParameterFloat("shiftranml", shiftML_range);
+ SetParameterFloat("shiftaziml", shiftML_azimut);
}
// Vector for filters
std::vector<itk::ProcessObject::Pointer> m_Ref;
diff --git a/include/otbSARTemporalCorrelationGridImageFilter.h b/include/otbSARTemporalCorrelationGridImageFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..f0085955f9b0275bbd3fe18e986da48925cfb9c0
--- /dev/null
+++ b/include/otbSARTemporalCorrelationGridImageFilter.h
@@ -0,0 +1,255 @@
+/*
+ * 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 otbSARTemporalCorrelationGridImageFilter_h
+#define otbSARTemporalCorrelationGridImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkSmartPointer.h"
+#include "itkPoint.h"
+
+#include "otbImageKeywordlist.h"
+#include "otbGenericRSTransform.h"
+#include "itkMinimumMaximumImageCalculator.h"
+#include "itkStatisticsImageFilter.h"
+#include "itkRescaleIntensityImageFilter.h"
+
+#include "itkFFTWCommon.h"
+
+namespace otb
+{
+ /** \class SARTemporalCorrelationGridImageFilter
+ * \brief Estimates a deformation grid between two images (master and slave) for the azimut and range
+ * dimension. A correlation rate is also calculated. The correlation is estimed into the temporal domain
+ *
+ * This filter performs the estimation of deformations between two images (the images are reals).
+ *
+ * \ingroup DiapOTBModule
+ */
+
+ template <typename TImageIn, typename TImageOut> class ITK_EXPORT SARTemporalCorrelationGridImageFilter :
+ public itk::ImageToImageFilter<TImageIn,TImageOut>
+ {
+ public:
+
+ // Standard class typedefs
+ typedef SARTemporalCorrelationGridImageFilter Self;
+ typedef itk::ImageToImageFilter<TImageIn,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(SARCorrelationGridFilter,ImageToImageFilter);
+
+ /** Typedef to image input type : otb::Image for master and slave images */
+ typedef TImageIn ImageInType;
+ /** Typedef to describe the inout image pointer type. */
+ typedef typename ImageInType::Pointer ImageInPointer;
+ typedef typename ImageInType::ConstPointer ImageInConstPointer;
+ /** Typedef to describe the inout image region type. */
+ typedef typename ImageInType::RegionType ImageInRegionType;
+ /** Typedef to describe the type of pixel and point for inout image. */
+ typedef typename ImageInType::PixelType ImageInPixelType;
+ typedef typename ImageInType::PointType ImageInPointType;
+ /** Typedef to describe the image index, size types and spacing for inout image. */
+ typedef typename ImageInType::IndexType ImageInIndexType;
+ typedef typename ImageInType::IndexValueType ImageInIndexValueType;
+ typedef typename ImageInType::SizeType ImageInSizeType;
+ typedef typename ImageInType::SizeValueType ImageInSizeValueType;
+ typedef typename ImageInType::SpacingType ImageInSpacingType;
+ typedef typename ImageInType::SpacingValueType ImageInSpacingValueType;
+
+ /** Typedef to image output type : otb::VectorImage with three components (range deformations,
+ azimut deformations and correlation rate) */
+ 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>;
+
+
+ // Define Filter used here (ie : RSTRansform from Master to Slave)
+ typedef typename itk::MinimumMaximumImageCalculator< ImageInType> MinMaxCalculatorType;
+ typedef typename itk::RescaleIntensityImageFilter< ImageInType, ImageInType > RescaleFilterType;
+
+ // ITK proxy to the fftw library
+ typedef typename itk::fftw::Proxy<ImageInPixelType> FFTWProxyType;
+ typedef typename FFTWProxyType::ComplexType FFTProxyComplexType;
+ typedef typename FFTWProxyType::PixelType FFTProxyPixelType;
+ typedef typename FFTWProxyType::PlanType FFTProxyPlanType;
+
+ // Setter/Getter
+ itkSetMacro(RoughShift_ran, int);
+ itkGetMacro(RoughShift_ran, int);
+ itkSetMacro(RoughShift_azi, int);
+ itkGetMacro(RoughShift_azi, int);
+
+ /** Set/Get fine coherency threshold */
+ itkSetMacro(CorrelationThreshold, double);
+ itkGetMacro(CorrelationThreshold, double);
+
+ /** Set/Get unsigned int grid step */
+ void SetGridStep(unsigned int stepRange, unsigned int stepAzimut)
+ {
+ m_GridStep[0] = stepRange;
+ m_GridStep[1] = stepAzimut;
+ }
+ unsigned int GetGridStepRange()
+ {
+ return m_GridStep[0];
+ }
+ unsigned int GetGridStepAzimut()
+ {
+ return m_GridStep[1];
+ }
+ /** Set/Get ML factors */
+ itkSetMacro(MLran, unsigned int);
+ itkGetMacro(MLran, unsigned int);
+ itkSetMacro(MLazi, unsigned int);
+ itkGetMacro(MLazi, unsigned int);
+
+
+ // Setter/Getter for master and slave images (inputs)
+ /** Connect one of the operands for registration */
+ void SetMasterInput( const ImageInType* image);
+
+ /** Connect one of the operands for registration */
+ void SetSlaveInput(const ImageInType* image);
+
+ /** Get the inputs */
+ const ImageInType * GetMasterInput() const;
+ const ImageInType * GetSlaveInput() const;
+
+ protected:
+ // Constructor
+ SARTemporalCorrelationGridImageFilter();
+
+ // Destructor
+ virtual ~SARTemporalCorrelationGridImageFilter() ITK_OVERRIDE {};
+
+ // Print
+ void PrintSelf(std::ostream & os, itk::Indent indent) const ITK_OVERRIDE;
+
+ /** SARTemporalCorrelationGridImageFilter produces an vector image to indicate the deformation for the two dimensions
+ * and the correlation rate. The differences between output and input images are the size of images and the
+ * dimensions. The input images a classic images and the output is a otb::VectorImage with three components.
+ * As such, SARTemporalCorrelationGridImageFilter needs to provide an implementation for
+ * GenerateOutputInformation() in order to inform the pipeline execution model.
+ */
+ virtual void GenerateOutputInformation() ITK_OVERRIDE;
+
+ /** SARTemporalCorrelationGridImageFilter needs input requested regions (for master and slave images) that
+ * corresponds to the projection into the requested region of the deformation grid (our output requested
+ * region).
+ * As such, SARQuadraticAveragingImageFilter needs to provide an implementation for
+ * GenerateInputRequestedRegion() in order to inform the pipeline execution model.
+ * \sa ProcessObject::GenerateInputRequestedRegion() */
+ virtual void GenerateInputRequestedRegion() ITK_OVERRIDE;
+
+ /**
+ * OutputRegionToInputRegion assigne master and slave regions
+ */
+ void OutputRegionToInputRegion(const ImageOutRegionType& outputRegion) const;
+
+ /** SARTemporalCorrelationGridImageFilter can be implemented as a multithreaded filter.
+ * Therefore, this implementation provides a ThreadedGenerateData() routine
+ * which is called for each processing thread. The 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() */
+ void ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread,
+ itk::ThreadIdType threadId ) override;
+
+
+ private:
+ SARTemporalCorrelationGridImageFilter(const Self&); // purposely not implemented
+ void operator=(const Self &); // purposely not
+
+ /** Function to interpolate the best interger shift */
+ void PIC(float TAB_PIC[3][3], float *DX, float *DY, float *VAL_MAX, int *CR) const;
+
+ inline double Correlation(const std::vector<double>& master,
+ const std::vector<double>& slave) const;
+
+ inline double CorrelationDiapason(const std::vector<double>& master,
+ const std::vector<double>& slave) const;
+
+ /** Coarse correlation threshold */
+ double m_CorrelationThreshold;
+
+ /** Grid step */
+ ImageInSizeType m_GridStep;
+
+ /** Grid step into ML geometry */
+ ImageInSizeType m_GridStep_ML;
+
+ /** ML factors (in range and in azimut) */
+ unsigned int m_MLran;
+ unsigned int m_MLazi;
+
+ /** Input rough shifts */
+ int m_RoughShift_ran;
+ int m_RoughShift_azi;
+
+ /** Search Window around rough shifts (into ML geometry) */
+ int m_WinAround_ShiftRan;
+ int m_WinAround_ShiftAzi;
+
+ /** Window size to estimate correlation */
+ int m_WinCor_ran;
+ int m_WinCor_azi;
+
+ // Extract Parameters
+ ImageInPointer * m_masterExtractPieceTab;
+ ImageInPointer * m_slaveExtractPieceTab;
+
+ };
+
+} // End namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSARTemporalCorrelationGridImageFilter.txx"
+#endif
+
+
+
+#endif
diff --git a/include/otbSARTemporalCorrelationGridImageFilter.txx b/include/otbSARTemporalCorrelationGridImageFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..d51e18eda9176506f8526ee8ba5a70b2dc14b1e0
--- /dev/null
+++ b/include/otbSARTemporalCorrelationGridImageFilter.txx
@@ -0,0 +1,797 @@
+/*
+ * 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 otbSARTemporalCorrelationGridImageFilter_txx
+#define otbSARTemporalCorrelationGridImageFilter_txx
+
+#include "otbSARTemporalCorrelationGridImageFilter.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkImageScanlineIterator.h"
+#include "itkProgressReporter.h"
+#include "itkNumericTraitsPointPixel.h"
+#include "itkContinuousIndex.h"
+
+#include "otbDEMHandler.h"
+
+#include "itkFFTWCommon.h"
+#include "itkFFTWGlobalConfiguration.h"
+
+#include <cmath>
+
+
+namespace otb
+{
+ /**
+ * Constructor with default initialization
+ */
+ template <class TImageIn, class TImageOut>
+ SARTemporalCorrelationGridImageFilter< TImageIn ,TImageOut >::SARTemporalCorrelationGridImageFilter()
+ {
+ this->SetNumberOfRequiredInputs(2);
+
+ // Default window sizes
+ m_WinAround_ShiftRan = 10;
+ m_WinAround_ShiftAzi = 20;
+ m_WinCor_ran = 10;
+ m_WinCor_azi = 10;
+
+ // Default rough shits (0)
+ m_RoughShift_ran = 0;
+ m_RoughShift_azi = 0;
+
+ // Grid Step
+ m_GridStep.Fill(1);
+ m_GridStep_ML.Fill(1);
+
+ m_CorrelationThreshold = 0.3;
+
+ m_MLran = 3;
+ m_MLazi = 3;
+ }
+
+ /**
+ * Set Master Image
+ */
+ template <class TImageIn, class TImageOut>
+ void
+ SARTemporalCorrelationGridImageFilter< TImageIn ,TImageOut >
+ ::SetMasterInput(const ImageInType* image )
+ {
+ // Process object is not const-correct so the const casting is required.
+ this->SetNthInput(0, const_cast<ImageInType *>(image));
+ }
+
+ /**
+ * Set Slave Image
+ */
+ template <class TImageIn, class TImageOut>
+ void
+ SARTemporalCorrelationGridImageFilter< TImageIn ,TImageOut >
+ ::SetSlaveInput(const ImageInType* image)
+ {
+ // Process object is not const-correct so the const casting is required.
+ this->SetNthInput(1, const_cast<ImageInType *>(image));
+ }
+
+ /**
+ * Get Master Image
+ */
+ template <class TImageIn, class TImageOut>
+ const typename SARTemporalCorrelationGridImageFilter< TImageIn ,TImageOut>::ImageInType *
+ SARTemporalCorrelationGridImageFilter< TImageIn ,TImageOut >
+ ::GetMasterInput() const
+ {
+ if (this->GetNumberOfInputs()<1)
+ {
+ return 0;
+ }
+ return static_cast<const ImageInType *>(this->itk::ProcessObject::GetInput(0));
+ }
+
+ /**
+ * Get Slave Image
+ */
+ template <class TImageIn, class TImageOut>
+ const typename SARTemporalCorrelationGridImageFilter< TImageIn ,TImageOut>::ImageInType *
+ SARTemporalCorrelationGridImageFilter< TImageIn ,TImageOut >
+ ::GetSlaveInput() const
+ {
+ if (this->GetNumberOfInputs()<2)
+ {
+ return 0;
+ }
+ return static_cast<const ImageInType *>(this->itk::ProcessObject::GetInput(1));
+ }
+
+ /**
+ * Print
+ */
+ template<class TImageIn, class TImageOut>
+ void
+ SARTemporalCorrelationGridImageFilter< TImageIn, TImageOut >
+ ::PrintSelf(std::ostream & os, itk::Indent indent) const
+ {
+ Superclass::PrintSelf(os, indent);
+
+ os << indent << "ML range factor : " << m_MLran << std::endl;
+ os << indent << "ML azimut factor : " << m_MLazi << std::endl;
+ os << indent << "Step for the grid : " << m_GridStep[0] << ", " << m_GridStep[1] << std::endl;
+ }
+
+ /**
+ * Method GenerateOutputInformaton()
+ **/
+ template<class TImageIn, class TImageOut>
+ void
+ SARTemporalCorrelationGridImageFilter< TImageIn, TImageOut >
+ ::GenerateOutputInformation()
+ {
+ // Call superclass implementation
+ Superclass::GenerateOutputInformation();
+
+ // Get pointers to the input and output
+ ImageInType * masterPtr = const_cast<ImageInType * >(this->GetMasterInput());
+ ImageInType * slavePtr = const_cast<ImageInType * >(this->GetSlaveInput());
+ ImageOutPointer outputPtr = this->GetOutput();
+
+ // Check GridSteps (must be a multiple of ML Factors)
+ if (m_GridStep[0] % m_MLran)
+ {
+ itkExceptionMacro(<<"GridSteps range mot a multiple of MLRan.");
+ }
+ if (m_GridStep[1] % m_MLazi)
+ {
+ itkExceptionMacro(<<"GridSteps azimut mot a multiple of MLAzi.");
+ }
+
+ m_GridStep_ML[0] = m_GridStep[0]/m_MLran;
+ m_GridStep_ML[1] = m_GridStep[1]/m_MLazi;
+
+ ///////// Checks (with input keywordlists/metadata) /////////////
+ // Check ML Factors (inputs of this filter with Master and Slave metadata)
+ // Get Master and Slave Keywordlist
+ ImageKeywordlist masterKWL = masterPtr->GetImageKeywordlist();
+ ImageKeywordlist slaveKWL = slavePtr->GetImageKeywordlist();
+
+ if (masterKWL.HasKey("support_data.ml_ran") && masterKWL.HasKey("support_data.ml_azi"))
+ {
+ // Get Master ML Factors
+ unsigned int master_MLRan = atoi(masterKWL.GetMetadataByKey("support_data.ml_ran").c_str());
+ unsigned int master_MLAzi = atoi(masterKWL.GetMetadataByKey("support_data.ml_azi").c_str());
+
+ if (slaveKWL.HasKey("support_data.ml_ran") && slaveKWL.HasKey("support_data.ml_azi"))
+ {
+ // Get Slave ML Factors
+ unsigned int slave_MLRan = atoi(slaveKWL.GetMetadataByKey("support_data.ml_ran").c_str());
+ unsigned int slave_MLAzi = atoi(slaveKWL.GetMetadataByKey("support_data.ml_azi").c_str());
+
+ if ((slave_MLRan != master_MLRan) || (slave_MLAzi != master_MLAzi))
+ {
+ itkExceptionMacro(<<"ML Factor betwwen master and slave are different.");
+ }
+ }
+
+ if ((master_MLRan != m_MLran) || (master_MLAzi != m_MLazi))
+ {
+ itkExceptionMacro(<<"ML Factor betwwen master and inputs of this application are different.");
+ }
+ }
+ else
+ {
+ if (slaveKWL.HasKey("support_data.ml_ran") && slaveKWL.HasKey("support_data.ml_azi"))
+ {
+ // Get Slave ML Factors
+ unsigned int slave_MLRan = atoi(slaveKWL.GetMetadataByKey("support_data.ml_ran").c_str());
+ unsigned int slave_MLAzi = atoi(slaveKWL.GetMetadataByKey("support_data.ml_azi").c_str());
+
+ if ((slave_MLRan != m_MLran) || (slave_MLAzi != m_MLazi))
+ {
+ itkExceptionMacro(<<"ML Factor betwwen slave and inputs of this application are different.");
+ }
+ }
+ }
+
+ // Set the number of Components for the Output VectorImage
+ // 3 Components :
+ // _ range shift between Master and Slave
+ // _ azimut shift between Master and Slave
+ // _ correlation rate
+ outputPtr->SetNumberOfComponentsPerPixel(3);
+
+ // Update size and spacing according to grid step
+ ImageOutRegionType largestRegion = static_cast<ImageOutRegionType>(masterPtr->GetLargestPossibleRegion());
+ ImageOutSizeType outputSize = static_cast<ImageOutSizeType>(largestRegion.GetSize());
+ ImageOutSpacingType outputSpacing = static_cast<ImageOutSpacingType>(masterPtr->GetSpacing());
+ ImageOutPointType outputOrigin = static_cast<ImageOutPointType>(masterPtr->GetOrigin());
+
+ // First/Last Colunm and row
+ int firstC = std::max (1 + m_WinCor_ran, 1 + m_WinCor_ran - m_RoughShift_ran + m_WinAround_ShiftRan);
+ int firstL = std::max (1 + m_WinCor_azi, 1 + m_WinCor_azi - m_RoughShift_azi + m_WinAround_ShiftAzi);
+ int lastC = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[0] - m_RoughShift_ran -
+ m_WinAround_ShiftRan - m_WinCor_ran, masterPtr->GetLargestPossibleRegion().GetSize()[0]
+ - m_WinCor_ran);
+ int lastL = std::min(slavePtr->GetLargestPossibleRegion().GetSize()[1] - m_RoughShift_azi -
+ m_WinAround_ShiftAzi- m_WinCor_azi, masterPtr->GetLargestPossibleRegion().GetSize()[1]
+ - m_WinCor_azi);
+
+ outputSize[0]= lastC - firstC + 1;
+ outputSize[1] = lastL - firstL + 1;
+
+
+ for(unsigned int dim = 0; dim < ImageOutType::ImageDimension; ++dim)
+ {
+ outputSize[dim] = (outputSize[dim] - outputSize[dim]%m_GridStep_ML[dim]) / m_GridStep_ML[dim];
+ outputSpacing[dim] *= m_GridStep_ML[dim];
+ }
+
+ // Set spacing
+ outputPtr->SetSpacing(outputSpacing);
+
+ // Set origin
+ outputPtr->SetOrigin(outputOrigin);
+
+ // Set largest region size
+ largestRegion.SetSize(outputSize);
+ outputPtr->SetLargestPossibleRegion(largestRegion);
+
+ // Add GridSteps into KeyWordList
+ ImageKeywordlist outputKWL;
+ outputKWL.AddKey("support_data.gridstep.range", std::to_string(m_GridStep[0]));
+ outputKWL.AddKey("support_data.gridstep.azimut", std::to_string(m_GridStep[1]));
+
+ outputPtr->SetImageKeywordList(outputKWL);
+ }
+
+ /**
+ * Method OutputRegionToInputRegion for GenerateInputRequestedRegion
+ */
+ template<class TImageIn, class TImageOut>
+ void
+ SARTemporalCorrelationGridImageFilter< TImageIn, TImageOut >
+ ::OutputRegionToInputRegion(const ImageOutRegionType& outputRegion) const
+ {
+ // Get pointers to inputs
+ ImageInType * masterPtr = const_cast<ImageInType * >(GetMasterInput());
+ ImageInType * slavePtr = const_cast<ImageInType * >(GetSlaveInput());
+
+ int nbColSlave = slavePtr->GetLargestPossibleRegion().GetSize()[0];
+ int nbLinesSlave = slavePtr->GetLargestPossibleRegion().GetSize()[1];
+
+ /////////////////////// Requested region for master ////////////////////
+ // Correlation grid is bonded to the master image => should represent with the
+ // m_Grid_step the output requested region
+ // get a copy of the master requested region (should equal the output
+ // requested region)
+ ImageInRegionType masterRequestedRegion, slaveRequestedRegion;
+ masterRequestedRegion = outputRegion;
+
+ // Apply grid step
+ ImageInSizeType masterRequestedSize = masterRequestedRegion.GetSize();
+ ImageInIndexType masterRequestedIndex = masterRequestedRegion.GetIndex();
+
+ for(unsigned int dim = 0; dim < ImageOutType::ImageDimension; ++dim)
+ {
+ masterRequestedSize [dim] *= m_GridStep_ML[dim];
+ masterRequestedIndex[dim] *= m_GridStep_ML[dim];
+ }
+
+ int firstC_Grid = std::max (1 + m_WinCor_ran, 1 + m_WinCor_ran - m_RoughShift_ran + m_WinAround_ShiftRan);
+ int firstL_Grid = std::max (1 + m_WinCor_azi, 1 + m_WinCor_azi - m_RoughShift_azi + m_WinAround_ShiftAzi);
+
+ // With correlation window and first row/colunm
+ masterRequestedIndex[0] -= (m_WinCor_ran - firstC_Grid);
+ masterRequestedIndex[1] -= (m_WinCor_azi - firstL_Grid);
+ masterRequestedSize[0] += (2*m_WinCor_ran+1) + firstC_Grid;
+ masterRequestedSize[1] += (2*m_WinCor_azi+1) + firstL_Grid;
+
+ masterRequestedRegion.SetSize(masterRequestedSize);
+ masterRequestedRegion.SetIndex(masterRequestedIndex);
+
+ masterPtr->SetRequestedRegion(masterRequestedRegion);
+
+ // Crop the fixed region at the fixed's largest possible region (or buffered)
+ if ( masterRequestedRegion.Crop(masterPtr->GetLargestPossibleRegion()))
+ {
+ masterPtr->SetRequestedRegion( masterRequestedRegion );
+ }
+ else
+ {
+ // Couldn't crop the region (requested region is outside the largest
+ // possible region). Throw an exception.
+ // store what we tried to request (prior to trying to crop)
+ masterPtr->SetRequestedRegion( masterRequestedRegion );
+
+ // Build an exception
+ itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
+ std::ostringstream msg;
+ msg << this->GetNameOfClass()
+ << "::GenerateInputRequestedRegion()";
+ e.SetLocation(msg.str().c_str());
+ e.SetDescription("Requested region is (at least partially) outside the largest possible region of image 1.");
+ e.SetDataObject(masterPtr);
+ throw e;
+ }
+
+ //////////////////// Requested region for slave ////////////////////
+ int firstL, firstC, lastL, lastC;
+ // Transform back into slave region index space with an maximum shift between master and slave
+ ImageInIndexType slaveIndex;
+
+ // Find requested region
+ ImageInSizeType slaveSize;
+
+ // Set the max_shift (margin)
+ firstL = masterRequestedIndex[1] - m_RoughShift_azi - m_WinAround_ShiftAzi - m_WinCor_azi;
+ firstC = masterRequestedIndex[0] - m_RoughShift_ran - m_WinAround_ShiftRan - m_WinCor_ran;
+ lastL = masterRequestedIndex[1] + masterRequestedSize[1] + m_RoughShift_azi + m_WinAround_ShiftAzi +
+ m_WinCor_azi;
+ lastC = masterRequestedIndex[0] + masterRequestedSize[0] + m_RoughShift_ran + m_WinAround_ShiftRan + m_WinCor_ran;
+
+ // Check the limits
+ if (firstC < 0)
+ {
+ firstC = 0;
+ }
+ if (firstL < 0)
+ {
+ firstL = 0;
+ }
+ if (lastC > nbColSlave-1)
+ {
+ lastC = nbColSlave-1;
+ }
+ if (lastL > nbLinesSlave-1)
+ {
+ lastL = nbLinesSlave-1;
+ }
+
+
+ slaveIndex[0] = static_cast<ImageInIndexValueType>(firstC);
+ slaveIndex[1] = static_cast<ImageInIndexValueType>(firstL);
+ slaveSize[0] = static_cast<ImageInIndexValueType>(lastC - firstC)+1;
+ slaveSize[1] = static_cast<ImageInIndexValueType>(lastL - firstL)+1;
+
+ slaveRequestedRegion.SetIndex(slaveIndex);
+ slaveRequestedRegion.SetSize(slaveSize);
+
+ // crop the moving region at the moving's largest possible region
+ if ( slaveRequestedRegion.Crop(slavePtr->GetLargestPossibleRegion()))
+ {
+ slavePtr->SetRequestedRegion( slaveRequestedRegion );
+ }
+ else
+ {
+ // Couldn't crop the region (requested region is outside the largest
+ // possible region). This case might happen so we do not throw any exception but
+ // request a null region instead
+ slaveSize.Fill(0);
+ slaveRequestedRegion.SetSize(slaveSize);
+ slaveIndex.Fill(0);
+ slaveRequestedRegion.SetIndex(slaveIndex);
+
+ // store what we tried to request (prior to trying to crop)
+ slavePtr->SetRequestedRegion(slaveRequestedRegion);
+ }
+ }
+
+
+ /**
+ * Method GenerateInputRequestedRegion
+ */
+ template<class TImageIn, class TImageOut>
+ void
+ SARTemporalCorrelationGridImageFilter< TImageIn, TImageOut >
+ ::GenerateInputRequestedRegion()
+ {
+ // Call the superclass' implementation of this method
+ Superclass::GenerateInputRequestedRegion();
+
+ ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
+ this->OutputRegionToInputRegion(outputRequestedRegion);
+ }
+
+
+ /**
+ * Method BeforeThreadedGenerateData
+ */
+ template<class TImageIn, class TImageOut>
+ void
+ SARTemporalCorrelationGridImageFilter< TImageIn, TImageOut >
+ ::ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread,
+ itk::ThreadIdType threadId)
+ {
+ // Get pointers to inputs
+ const ImageInType * masterPtr = this->GetMasterInput();
+ const ImageInType * slavePtr = this->GetSlaveInput();
+
+ // Typedef for iterators
+ typedef itk::ImageScanlineIterator< ImageOutType > OutputIterator;
+ typedef itk::ImageRegionConstIterator<ImageInType> MasterSlaveIt;
+ typedef itk::ImageRegionIterator<ImageInType> ExtractIt;
+
+ // Iterator on output (Grid geometry)
+ OutputIterator OutIt(this->GetOutput(), outputRegionForThread);
+ OutIt.GoToBegin();
+
+ ImageOutPixelType pixelCorrelationGrid;
+ pixelCorrelationGrid.Reserve(3);
+
+ // Support progress methods/callbacks
+ itk::ProgressReporter progress(this, threadId, this->GetOutput()->GetRequestedRegion().GetNumberOfPixels());
+
+ int nbColSlave = slavePtr->GetLargestPossibleRegion().GetSize()[0];
+ int nbLinesSlave = slavePtr->GetLargestPossibleRegion().GetSize()[1];
+
+ // Array to store correlation rate for each pixel and each shift
+ double arrayRate [2*m_WinAround_ShiftAzi+1][2*m_WinAround_ShiftRan+1];
+ float arrayPic[3][3];
+
+ float dx, dy, VAL_MAX;
+ int checkInterpolation;
+
+ // Size for correlation estimation
+ ImageInSizeType masterSize;
+ masterSize[0] = 2*m_WinCor_ran+1;
+ masterSize[1] = 2*m_WinCor_azi+1;
+
+ ImageInSizeType slaveSize;
+ slaveSize[0] = 2*m_WinCor_ran+1;
+ slaveSize[1] = 2*m_WinCor_azi+1;
+
+ // vector master and slave
+ std::vector<double> master;
+ std::vector<double> slave;
+
+ // First L/C for the grid
+ int firstC = std::max (1 + m_WinCor_ran, 1 + m_WinCor_ran - m_RoughShift_ran + m_WinAround_ShiftRan);
+ int firstL = std::max (1 + m_WinCor_azi, 1 + m_WinCor_azi - m_RoughShift_azi + m_WinAround_ShiftAzi);
+
+
+ // For each line of output
+ while ( !OutIt.IsAtEnd())
+ {
+ OutIt.GoToBeginOfLine();
+
+ // For each colunm of the grid
+ while (!OutIt.IsAtEndOfLine())
+ {
+ // Get index
+ ImageInIndexType masterIndex;
+
+ for(unsigned int dim = 0; dim < ImageInType::ImageDimension; ++dim)
+ {
+ masterIndex[dim] = OutIt.GetIndex()[dim] * m_GridStep_ML[dim]; // Grid geometry
+ }
+
+ // To center correlation estimation and shift to the first L/C
+ masterIndex[0] -= (m_WinCor_ran-firstC);
+ masterIndex[1] -= (m_WinCor_azi-firstL);
+
+ // Build the master region (for correlation estimation)
+ ImageInRegionType masterCurrentRegion;
+ masterCurrentRegion.SetIndex(masterIndex);
+ masterCurrentRegion.SetSize(masterSize);
+ masterCurrentRegion.Crop(masterPtr->GetLargestPossibleRegion());
+
+ // Iterators on master Ptr
+ MasterSlaveIt masterIt(masterPtr, masterCurrentRegion);
+ masterIt.GoToBegin();
+
+
+ // Fill master patch
+ master.clear();
+ while( !masterIt.IsAtEnd())
+ {
+ // Add to the master vector
+ double value = masterIt.Get();
+ master.push_back(value);
+
+ ++masterIt;
+ }
+
+ // Best integer shifts and correlation rate
+ int bestShift_ran = 0;
+ int bestShift_azi = 0;
+ double bestRate = 0;
+
+ // Loop on shift around rough ones
+ for (int l = -m_WinAround_ShiftAzi; l <= m_WinAround_ShiftAzi; l++)
+ {
+ for (int c = -m_WinAround_ShiftRan; c <= m_WinAround_ShiftRan; c++)
+ {
+ // Find slave index for current shifts
+ ImageInIndexType slaveIndex;
+ slaveIndex[0] = masterIndex[0] + c + m_RoughShift_ran;
+ slaveIndex[1] = masterIndex[1] + l + m_RoughShift_azi;
+ //slaveIndex[0] = masterIndex[0];
+ //slaveIndex[1] = masterIndex[1];
+
+ // Check Index
+ if (slaveIndex[0] < 0)
+ {
+ slaveIndex[0] = 0;
+ }
+ if (slaveIndex[1] < 0)
+ {
+ slaveIndex[1] = 0;
+ }
+ if (slaveIndex[0] > (nbColSlave-(2*m_WinCor_ran)))
+ {
+ slaveIndex[0] = nbColSlave-2*m_WinCor_ran;
+ }
+ if (slaveIndex[1] > (nbLinesSlave-(2*m_WinCor_azi)))
+ {
+ slaveIndex[1] = nbLinesSlave-2*m_WinCor_azi;
+ }
+
+
+ // Build the slave region (for correlation estimation)
+ ImageInRegionType slaveCurrentRegion;
+ slaveCurrentRegion.SetIndex(slaveIndex);
+ slaveCurrentRegion.SetSize(slaveSize);
+ slaveCurrentRegion.Crop(slavePtr->GetLargestPossibleRegion());
+
+ // Iterators on slave Ptr
+ MasterSlaveIt slaveIt(slavePtr, slaveCurrentRegion);
+ slaveIt.GoToBegin();
+
+ // Fill slave patch
+ slave.clear();
+ while (!slaveIt.IsAtEnd())
+ {
+ // Add to the slave vector
+ double value = slaveIt.Get();
+ slave.push_back(value);
+
+ ++slaveIt;
+ }
+
+ //arrayRate[l][c] = Correlation(master, slave);
+
+ arrayRate[l + m_WinAround_ShiftAzi][c + m_WinAround_ShiftRan] =
+ CorrelationDiapason(master, slave);
+
+ if (arrayRate[l+ m_WinAround_ShiftAzi][c+ m_WinAround_ShiftRan] > bestRate)
+ {
+ bestRate = arrayRate[l + m_WinAround_ShiftAzi][c + m_WinAround_ShiftRan];
+ bestShift_ran = c;
+ bestShift_azi = l;
+ }
+ }
+ }
+
+
+ // Function PIC on the best shifts (interpolation)
+ if ((std::abs(bestShift_ran) < m_WinAround_ShiftRan) &&
+ (std::abs(bestShift_azi) < m_WinAround_ShiftAzi) &&
+ (bestRate >= m_CorrelationThreshold))
+ {
+ // Built array_pic
+ for (int k=-1;k<2;k++)
+ for (int j=-1;j<2;j++)
+ arrayPic[k+1][j+1] = static_cast<float>(arrayRate[bestShift_azi + k + m_WinAround_ShiftAzi]
+ [bestShift_ran + j + m_WinAround_ShiftRan]);
+
+ // Launch interpolation
+ PIC(arrayPic, &dx, &dy, &VAL_MAX, &checkInterpolation);
+
+ pixelCorrelationGrid[0] = m_RoughShift_ran + bestShift_ran + dx;
+ pixelCorrelationGrid[1] = m_RoughShift_azi + bestShift_azi + dy;
+ pixelCorrelationGrid[2] = VAL_MAX;
+
+ //std::cout << arrayPic[1][1] << " at " << bestShift_ran << " , " << bestShift_azi << std::endl;
+
+ }
+ else
+ {
+ pixelCorrelationGrid[0] = m_RoughShift_ran + bestShift_ran;
+ pixelCorrelationGrid[1] = m_RoughShift_azi + bestShift_azi;
+ pixelCorrelationGrid[2] = 0;
+ }
+
+
+ // Correlation into SLC geometry and not ML geometry
+ pixelCorrelationGrid[0] *= m_MLran;
+ pixelCorrelationGrid[1] *= m_MLazi;
+
+ // Assign Output
+ OutIt.Set(pixelCorrelationGrid);
+ progress.CompletedPixel();
+
+ // Next Colunm for output
+ ++OutIt;
+ }
+
+ // Next Line for output
+ OutIt.NextLine();
+ }
+ }
+
+
+ /**
+ * Method PIC
+ */
+ template<class TImageIn, class TImageOut>
+ void
+ SARTemporalCorrelationGridImageFilter< TImageIn, TImageOut >
+ ::PIC(float TAB_PIC[3][3], float *DX, float *DY, float *VAL_MAX, int *CR) const
+ {
+ float TAB_COS[3]= {-0.5,1.,-0.5};
+ float TAB_SIN[3]= {-0.8660254,0.,+0.8660254};
+ float SXF,SYF,SF,N,A1,B1,SDXF,SDYF,A2,B2;
+ float AX, AY;
+ int i, j;
+
+ // The method of least squares in cosinus
+ *CR = 1;
+ SF = 0.;
+ SXF = 0.;
+ SYF = 0.;
+
+ for (i=-1;i<2;i++)
+ {
+ for (j=-1;j<2;j++)
+ {
+ SF = SF + TAB_PIC[i+1][j+1];
+ SXF = SXF + TAB_PIC[i+1][j+1]*TAB_COS[j+1];
+ SYF = SYF + TAB_PIC[i+1][j+1]*TAB_COS[i+1];
+ }
+ }
+
+ N = SF / 9.;
+ A1 = SXF/ 4.5;
+ B1 = SYF/ 4.5;
+
+ // The method of least squares in sinus
+ SDXF = 0.;
+ SDYF = 0.;
+
+ for (i=-1;i<2;i++)
+ {
+ for (j=-1;j<2;j++)
+ {
+ SDXF = SDXF + TAB_PIC[i+1][j+1]*TAB_SIN[j+1];
+ SDYF = SDYF + TAB_PIC[i+1][j+1]*TAB_SIN[i+1];
+ }
+ }
+
+ A2 = SDXF/3.; // 4*.75
+ B2 = SDYF/3.; // 4*.75
+
+ A2 = A2/1.5;
+ B2 = B2/1.5;
+
+ // Max
+ AX = std::sqrt (A1*A1+A2*A2);
+ AY = std::sqrt (B1*B1+B2*B2);
+
+ if (AX*AY == 0.)
+ {
+ *VAL_MAX = -1.;
+ *DX = 0.;
+ *DY = 0.;
+ return;
+ }
+
+ *DX = std::atan2 (A2,A1) /(2*M_PI/3.);
+ *DY = std::atan2 (B2,B1) /(2*M_PI/3.);
+
+ *VAL_MAX = N + (AX + AY)* 0.5;
+ *VAL_MAX = std::max(*VAL_MAX,TAB_PIC[1][1]);
+ *VAL_MAX = std::min(static_cast<float>(1.),*VAL_MAX);
+
+ // Check values
+ if ((fabs((double)*DX)>.75) || (fabs((double)*DY)>.75))
+ {
+ *VAL_MAX = -3.;
+ *DX = 0.;
+ *DY = 0.;
+ }
+
+ return;
+
+ }
+
+ template<class TImageIn, class TImageOut>
+ double
+ SARTemporalCorrelationGridImageFilter< TImageIn, TImageOut >
+ ::Correlation(const std::vector<double>& master, const std::vector<double>& slave) const
+ {
+ double meanSlave = 0;
+ double meanMaster = 0;
+ double correlationValue = 0;
+
+ // Compute means
+ for(unsigned int i = 0; i < master.size(); ++i)
+ {
+ meanSlave += slave[i];
+ meanMaster += master[i];
+ }
+ meanSlave /= slave.size();
+ meanMaster /= master.size();
+
+ double crossProd = 0.;
+ double squareSumMaster = 0.;
+ double squareSumSlave = 0.;
+
+
+ // Compute correlation
+ for(unsigned int i = 0; i < master.size(); ++i)
+ {
+ crossProd += (slave[i]-meanSlave) * (master[i]-meanMaster);
+ squareSumSlave += (slave[i]-meanSlave) * (slave[i]-meanSlave);
+ squareSumMaster += (master[i]-meanMaster) * (master[i]-meanMaster);
+ }
+
+ correlationValue = std::abs(crossProd/(std::sqrt(squareSumSlave)*std::sqrt(squareSumMaster)));
+
+ return correlationValue;
+ }
+
+ template<class TImageIn, class TImageOut>
+ double
+ SARTemporalCorrelationGridImageFilter< TImageIn, TImageOut >
+ ::CorrelationDiapason(const std::vector<double>& master, const std::vector<double>& slave) const
+ {
+ double sumSlave = 0;
+ double sumMaster = 0;
+ double sumMasterSlave = 0;
+ double sumSlaveSquare = 0;
+ double sumMasterSquare = 0;
+
+ double correlationValue = 0;
+
+ // Compute sum, sum square and cross prod
+ for(unsigned int i = 0; i < master.size(); ++i)
+ {
+ sumSlave += slave[i];
+ sumMaster += master[i];
+ sumMasterSlave += slave[i]*master[i];
+ sumSlaveSquare += slave[i]*slave[i];
+ sumMasterSquare += master[i]*master[i];
+ }
+
+ // Compute means
+ double meanMaster = sumMaster/master.size();
+ double meanSlave = sumSlave/slave.size();
+
+ // Compute sigma and covariance
+ double sigmaMaster = sumMasterSquare/master.size() - (meanMaster*meanMaster);
+ double sigmaSlave = sumSlaveSquare/slave.size() - (meanSlave*meanSlave);
+ double covariance = sumMasterSlave/master.size();
+
+
+
+ if (sigmaMaster > 0 && sigmaSlave > 0)
+ {
+ correlationValue = (covariance - (meanMaster*meanSlave))/(std::sqrt(sigmaMaster*sigmaSlave));
+ }
+
+ return correlationValue;
+ }
+
+
+
+} /*namespace otb*/
+
+#endif