ENH: New application for precise correlation (part 1)

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}
+)

app/otbSARCorrelationPrecise.cxx

+/*
+ * 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)

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;

include/otbSARTemporalCorrelationGridImageFilter.h

+/*
+ * 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