diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 790164a1cf69a3f50aadafc798247c358a1cbb8d..dd339f5ab7bfa4d76c5e1ea0d35c64bb2a34d03f 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -98,3 +98,9 @@ OTB_CREATE_APPLICATION(NAME SARESD
SOURCES otbSARESD.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES}
)
+
+OTB_CREATE_APPLICATION(NAME SARCorrelationRough
+ SOURCES otbSARCorrelationRough.cxx
+ LINK_LIBRARIES ${${otb-module}_LIBRARIES}
+)
+
diff --git a/app/otbSARCorrectionGrid.cxx b/app/otbSARCorrectionGrid.cxx
index 5547de7a3f54a7ddba01d05308bcf2eefd40db6c..75af1927780aa3844cde33abc052e32262b4aab6 100644
--- a/app/otbSARCorrectionGrid.cxx
+++ b/app/otbSARCorrectionGrid.cxx
@@ -24,6 +24,11 @@
#include "otbSARStreamingGridInformationImageFilter.h"
#include "otbSARCorrectionGridImageFilter.h"
+#include "otbMultiToMonoChannelExtractROI.h"
+#include "otbImageList.h"
+#include "otbImageListToVectorImageFilter.h"
+#include "itkAddImageFilter.h"
+
#include <iostream>
#include <string>
#include <fstream>
@@ -44,7 +49,16 @@ public:
// Filters
typedef otb::SARCorrectionGridImageFilter<FloatVectorImageType, FloatVectorImageType> CorrectionGridFilterType;
- typedef otb::SARStreamingGridInformationImageFilter<FloatVectorImageType> MeanFilterType;
+ typedef otb::SARStreamingGridInformationImageFilter<FloatVectorImageType> MeanFilterType;
+ typedef otb::ImageList<FloatImageType> ImageListType;
+ typedef otb::ImageListToVectorImageFilter<ImageListType,
+ FloatVectorImageType > ListConcatenerFilterType;
+ typedef otb::MultiToMonoChannelExtractROI<FloatVectorImageType::InternalPixelType,
+ FloatImageType::PixelType> ExtractROIFilterType;
+ typedef itk::AddImageFilter <FloatImageType, FloatImageType, FloatImageType> AddFilterType;
+
+ typedef ObjectList<ExtractROIFilterType> ExtractROIFilterListType;
+ typedef ObjectList<AddFilterType> AddFilterListType;
private:
@@ -87,6 +101,10 @@ private:
SetMinimumParameterFloatValue("gap", 0.01);
MandatoryOff("gap");
+ AddParameter(ParameterType_String, "advantage", "Give an advantage to DEM or Correlation Grid");
+ SetParameterDescription("advantage","Give an advantage to DEM or Correlation Grid.");
+ MandatoryOff("advantage");
+
AddRAMParameter();
SetDocExampleParameterValue("indemgrid","./demGrid.tiff");
@@ -109,39 +127,156 @@ void DoExecute() override
otbAppLogINFO(<<"Threshold on correlation rate :"<<threshold);
otbAppLogINFO(<<"Gap on the difference between grid values and mean values: "<<gap);
+ std::string advantage = GetParameterString("advantage");
+ bool demGridToFavor = true;
+
+ if (advantage == "correlation" || advantage == "cor" || advantage == "corgrid")
+ {
+ demGridToFavor = false;
+ }
+
+ if (demGridToFavor)
+ {
+ otbAppLogINFO(<<"Opt in favor of projection (DEMGrid)");
+ }
+ else
+ {
+ otbAppLogINFO(<<"Opt in favor of correlation (CorGrid)");
+ }
+
+
// Get input grids
FloatVectorImageType::Pointer DEMGridPtr = GetParameterFloatVectorImage("indemgrid");
FloatVectorImageType::Pointer CorGridPtr = GetParameterFloatVectorImage("incorgrid");
- // Mean Persistent Filter
- MeanFilterType::Pointer filterMean = MeanFilterType::New();
- m_Ref.push_back(filterMean.GetPointer());
-
- // Configure Mean Filter
- filterMean->SetThreshold(threshold);
- filterMean->SetInput(DEMGridPtr);
- filterMean->SetEstimateMean(true);
-
- // Get mean values
- filterMean->Update();
- double meanRange, meanAzimut;
- filterMean->GetMean(meanRange, meanAzimut);
-
- // Correction Filter
- CorrectionGridFilterType::Pointer filterCorrectionGrid = CorrectionGridFilterType::New();
- m_Ref.push_back(filterCorrectionGrid.GetPointer());
+ if (demGridToFavor)
+ {
+ // Mean Persistent Filter
+ MeanFilterType::Pointer filterMean = MeanFilterType::New();
+ m_Ref.push_back(filterMean.GetPointer());
+
+ // Configure Mean Filter
+ filterMean->SetThreshold(threshold);
+ filterMean->SetInput(DEMGridPtr);
+ filterMean->SetEstimateMean(true);
+
+ // Get mean values
+ filterMean->Update();
+ double meanRange, meanAzimut;
+ filterMean->GetMean(meanRange, meanAzimut);
+
+ // Correction Filter
+ CorrectionGridFilterType::Pointer filterCorrectionGrid = CorrectionGridFilterType::New();
+ m_Ref.push_back(filterCorrectionGrid.GetPointer());
+
+ // Configure CorrectionGrid Filter
+ filterCorrectionGrid->SetThreshold(threshold);
+ filterCorrectionGrid->SetGap(gap);
+ filterCorrectionGrid->SetMeanRange(meanRange);
+ filterCorrectionGrid->SetMeanAzimut(meanAzimut);
- // Configure CorrectionGrid Filter
- filterCorrectionGrid->SetThreshold(threshold);
- filterCorrectionGrid->SetGap(gap);
- filterCorrectionGrid->SetMeanRange(meanRange);
- filterCorrectionGrid->SetMeanAzimut(meanAzimut);
+ // Define the main pipeline
+ filterCorrectionGrid->SetDEMGridInput(DEMGridPtr);
+ filterCorrectionGrid->SetCorGridInput(CorGridPtr);
- // Define the main pipeline
- filterCorrectionGrid->SetDEMGridInput(DEMGridPtr);
- filterCorrectionGrid->SetCorGridInput(CorGridPtr);
+ SetParameterOutputImage("out", filterCorrectionGrid->GetOutput());
+ }
+ else
+ {
+ // Mean Persistent Filters (for the two grids)
+ MeanFilterType::Pointer filterMeanDEM = MeanFilterType::New();
+ m_Ref.push_back(filterMeanDEM.GetPointer());
+
+ MeanFilterType::Pointer filterMeanCor = MeanFilterType::New();
+ m_Ref.push_back(filterMeanCor.GetPointer());
+
+ // Configure Mean Filters
+ filterMeanDEM->SetThreshold(threshold);
+ filterMeanDEM->SetInput(DEMGridPtr);
+ filterMeanDEM->SetEstimateMean(true);
+
+ filterMeanCor->SetThreshold(threshold);
+ filterMeanCor->SetInput(CorGridPtr);
+ filterMeanCor->SetEstimateMean(true);
+
+ // Get mean values
+ filterMeanDEM->Update();
+ double meanRangeDEM, meanAzimutDEM;
+ filterMeanDEM->GetMean(meanRangeDEM, meanAzimutDEM);
+
+ filterMeanCor->Update();
+ double meanRangeCor, meanAzimutCor;
+ filterMeanCor->GetMean(meanRangeCor, meanAzimutCor);
+
+ // Difference between means
+ double diffMeanRange = meanRangeCor - meanRangeDEM;
+ double diffMeanAzimut = meanAzimutCor - meanAzimutDEM;
+
+ std::cout << "diffMeanRange = " << diffMeanRange << std::endl;
+ std::cout << "diffMeanAzimut = " << diffMeanAzimut << std::endl;
+
+ // Offset for DEMGrid (Correction with the difference)
+ // Instanciate all filters
+ ListConcatenerFilterType::Pointer m_Concatener =
+ ListConcatenerFilterType::New();
+ ExtractROIFilterListType::Pointer m_ExtractorList =
+ ExtractROIFilterListType::New();
+ AddFilterListType::Pointer m_AddOffsetList =
+ AddFilterListType::New();
+ ImageListType::Pointer m_ImageList =
+ ImageListType::New();
+
+ // Check number of components (at least two : range and azimut)
+ unsigned int nbComponents = DEMGridPtr->GetNumberOfComponentsPerPixel();
+ if (nbComponents < 2)
+ {
+ otbAppLogFATAL(<< "Input grid has to contain at least two components");
+ }
+
+ // Extract eah band and apply the corresponding offset
+ for (unsigned int i = 0; i < nbComponents; ++i)
+ {
+ // Instanciate filters
+ ExtractROIFilterType::Pointer extractor =
+ ExtractROIFilterType::New();
+ AddFilterType::Pointer addOffset =
+ AddFilterType::New();
+
+ // Set the channel to extract
+ extractor->SetInput(DEMGridPtr);
+ extractor->SetChannel(i+1);
+ extractor->UpdateOutputInformation();
+ m_ExtractorList->PushBack(extractor);
+
+ addOffset->SetInput(extractor->GetOutput());
+
+ // Adapt offset according to the numnber of band
+ if (i == 0)
+ {
+ addOffset->SetConstant2(diffMeanRange);
+ }
+ else if (i == 1)
+ {
+ addOffset->SetConstant2(diffMeanAzimut);
+ }
+ else
+ {
+ addOffset->SetConstant2(0);
+ }
+
+
+ m_AddOffsetList->PushBack(addOffset);
+ m_ImageList->PushBack(addOffset->GetOutput() );
+
+ }
- SetParameterOutputImage("out", filterCorrectionGrid->GetOutput());
+ m_Concatener->SetInput( m_ImageList );
+
+ SetParameterOutputImage("out", m_Concatener->GetOutput());
+ RegisterPipeline();
+ //SetParameterOutputImage("out", DEMGridPtr);
+
+ }
}
// Vector for filters
diff --git a/app/otbSARCorrelationGrid.cxx b/app/otbSARCorrelationGrid.cxx
index b908e2946afdefd6354871efff10aff6f913ba41..7976294dc787d51ade71678a69debb13785e7598 100644
--- a/app/otbSARCorrelationGrid.cxx
+++ b/app/otbSARCorrelationGrid.cxx
@@ -21,7 +21,10 @@
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
-#include "otbSARCorrelationGridImageFilter.h"
+#include "itkFFTNormalizedCorrelationImageFilter.h"
+#include "itkMinimumMaximumImageCalculator.h"
+
+#include "otbSARTemporalCorrelationGridImageFilter.h"
#include <iostream>
#include <string>
@@ -42,22 +45,23 @@ public:
itkTypeMacro(SARCorrelationGrid, otb::Wrapper::Application);
// Filters
- typedef otb::SARCorrelationGridImageFilter<FloatImageType, FloatVectorImageType> CorGridFilterType;
+ typedef itk::FFTNormalizedCorrelationImageFilter<FloatImageType, FloatImageType> CorFilterType;
+ typedef itk::MinimumMaximumImageCalculator< FloatImageType> MinMaxCalculatorType;
+ typedef otb::SARTemporalCorrelationGridImageFilter<FloatImageType, FloatVectorImageType> CorGridFilterType;
private:
void DoInit() override
{
SetName("SARCorrelationGrid");
- SetDescription("Computes SAR correlation grid.");
+ SetDescription("Computes SAR correlation shift (into temporal domain).");
- SetDocLongDescription("This application creates a deformation grid based "
- "on correlation estimation.\nThe output grid is a VectorImage composed of"
- " three values : shift in range, shift in azimut and a correlation rate. "
- "The inputs of this application are MultiLooked images (real images).");
+ 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 are supported for now.");
+ SetDocLimitations("Only Sentinel 1 products and Cosmo are supported for now.");
SetDocAuthors("OTB-Team");
SetDocSeeAlso(" ");
AddDocTag(Tags::SAR);
@@ -69,13 +73,6 @@ private:
AddParameter(ParameterType_InputImage, "inslave", "Input Slave image (real image)");
SetParameterDescription("inslave", "Slave Image (real image).");
- AddParameter(ParameterType_InputFilename, "indem", "Input DEM directory");
- SetParameterDescription("indem", "Input DEM directory to help projection.");
- MandatoryOff("indem");
-
- AddParameter(ParameterType_OutputImage, "out", "Output Correlation grid (Vector Image)");
- SetParameterDescription("out","Output Correlation Grid Vector Image (Shift_ran, Shift_azi, Correlation_rate).");
-
AddParameter(ParameterType_Int, "mlran", "MultiLook factor on distance");
SetParameterDescription("mlran","MultiLook factor on distance.");
SetDefaultParameterInt("mlran", 3);
@@ -87,12 +84,6 @@ private:
SetDefaultParameterInt("mlazi", 3);
SetMinimumParameterIntValue("mlazi", 1);
MandatoryOff("mlazi");
-
- AddParameter(ParameterType_Int, "patchsize", "Size of patch into each dimension for correlation estimation");
- SetParameterDescription("patchsize","Size of patch into each dimension for correlation estimation.");
- SetDefaultParameterInt("patchsize", 20);
- SetMinimumParameterIntValue("patchsize", 1);
- MandatoryOff("patchsize");
AddParameter(ParameterType_Int, "gridsteprange", "Grid step for range dimension (into SLC/SAR geometry)");
SetParameterDescription("gridsteprange","Grid step for range dimension (into SLC/SAR geometry).");
@@ -106,18 +97,21 @@ private:
SetMinimumParameterIntValue("gridstepazimut", 1);
MandatoryOff("gridstepazimut");
- AddParameter(ParameterType_Int,"subpixel","Subpixel measure");
- SetParameterDescription("subpixel", "Subpixel measure.");
- SetDefaultParameterInt("subpixel", 10);
- SetMinimumParameterIntValue("subpixel", 1);
- MandatoryOff("subpixel");
+ AddParameter(ParameterType_Int, "nooffset", "Set 0 to offset of first Line and Colunm of output grid");
+ SetParameterDescription("nooffset", "If 1, then no offset for the first L/C if output grid.");
+ SetDefaultParameterInt("nooffset", 0);
+ SetMinimumParameterIntValue("nooffset", 0);
+ MandatoryOff("nooffset");
+ // 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("indem","./DEMDir");
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","corGrid.tif");
+ SetDocExampleParameterValue("out", "out_CorrelationGrid.tiff");
}
void DoUpdateParameters() override
@@ -132,38 +126,57 @@ void DoExecute() override
int factorML_ran = GetParameterInt("mlran");
int grid_step_azi = GetParameterInt("gridstepazimut");
int grid_step_ran = GetParameterInt("gridsteprange");
- int patch_size = GetParameterInt("patchsize");
- int subpixel = GetParameterInt("subpixel");
+ int nooffset = GetParameterInt("nooffset");
// 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);
- otbAppLogINFO(<<"Patch Size : "<<patch_size);
- otbAppLogINFO(<<"SubPixel Measure : "<<subpixel);
// Get master and slave image
FloatImageType::Pointer MasterPtr = GetParameterFloatImage("inmaster");
FloatImageType::Pointer SlavePtr = GetParameterFloatImage("inslave");
+
- // CorrelationGrid Filter
+ // 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->SetPatchSizePerDim(patch_size);
filterCorrelationGrid->SetGridStep(grid_step_ran, grid_step_azi);
- filterCorrelationGrid->SetSubPixelAccuracy(subpixel);
- // Set DEM Directory if required
- if (GetParameterByKey("indem")->HasValue())
+ filterCorrelationGrid->SetRoughShift_ran(shiftML_range);
+ filterCorrelationGrid->SetRoughShift_azi(shiftML_azimut);
+ if (nooffset)
{
- std::string demDir = GetParameterString("indem");
- filterCorrelationGrid->activateDEM(demDir);
+ filterCorrelationGrid->SetFirstLCOffset(false);
}
+
// Define the main pipeline
filterCorrelationGrid->SetMasterInput(MasterPtr);
filterCorrelationGrid->SetSlaveInput(SlavePtr);
diff --git a/app/otbSARCorrelationRough.cxx b/app/otbSARCorrelationRough.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..39f9e00c0a841fa6909933b3e3e51917261dbafb
--- /dev/null
+++ b/app/otbSARCorrelationRough.cxx
@@ -0,0 +1,169 @@
+/*
+ * 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 <iostream>
+#include <string>
+#include <fstream>
+
+namespace otb
+{
+namespace Wrapper
+{
+
+class SARCorrelationRough : public Application
+{
+public:
+ typedef SARCorrelationRough Self;
+ typedef itk::SmartPointer<Self> Pointer;
+
+ itkNewMacro(Self);
+ itkTypeMacro(SARCorrelationRough, otb::Wrapper::Application);
+
+ // Filters
+typedef itk::FFTNormalizedCorrelationImageFilter<FloatImageType, FloatImageType> CorFilterType;
+typedef itk::MinimumMaximumImageCalculator< FloatImageType> MinMaxCalculatorType;
+
+private:
+ void DoInit() override
+ {
+ SetName("SARCorrelationRough");
+ SetDescription("Computes SAR correlation shift (rough shifts).");
+
+ SetDocLongDescription("This application comutes rough 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");
+
+ // Parameter Output
+ 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();
+
+ 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");
+ }
+
+ 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");
+
+ // Log information
+ otbAppLogINFO(<<"ML Factor on azimut :"<<factorML_azi);
+ otbAppLogINFO(<<"ML Factor on range : "<<factorML_ran);
+
+ // 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]);
+
+ // Assigne Outputs
+ 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;
+};
+
+
+}
+
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::SARCorrelationRough)
diff --git a/app/otbSARFineDeformationGrid.cxx b/app/otbSARFineDeformationGrid.cxx
index 3353a57490f688220433d82001c6322ae3c7dcf3..4686a1eeddb8b8fb26214860b51de9624d1b76ae 100644
--- a/app/otbSARFineDeformationGrid.cxx
+++ b/app/otbSARFineDeformationGrid.cxx
@@ -107,12 +107,6 @@ private:
SetMinimumParameterIntValue("mlazi", 1);
MandatoryOff("mlazi");
- AddParameter(ParameterType_Int, "patchsize", "Size of patch into each dimension for correlation estimation");
- SetParameterDescription("patchsize","Size of patch into each dimension for correlation estimation.");
- SetDefaultParameterInt("patchsize", 20);
- SetMinimumParameterIntValue("patchsize", 1);
- MandatoryOff("patchsize");
-
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);
@@ -125,12 +119,6 @@ private:
SetMinimumParameterIntValue("gridstepazimut", 1);
MandatoryOff("gridstepazimut");
- AddParameter(ParameterType_Int,"subpixel","Subpixel measure");
- SetParameterDescription("subpixel", "Subpixel measure.");
- SetDefaultParameterInt("subpixel", 10);
- SetMinimumParameterIntValue("subpixel", 1);
- MandatoryOff("subpixel");
-
AddParameter(ParameterType_Float, "threshold", "Threshold for correlation rate");
SetParameterDescription("threshold","Threshold for correlation rate.");
SetDefaultParameterFloat("threshold", 0.3);
@@ -142,6 +130,10 @@ private:
SetDefaultParameterFloat("gap", 0.7);
SetMinimumParameterFloatValue("gap", 0.01);
MandatoryOff("gap");
+
+ AddParameter(ParameterType_String, "advantage", "Give an advantage to DEM or Correlation Grid");
+ SetParameterDescription("advantage","Give an advantage to DEM or Correlation Grid.");
+ MandatoryOff("advantage");
SetDocExampleParameterValue("indem","S21E055.hgt");
SetDocExampleParameterValue("insarmaster","s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002_SLC.tiff");
@@ -171,21 +163,39 @@ private:
int factor_ran = GetParameterInt("mlran");
int grid_step_azi = GetParameterInt("gridstepazimut");
int grid_step_ran = GetParameterInt("gridsteprange");
- int patch_size = GetParameterInt("patchsize");
- int subpixel = GetParameterInt("subpixel");
float threshold = GetParameterFloat("threshold");
float gap = GetParameterFloat("gap");
+
+ std::string advantage = GetParameterString("advantage");
+ bool demGridToFavor = true;
+
+ if (advantage == "correlation" || advantage == "cor" || advantage == "corgrid")
+ {
+ demGridToFavor = false;
+ }
+ if (demGridToFavor)
+ {
+ otbAppLogINFO(<<"Opt in favor of projection (DEMGrid)");
+ }
+ else
+ {
+ otbAppLogINFO(<<"Opt in favor of correlation (CorGrid)");
+ }
+
+
// Execute the Pipelone with internal applications
// First application : SARCorrelationGrid
GetInternalApplication("CorGridApp")->SetParameterInputImage("inmaster", inputMasterML);
GetInternalApplication("CorGridApp")->SetParameterInputImage("inslave", inputSlaveML);
GetInternalApplication("CorGridApp")->SetParameterInt("mlran", factor_ran);
GetInternalApplication("CorGridApp")->SetParameterInt("mlazi", factor_azi);
- GetInternalApplication("CorGridApp")->SetParameterInt("patchsize", patch_size);
GetInternalApplication("CorGridApp")->SetParameterInt("gridsteprange", grid_step_ran);
GetInternalApplication("CorGridApp")->SetParameterInt("gridstepazimut", grid_step_azi);
- GetInternalApplication("CorGridApp")->SetParameterInt("subpixel", subpixel);
+ if (demGridToFavor)
+ {
+ GetInternalApplication("CorGridApp")->SetParameterInt("nooffset", 1);
+ }
ExecuteInternal("CorGridApp");
// Second application : SARDEMGrid
@@ -206,13 +216,14 @@ private:
GetInternalApplication("DEMGridApp")->SetParameterInt("gridstepazimut", grid_step_azi);
ExecuteInternal("DEMGridApp");
- // Third (and last) application : SARMultiLook
+ // Third (and last) application : SARMultiLook
GetInternalApplication("CorrectionGridApp")->SetParameterInputImage("indemgrid",
GetInternalApplication("DEMGridApp")->GetParameterOutputImage("out"));
GetInternalApplication("CorrectionGridApp")->SetParameterInputImage("incorgrid",
GetInternalApplication("CorGridApp")->GetParameterOutputImage("out"));
GetInternalApplication("CorrectionGridApp")->SetParameterFloat("threshold", threshold);
GetInternalApplication("CorrectionGridApp")->SetParameterFloat("gap", gap);
+ GetInternalApplication("CorrectionGridApp")->SetParameterString("advantage", advantage);
ExecuteInternal("CorrectionGridApp");
// Retrive the output of CorrectionGridApp
diff --git a/app/otbSARFineMetadata.cxx b/app/otbSARFineMetadata.cxx
index fb6069706f7c9727ad4b206a907ad6576b20097a..308af0d93a8b3e4399eaec32bed543f2ee724f74 100644
--- a/app/otbSARFineMetadata.cxx
+++ b/app/otbSARFineMetadata.cxx
@@ -22,6 +22,7 @@
#include "otbWrapperApplicationFactory.h"
#include "otbSARStreamingGridHistogramFilter.h"
+#include "otbSARUpdateMetadataImageFilter.h"
#include <iostream>
#include <string>
@@ -30,11 +31,20 @@
#include <cmath>
// include ossim
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
#include "ossim/ossimTimeUtilities.h"
#include "ossim/base/ossimKeywordlist.h"
#include "ossim/base/ossimString.h"
#include "ossim/base/ossimDate.h"
-
+#pragma GCC diagnostic pop
+#else
+#include "ossim/ossimTimeUtilities.h"
+#include "ossim/base/ossimKeywordlist.h"
+#include "ossim/base/ossimString.h"
+#include "ossim/base/ossimDate.h"
+#endif
namespace otb
{
@@ -52,6 +62,7 @@ public:
// Filter
typedef otb::SARStreamingGridHistogramFilter<FloatVectorImageType> GridHistogramType;
+ typedef otb::SARUpdateMetadataImageFilter<ComplexFloatImageType> UpdateMetadataeFilter;
// Time/Date ossim
typedef ossimplugins::time::ModifiedJulianDate TimeType;
@@ -101,6 +112,9 @@ private:
AddParameter(ParameterType_Float,"slantrange","Slant near range (output of this application)");
SetParameterDescription("slantrange", "Slant near range (output of this application).");
SetParameterRole("slantrange", Role_Output);
+
+ AddParameter(ParameterType_OutputImage, "out", "Output image with precise metadata");
+ SetParameterDescription("out","Output image with precise metadata.");
AddRAMParameter();
@@ -214,7 +228,15 @@ private:
SetParameterString("timefirstline", fineTimeFirstLine);
SetParameterFloat("slantrange", fineSlantRange);
-
+
+ // Copie input image and update the metadata into a naw geom file
+ UpdateMetadataeFilter::Pointer updateMetadataFilter = UpdateMetadataeFilter::New();
+ m_Ref.push_back(updateMetadataFilter.GetPointer());
+ updateMetadataFilter->SetTimeFirstLine(fineTimeFirstLine);
+ updateMetadataFilter->SetSlantRange(fineSlantRange);
+ updateMetadataFilter->SetInput(SARPtr);
+
+ SetParameterOutputImage("out", updateMetadataFilter->GetOutput());
}
// Vector for filters
diff --git a/include/otbSARDEMPolygonsAnalysisImageFilter.txx b/include/otbSARDEMPolygonsAnalysisImageFilter.txx
index 388997225c0fde787827516834c11662431cda93..34229bed64a14ead442d4401d8db7476cecbcf27 100644
--- a/include/otbSARDEMPolygonsAnalysisImageFilter.txx
+++ b/include/otbSARDEMPolygonsAnalysisImageFilter.txx
@@ -179,9 +179,13 @@ namespace otb
// Set the margin
int nbLinesSAR = m_SarImagePtr->GetLargestPossibleRegion().GetSize()[1];
int nbLinesDEM = m_DemImagePtr->GetLargestPossibleRegion().GetSize()[1];
- if (nbLinesSAR > nbLinesDEM)
+ if (nbLinesSAR/2 > nbLinesDEM)
{
- m_Margin = (nbLinesSAR/nbLinesDEM)*5;
+ m_Margin = (nbLinesSAR/nbLinesDEM)*30;
+ }
+ else if (nbLinesSAR > nbLinesDEM)
+ {
+ m_Margin = (nbLinesSAR/nbLinesDEM)*60;
}
else
{
@@ -1102,6 +1106,7 @@ namespace otb
while (!InIt_UR.IsAtEndOfLine() && !InIt_LR.IsAtEndOfLine())
{
UL_line = InIt_UL.Get()[1];
+
// Protect Memory access and input access
if (UL_line > ind_Line_Into_SARGeo - m_Margin && UL_line < ind_Line_Into_SARGeo + m_Margin)
{
@@ -1109,7 +1114,7 @@ namespace otb
UR_col = InIt_UR.Get()[0];
LR_col = InIt_LR.Get()[0];
LL_col = InIt_LL.Get()[0];
-
+
// Check if No Data
if (!(UL_col == m_NoData || UR_col == m_NoData || LR_col == m_NoData
|| LL_col == m_NoData))
diff --git a/include/otbSARDEMProjectionImageFilter.txx b/include/otbSARDEMProjectionImageFilter.txx
index da9ba196b899bef948fae568b2b09453f4ead39c..d70299e22e0f2620fcf878fc91f5f83f35d80945 100644
--- a/include/otbSARDEMProjectionImageFilter.txx
+++ b/include/otbSARDEMProjectionImageFilter.txx
@@ -29,8 +29,16 @@
#include "itkProgressReporter.h"
#include "itkNumericTraitsPointPixel.h"
+#if defined(__GNUC__) || defined(__clang__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
#include "ossim/base/ossimFilename.h"
#include "ossim/base/ossimGpt.h"
+#pragma GCC diagnostic pop
+#else
+#include "ossim/base/ossimFilename.h"
+#include "ossim/base/ossimGpt.h"
+#endif
#include <cmath>
diff --git a/include/otbSARDopplerCentroidFreqImageFilter.txx b/include/otbSARDopplerCentroidFreqImageFilter.txx
index 593e81b4ad02821f88bc4daab3ed0ea5d999cc63..a401dfb3f330c67b3e8a4e2e3a9bf77b328c5bb2 100644
--- a/include/otbSARDopplerCentroidFreqImageFilter.txx
+++ b/include/otbSARDopplerCentroidFreqImageFilter.txx
@@ -90,7 +90,7 @@ namespace otb
for (std::size_t listId=0; listId!= nbLists ; ++listId)
{
- const int pos = sprintf(fmRatePrefix_, "%s%d.", FM_PREFIX.c_str(), (listId+1));
+ const int pos = sprintf(fmRatePrefix_, "%s%zu.", FM_PREFIX.c_str(), (listId+1));
assert(pos > 0 && pos < sizeof(fmRatePrefix_));
const std::string FMPrefix(fmRatePrefix_, pos);
@@ -120,7 +120,7 @@ namespace otb
for (std::size_t listId=0; listId!= nbLists ; ++listId)
{
- const int pos = sprintf(dcfPrefix_, "%s%d.", DCF_PREFIX.c_str(), (listId+1));
+ const int pos = sprintf(dcfPrefix_, "%s%zu.", DCF_PREFIX.c_str(), (listId+1));
assert(pos > 0 && pos < sizeof(dcfPrefix_));
const std::string DCFPrefix(dcfPrefix_, pos);
diff --git a/include/otbSARESDOffsetsImageFilter.txx b/include/otbSARESDOffsetsImageFilter.txx
index 90dbcab741dc7ee29c8b17116d93a80fd91b9756..fc47727cb249bfb0b4b64ec99884a04e85c611d3 100644
--- a/include/otbSARESDOffsetsImageFilter.txx
+++ b/include/otbSARESDOffsetsImageFilter.txx
@@ -184,9 +184,6 @@ namespace otb
// For each colunm
while (!OutIt.IsAtEndOfLine() && !In1It.IsAtEndOfLine() && !In2It.IsAtEndOfLine())
{
- // Output index
- ImageIndexType index = OutIt.GetIndex();
-
double diff_DopCentroid = In1It.Get() - In2It.Get();
// Shift for azimut dimension
diff --git a/include/otbSARCorrelationGridImageFilter.h b/include/otbSARFFTCorrelationGridImageFilter.h
similarity index 91%
rename from include/otbSARCorrelationGridImageFilter.h
rename to include/otbSARFFTCorrelationGridImageFilter.h
index 2beff8226b3fbd291f17405846f70500cc20f5ed..35533c5c206522e0811fb578faabd5d09e6cc54d 100644
--- a/include/otbSARCorrelationGridImageFilter.h
+++ b/include/otbSARFFTCorrelationGridImageFilter.h
@@ -18,8 +18,8 @@
* limitations under the License.
*/
-#ifndef otbSARCorrelationGridImageFilter_h
-#define otbSARCorrelationGridImageFilter_h
+#ifndef otbSARFFTCorrelationGridImageFilter_h
+#define otbSARFFTCorrelationGridImageFilter_h
#include "itkImageToImageFilter.h"
#include "itkSmartPointer.h"
@@ -35,22 +35,23 @@
namespace otb
{
- /** \class SARCorrelationGridImageFilter
+ /** \class SARFFTCorrelationGridImageFilter
* \brief Estimates a deformation grid between two images (master and slave) for the azimut and range
* dimension. A correlation rate is also calculated. The DEM can be used.
*
- * This filter performs the estimation of deformations between two images (the images are reals).
+ * This filter performs the estimation of deformations between two images (the images are reals) in spectral
+ * domain.
*
* \ingroup DiapOTBModule
*/
- template <typename TImageIn, typename TImageOut> class ITK_EXPORT SARCorrelationGridImageFilter :
+ template <typename TImageIn, typename TImageOut> class ITK_EXPORT SARFFTCorrelationGridImageFilter :
public itk::ImageToImageFilter<TImageIn,TImageOut>
{
public:
// Standard class typedefs
- typedef SARCorrelationGridImageFilter Self;
+ typedef SARFFTCorrelationGridImageFilter Self;
typedef itk::ImageToImageFilter<TImageIn,TImageOut> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
@@ -169,23 +170,23 @@ namespace otb
protected:
// Constructor
- SARCorrelationGridImageFilter();
+ SARFFTCorrelationGridImageFilter();
// Destructor
- virtual ~SARCorrelationGridImageFilter() ITK_OVERRIDE {};
+ virtual ~SARFFTCorrelationGridImageFilter() ITK_OVERRIDE {};
// Print
void PrintSelf(std::ostream & os, itk::Indent indent) const ITK_OVERRIDE;
- /** SARCorrelationGridImageFilter produces an vector image to indicate the deformation for the two dimensions
+ /** SARFFTCorrelationGridImageFilter 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, SARCorrelationGridImageFilter needs to provide an implementation for
+ * As such, SARFFTCorrelationGridImageFilter needs to provide an implementation for
* GenerateOutputInformation() in order to inform the pipeline execution model.
*/
virtual void GenerateOutputInformation() ITK_OVERRIDE;
- /** SARCorrelationGridImageFilter needs input requested regions (for master and slave images) that
+ /** SARFFTCorrelationGridImageFilter 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
@@ -198,7 +199,7 @@ namespace otb
*/
void OutputRegionToInputRegion(const ImageOutRegionType& outputRegion) const;
- /** SARCorrelationGridImageFilter can be implemented as a multithreaded filter.
+ /** SARFFTCorrelationGridImageFilter 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
@@ -218,7 +219,7 @@ namespace otb
private:
- SARCorrelationGridImageFilter(const Self&); // purposely not implemented
+ SARFFTCorrelationGridImageFilter(const Self&); // purposely not implemented
void operator=(const Self &); // purposely not
/** Coarse correlation threshold */
@@ -268,7 +269,7 @@ namespace otb
} // End namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
-#include "otbSARCorrelationGridImageFilter.txx"
+#include "otbSARFFTCorrelationGridImageFilter.txx"
#endif
diff --git a/include/otbSARCorrelationGridImageFilter.txx b/include/otbSARFFTCorrelationGridImageFilter.txx
similarity index 93%
rename from include/otbSARCorrelationGridImageFilter.txx
rename to include/otbSARFFTCorrelationGridImageFilter.txx
index 6300fe1c5597ecaea224f490b59c925a2bf62bac..9b064c17a118a03ffaf1e157b2f6500c3887c5a7 100644
--- a/include/otbSARCorrelationGridImageFilter.txx
+++ b/include/otbSARFFTCorrelationGridImageFilter.txx
@@ -18,10 +18,10 @@
* limitations under the License.
*/
-#ifndef otbSARCorrelationGridImageFilter_txx
-#define otbSARCorrelationGridImageFilter_txx
+#ifndef otbSARFFTCorrelationGridImageFilter_txx
+#define otbSARFFTCorrelationGridImageFilter_txx
-#include "otbSARCorrelationGridImageFilter.h"
+#include "otbSARFFTCorrelationGridImageFilter.h"
#include "itkImageScanlineConstIterator.h"
#include "itkImageRegionIterator.h"
@@ -45,7 +45,7 @@ namespace otb
* Constructor with default initialization
*/
template <class TImageIn, class TImageOut>
- SARCorrelationGridImageFilter< TImageIn ,TImageOut >::SARCorrelationGridImageFilter()
+ SARFFTCorrelationGridImageFilter< TImageIn ,TImageOut >::SARFFTCorrelationGridImageFilter()
{
this->SetNumberOfRequiredInputs(2);
@@ -70,7 +70,7 @@ namespace otb
*/
template <class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn ,TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn ,TImageOut >
::SetMasterInput(const ImageInType* image )
{
// Process object is not const-correct so the const casting is required.
@@ -82,7 +82,7 @@ namespace otb
*/
template <class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn ,TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn ,TImageOut >
::SetSlaveInput(const ImageInType* image)
{
// Process object is not const-correct so the const casting is required.
@@ -93,8 +93,8 @@ namespace otb
* Get Master Image
*/
template <class TImageIn, class TImageOut>
- const typename SARCorrelationGridImageFilter< TImageIn ,TImageOut>::ImageInType *
- SARCorrelationGridImageFilter< TImageIn ,TImageOut >
+ const typename SARFFTCorrelationGridImageFilter< TImageIn ,TImageOut>::ImageInType *
+ SARFFTCorrelationGridImageFilter< TImageIn ,TImageOut >
::GetMasterInput() const
{
if (this->GetNumberOfInputs()<1)
@@ -108,8 +108,8 @@ namespace otb
* Get Slave Image
*/
template <class TImageIn, class TImageOut>
- const typename SARCorrelationGridImageFilter< TImageIn ,TImageOut>::ImageInType *
- SARCorrelationGridImageFilter< TImageIn ,TImageOut >
+ const typename SARFFTCorrelationGridImageFilter< TImageIn ,TImageOut>::ImageInType *
+ SARFFTCorrelationGridImageFilter< TImageIn ,TImageOut >
::GetSlaveInput() const
{
if (this->GetNumberOfInputs()<2)
@@ -124,7 +124,7 @@ namespace otb
*/
template<class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn, TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn, TImageOut >
::PrintSelf(std::ostream & os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
@@ -140,7 +140,7 @@ namespace otb
*/
template<class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn, TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn, TImageOut >
::activateDEM(const std::string DEMDirectory)
{
DEMHandlerPointerType DEMHandler = DEMHandler::Instance();
@@ -153,7 +153,7 @@ namespace otb
**/
template<class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn, TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn, TImageOut >
::GenerateOutputInformation()
{
// Call superclass implementation
@@ -269,7 +269,7 @@ namespace otb
*/
template<class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn, TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn, TImageOut >
::OutputRegionToInputRegion(const ImageOutRegionType& outputRegion) const
{
// Get pointers to inputs
@@ -302,7 +302,7 @@ namespace otb
//////////////////// Requested region for slave ////////////////////
//int max_shift = 10;
- int max_shift = 50;
+ int max_shift = 200;
int firstL, firstC, lastL, lastC;
// Transform back into slave region index space with an maximum shift between master and slave
ImageInIndexType slaveIndex;
@@ -350,7 +350,7 @@ namespace otb
masterRequestedRegion.PadByRadius( m_PatchSizePerDim / 2 );
slaveRequestedRegion.PadByRadius( m_PatchSizePerDim / 2 );
- // Crop the fixed region at the fixed's largest possible region
+ // Crop the fixed region at the fixed's largest possible region (or buffered)
if ( masterRequestedRegion.Crop(masterPtr->GetLargestPossibleRegion()))
{
masterPtr->SetRequestedRegion( masterRequestedRegion );
@@ -400,7 +400,7 @@ namespace otb
*/
template<class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn, TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn, TImageOut >
::GenerateInputRequestedRegion()
{
// Call the superclass' implementation of this method
@@ -416,7 +416,7 @@ namespace otb
*/
template<class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn, TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn, TImageOut >
::BeforeThreadedGenerateData()
{
// Get the number of threads
@@ -511,6 +511,7 @@ namespace otb
m_slaveExtractPieceTab[i]->SetRegions(slaveRegion);
m_slaveExtractPieceTab[i]->Allocate();
}
+
}
/**
@@ -518,10 +519,10 @@ namespace otb
*/
template<class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn, TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn, TImageOut >
::AfterThreadedGenerateData()
{
- // Get the number of threads
+ // Get the number of threads
unsigned int nbThreads = this->GetNumberOfThreads();
for (unsigned int i = 0; i < nbThreads; i++)
@@ -558,7 +559,7 @@ namespace otb
*/
template<class TImageIn, class TImageOut>
void
- SARCorrelationGridImageFilter< TImageIn, TImageOut >
+ SARFFTCorrelationGridImageFilter< TImageIn, TImageOut >
::ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread,
itk::ThreadIdType threadId)
{
@@ -670,6 +671,25 @@ namespace otb
slaveCurrentRegion.SetIndex(currentSlaveIndex);
}
+ // Check Size
+ if (slaveCurrentRegion.GetIndex()[0] + slaveCurrentRegion.GetSize()[0] > slavePtr->GetLargestPossibleRegion().GetSize()[0]
+ || slaveCurrentRegion.GetIndex()[1] + slaveCurrentRegion.GetSize()[1] > slavePtr->GetLargestPossibleRegion().GetSize()[1])
+ {
+ for(unsigned int i = 0; i < ImageInType::ImageDimension; i++)
+ {
+ if (slaveCurrentRegion.GetIndex()[i] + slaveCurrentRegion.GetSize()[i] > slavePtr->GetLargestPossibleRegion().GetSize()[i])
+ {
+ currentSlaveSize[i] = slavePtr->GetLargestPossibleRegion().GetSize()[i] - slaveCurrentRegion.GetIndex()[i];
+
+ if (slaveCurrentRegion.GetIndex()[i] >= slavePtr->GetLargestPossibleRegion().GetSize()[i])
+ {
+ currentSlaveSize[i] = 0;
+ }
+ }
+ }
+ slaveCurrentRegion.SetSize(currentSlaveSize);
+ }
+
// Rescale values to adapt the deformation
double rescale_range = masterCurrentRegion.GetSize()[0];
double rescale_azimut = masterCurrentRegion.GetSize()[1];
diff --git a/include/otbSARTemporalCorrelationGridImageFilter.h b/include/otbSARTemporalCorrelationGridImageFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..76002481bdb0f44c8c084a4a0ebce6499b362158
--- /dev/null
+++ b/include/otbSARTemporalCorrelationGridImageFilter.h
@@ -0,0 +1,260 @@
+/*
+ * 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);
+
+ itkSetMacro(FirstLCOffset, bool);
+ itkGetMacro(FirstLCOffset, bool);
+
+ // 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);
+
+ /** 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;
+
+ // Boolean to add an offset on the first L/C on the grid
+ bool m_FirstLCOffset;
+
+ };
+
+} // 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..38191291b28e81fe7b017e0fef437c6d08ba991d
--- /dev/null
+++ b/include/otbSARTemporalCorrelationGridImageFilter.txx
@@ -0,0 +1,864 @@
+/*
+ * 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_FirstLCOffset = true;
+
+ 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 with an offet
+ int firstC_master = 0;
+ int firstL_master = 0;
+ if (m_FirstLCOffset)
+ {
+ 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];
+ }
+
+ firstL_master = firstL + (lastL - firstL + 1 - m_GridStep_ML[1] * (outputSize[1])) / 2;
+ firstC_master = firstC + (lastC - firstC + 1 - m_GridStep_ML[0] * (outputSize[0])) / 2;
+ }
+ else
+ {
+ 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]));
+
+ if (m_FirstLCOffset)
+ {
+ // Add first L/C into SLC geometry
+ float firstC_SLC = (firstC_master)*m_MLran;
+ float firstL_SLC = (firstL_master)*m_MLazi;
+
+ outputKWL.AddKey("support_data.grid.first_row", std::to_string(firstL_SLC));
+ outputKWL.AddKey("support_data.grid.first_column", std::to_string(firstC_SLC));
+ }
+
+ outputPtr->SetImageKeywordList(outputKWL);
+ }
+
+ /**
+ * Method OutputRegionToInputRegion for GenerateInputRequestedRegion
+ */
+ template<class TImageIn, class TImageOut>
+ void
+ SARTemporalCorrelationGridImageFilter< TImageIn, TImageOut >
+ ::OutputRegionToInputRegion(const ImageOutRegionType& outputRegion)
+ {
+ // Get pointers to inputs
+ ImageInType * masterPtr = const_cast<ImageInType * >(GetMasterInput());
+ ImageInType * slavePtr = const_cast<ImageInType * >(GetSlaveInput());
+
+ ImageOutPointer outputPtr = this->GetOutput();
+
+ 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 = 0;
+ int firstL_Grid = 0;
+
+ if (m_FirstLCOffset)
+ {
+ 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);
+
+
+ firstC_Grid = firstC + (lastC - firstC + 1 - m_GridStep_ML[0] *
+ (outputPtr->GetLargestPossibleRegion().GetSize()[0])) / 2;
+ firstL_Grid = firstL + (lastL - firstL + 1 - m_GridStep_ML[1] *
+ (outputPtr->GetLargestPossibleRegion().GetSize()[1])) / 2;
+ }
+
+ // 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);
+ masterRequestedSize[1] += (2*m_WinCor_azi+1);
+
+ 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();
+
+ ImageOutPointer outPtr = this->GetOutput();
+
+ // 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_Grid = 0;
+ int firstL_Grid = 0;
+
+ if (m_FirstLCOffset)
+ {
+ 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);
+
+ firstC_Grid = firstC + (lastC - firstC + 1 - m_GridStep_ML[0] *
+ (outPtr->GetLargestPossibleRegion().GetSize()[0])) / 2;
+ firstL_Grid = firstL + (lastL - firstL + 1 - m_GridStep_ML[1] *
+ (outPtr->GetLargestPossibleRegion().GetSize()[1])) / 2;
+ }
+
+ // 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_Grid);
+ masterIndex[1] -= (m_WinCor_azi-firstL_Grid);
+
+ // 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
diff --git a/include/otbSARUpdateMetadataImageFilter.h b/include/otbSARUpdateMetadataImageFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..1784d496d7641209e94f3b33fe51a8bb7564aa0f
--- /dev/null
+++ b/include/otbSARUpdateMetadataImageFilter.h
@@ -0,0 +1,135 @@
+/*
+ * 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 otbSARUpdateMetadataImageFilter_h
+#define otbSARUpdateMetadataImageFilter_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 "-Wunused-parameter"
+# pragma GCC diagnostic ignored "-Woverloaded-virtual"
+# pragma GCC diagnostic ignored "-Wshadow"
+# include "ossim/ossimTimeUtilities.h"
+
+# pragma GCC diagnostic pop
+
+#else
+# include "ossim/ossimTimeUtilities.h"
+#endif
+
+#if defined(USE_BOOST_TIME)
+# include <boost/date_time/posix_time/posix_time.hpp>
+#include <ostream>
+#endif
+
+namespace otb
+{
+/** \class SARUpdateMetadataImageFilter
+ * \brief Update some metadata after a correction chain.
+ *
+ * This filter updates some metadata (only metadata, input image remains the same).
+ *
+ * \ingroup DiapOTBModule
+ */
+
+ template <typename TImage>
+ class ITK_EXPORT SARUpdateMetadataImageFilter :
+ public itk::ImageToImageFilter<TImage,TImage>
+{
+public:
+
+ // Standard class typedefs
+ typedef SARUpdateMetadataImageFilter Self;
+ typedef itk::ImageToImageFilter<TImage,TImage> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Typedef to image input/output type */
+ typedef TImage ImageType;
+ /** Typedef to describe the inout image pointer type. */
+ typedef typename ImageType::Pointer ImagePointer;
+ typedef typename ImageType::RegionType ImageRegionType;
+
+ // Method for creation through object factory
+ itkNewMacro(Self);
+ // Run-time type information
+ itkTypeMacro(SARUpdateMetadataImageFilter,ImageToImageFilter);
+
+ // Setter
+ itkSetMacro(TimeFirstLine, std::string);
+ itkSetMacro(SlantRange, double);
+
+ // Getter
+ itkGetMacro(TimeFirstLine, std::string);
+ itkGetMacro(SlantRange, double);
+
+
+protected:
+ // Constructor
+ SARUpdateMetadataImageFilter()
+ {
+ m_TimeFirstLine = "";
+ m_SlantRange = 0.;
+ };
+
+ // Destructor
+ virtual ~SARUpdateMetadataImageFilter() {};
+
+ // Print
+ void PrintSelf(std::ostream & os, itk::Indent indent) const ITK_OVERRIDE;
+
+ /** SARUpdateMetadataImageFilter updates the imagekeywordlist for new output metadata.
+ * As such, SARUpdateMetadataImageFilter needs to provide an implementation for
+ * GenerateOutputInformation() in order to inform the pipeline execution model.
+ */
+ virtual void GenerateOutputInformation() ITK_OVERRIDE;
+
+ void ThreadedGenerateData(const ImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
+ ITK_OVERRIDE;
+
+ private:
+ SARUpdateMetadataImageFilter(const Self&); // purposely not implemented
+ void operator=(const Self &); // purposely not
+
+ // Precise metadata
+ std::string m_TimeFirstLine;
+ double m_SlantRange;
+
+};
+
+} // End namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSARUpdateMetadataImageFilter.txx"
+#endif
+
+
+
+#endif
diff --git a/include/otbSARUpdateMetadataImageFilter.txx b/include/otbSARUpdateMetadataImageFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..13bbffbe204bb4d2da2195555902fe767ae19ff0
--- /dev/null
+++ b/include/otbSARUpdateMetadataImageFilter.txx
@@ -0,0 +1,115 @@
+/*
+ * 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 otbSARUpdateMetadataImageFilter_txx
+#define otbSARUpdateMetadataImageFilter_txx
+
+#include "otbSARUpdateMetadataImageFilter.h"
+
+#include "itkImageRegionConstIterator.h"
+#include "itkImageRegionIterator.h"
+
+#include <cmath>
+#include <algorithm>
+
+namespace otb
+{
+ /**
+ * Print
+ */
+ template<class TImage>
+ void
+ SARUpdateMetadataImageFilter< TImage >
+ ::PrintSelf(std::ostream & os, itk::Indent indent) const
+ {
+ Superclass::PrintSelf(os, indent);
+
+ os << indent << "Precise first line time : " << m_TimeFirstLine << std::endl;
+ os << indent << "Precise Slant range : " << m_SlantRange << std::endl;
+ }
+
+
+ /**
+ * Method GenerateOutputInformaton()
+ **/
+ template<class TImage>
+ void
+ SARUpdateMetadataImageFilter< TImage >
+ ::GenerateOutputInformation()
+ {
+ // Call superclass implementation
+ Superclass::GenerateOutputInformation();
+
+ // Get Input and Output ptr
+ ImagePointer outputPtr = this->GetOutput();
+ ImagePointer inputPtr = const_cast< ImageType * >(this->GetInput());
+
+ // Set new keyword list to output image with precise metadata
+ ImageKeywordlist outputKWL = inputPtr->GetImageKeywordlist();
+ //outputKWL.SetMetadataByKey("support_data.first_line_time", m_TimeFirstLine);
+ //outputKWL.SetMetadataByKey("support_data.slant_range_to_first_pixel", m_SlantRange);
+
+ std::ostringstream str_slantRange;
+ str_slantRange.precision(14);
+ str_slantRange << std::fixed << m_SlantRange;
+
+ std::cout << "str_slantRange = " << str_slantRange.str() << std::endl;
+
+ outputKWL.AddKey("support_data.first_line_time", m_TimeFirstLine);
+ outputKWL.AddKey("support_data.slant_range_to_first_pixel", str_slantRange.str());
+ outputPtr->SetImageKeywordList(outputKWL);
+ }
+
+
+ template<class TImage>
+ void
+ SARUpdateMetadataImageFilter< TImage >
+ ::ThreadedGenerateData(const ImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
+ {
+ // Get Input and Output ptr
+ ImagePointer outputPtr = this->GetOutput();
+ ImagePointer inputPtr = const_cast< ImageType * >(this->GetInput());
+
+ // Same region
+ ImageRegionType inputRegion = outputRegionForThread;
+
+ // Iterators
+ typedef itk::ImageRegionConstIterator<ImageType> InItType;
+ typedef itk::ImageRegionIterator<ImageType> OutItType;
+
+ InItType InIt(inputPtr, inputRegion);
+ InIt.GoToBegin();
+
+ OutItType OutIt(outputPtr, outputRegionForThread);
+ OutIt.GoToBegin();
+
+ // Copy input image into output
+ while (!InIt.IsAtEnd())
+ {
+ OutIt.Set(InIt.Get());
+
+ ++OutIt;
+ ++InIt;
+ }
+ }
+
+} /*namespace otb*/
+
+#endif
diff --git a/include/otbTilesExampleFunctor.h b/include/otbTilesExampleFunctor.h
index 4c9f7fe4ced99d5e2994307151596f483adde26f..a7d9f477a386593c6a96fe9ee367b24959e7764b 100644
--- a/include/otbTilesExampleFunctor.h
+++ b/include/otbTilesExampleFunctor.h
@@ -111,7 +111,7 @@ public:
// Process
void Process(ImageInPointer inputPtr, ImageOutRegionType outputRegion, ImageOutPixelType * TilesResult,
- int threadId=0)
+ int /*threadId=0*/)
{
// Transpose outputRegion (current tile) into input image
// Input (= output now)
diff --git a/json_schemas/schema_S1SM.json b/json_schemas/schema_S1SM.json
index 8605109824bb68cc0578e3772dcbe8e5c6b5eac4..e649de737ca723bdd7e2df4b51ef9368f8c5f792 100644
--- a/json_schemas/schema_S1SM.json
+++ b/json_schemas/schema_S1SM.json
@@ -30,6 +30,14 @@
"additionalProperties": false,
"properties": {"output_dir": {"type": "string"}}
}
+ "sensor":
+ {
+ "type": "object",
+ "required": ["satellite", "mode"],
+ "additionalProperties": false,
+ "properties": {"satellite": {"type": "string"},
+ "mode": {"type": "string"}}
+ }
},
"additionalProperties": false,
"required": ["in", "out"]
diff --git a/python_src/diapOTB.py b/python_src/diapOTB.py
index f99ef79afeaab02462dff2f901d682bd2e18b990..e59c71d7d96435155d164c365d9cd20961bc63a6 100644
--- a/python_src/diapOTB.py
+++ b/python_src/diapOTB.py
@@ -35,7 +35,7 @@ from jsonschema import validate
import os
import sys
import argparse
-
+import h5py
import otbApplication as otb
@@ -108,6 +108,13 @@ if __name__ == "__main__":
dem = dict_Global['in']['DEM_Path']
output_dir = dict_Global['out']['output_dir']
+ satellite = "default"
+ mode = "default"
+
+ if 'sensor' in dict_Global:
+ satellite = dict_Global['sensor']['satellite']
+ mode = dict_Global['sensor']['mode']
+
# Pre_Processing
ml_range = dict_PreProcessing['parameter']['ML_range']
ml_azimut = dict_PreProcessing['parameter']['ML_azimut']
@@ -158,7 +165,60 @@ if __name__ == "__main__":
master_Image_base = os.path.basename(master_Image)
slave_Image_base = os.path.basename(slave_Image)
-
+ # Check extension (if .h5 => HDF5 file => Cosmo Sensor)
+ master_ext = master_Image.split(".")[-1:]
+ slave_ext = slave_Image.split(".")[-1:]
+
+ print ("master_ext = " + master_ext[0] + "\n")
+ print ("salve_ext = " + slave_ext[0] + "\n")
+
+ if master_ext[0] == "h5" :
+ master_H5 = h5py.File(master_Image, 'r')
+ lDataSet_master = list(master_H5.keys())
+
+
+ if len(lDataSet_master) != 1 :
+ print("H5 input files does not contain the expected dataset \n")
+ quit()
+
+ if lDataSet_master[0] != "S01" :
+ print("H5 input files does not contain the expected dataset \n")
+ quit()
+
+ master_S01 = dict(master_H5['S01'])
+
+ if not 'SBI' in master_S01:
+ print("H5 input files does not contain the expected dataset \n")
+ quit()
+
+ # Change the name of master and slave image to read directly the //S01/SBI
+ master_Image = "HDF5:" + master_Image + "://S01/SBI"
+ # Adapt sattelite
+ satellite = "cosmo"
+
+
+ if slave_ext[0] == "h5" :
+ slave_H5 = h5py.File(slave_Image, 'r')
+ lDataSet_slave = list(slave_H5.keys())
+
+ if len(lDataSet_slave) != 1 :
+ print("H5 input files does not contain the expected dataset \n")
+ quit()
+
+ if lDataSet_slave[0] != "S01" :
+ print("H5 input files does not contain the expected dataset \n")
+ quit()
+
+ slave_S01 = dict(slave_H5['S01'])
+
+ if not 'SBI' in slave_S01 :
+ print("H5 input files does not contain the expected dataset \n")
+ quit()
+
+ slave_Image = "HDF5:" + slave_Image + "://S01/SBI"
+
+ print("master_Image = " + master_Image + "\n")
+ print("slave_Image = " + slave_Image + "\n")
####################### Pre Processing Chain ##########################
######## SARDoppler Application #######
@@ -241,8 +301,6 @@ if __name__ == "__main__":
appCorGrid.SetParameterInt("mlazi",ml_correlSimu_azimut)
appCorGrid.SetParameterInt("gridsteprange", correlSimu_gridstep_range)
appCorGrid.SetParameterInt("gridstepazimut", correlSimu_gridstep_azimut)
- appCorGrid.SetParameterInt("patchsize", 20)
- appCorGrid.SetParameterInt("subpixel", 10)
appCorGrid.SetParameterString("ram", "4000")
#appCorGrid.Execute()
appCorGrid.ExecuteAndWriteOutput()
@@ -306,11 +364,12 @@ if __name__ == "__main__":
appFineDeformationGrid.SetParameterInt("mlazi", ml_geoGrid_azimut)
appFineDeformationGrid.SetParameterInt("gridsteprange", geoGrid_gridstep_range)
appFineDeformationGrid.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut)
- appFineDeformationGrid.SetParameterInt("patchsize", 20)
- appFineDeformationGrid.SetParameterInt("subpixel", 10)
appFineDeformationGrid.SetParameterFloat("threshold", geoGrid_threshold)
appFineDeformationGrid.SetParameterFloat("gap", geoGrid_gap)
appFineDeformationGrid.SetParameterString("out", os.path.join(output_dir, fine_grid))
+ # Projection is not reliable for cosmo sensor => correction of all values with correlation grid
+ if satellite == "cosmo" or satellite == "CSK" :
+ appFineDeformationGrid.SetParameterString("advantage", "correlation")
appFineDeformationGrid.SetParameterString("ram", "4000")
appFineDeformationGrid.ExecuteAndWriteOutput()
diff --git a/python_src/diapOTB_S1IW.py b/python_src/diapOTB_S1IW.py
index 2dee6f144d5198ff1626bdf374b8570c7ce53429..7639b60c8075572d15fa70b26eadc6cef486d66b 100644
--- a/python_src/diapOTB_S1IW.py
+++ b/python_src/diapOTB_S1IW.py
@@ -449,8 +449,6 @@ if __name__ == "__main__":
appFineDeformationGrid.SetParameterInt("mlazi", ml_geoGrid_azimut)
appFineDeformationGrid.SetParameterInt("gridsteprange", geoGrid_gridstep_range)
appFineDeformationGrid.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut)
- appFineDeformationGrid.SetParameterInt("patchsize", 20)
- appFineDeformationGrid.SetParameterInt("subpixel", 10)
appFineDeformationGrid.SetParameterFloat("threshold", geoGrid_threshold)
appFineDeformationGrid.SetParameterFloat("gap", geoGrid_gap)
appFineDeformationGrid.SetParameterString("out", os.path.join(burst_dir, fine_grid))
diff --git a/test/otbSARDopplerCentroidFreqImageFilterTest.cxx b/test/otbSARDopplerCentroidFreqImageFilterTest.cxx
index a983cdd2413b6d61ebb4e963bb502c7cf9f10511..cb2844db5724119a16eaf49d20794ee2f2b30b2d 100644
--- a/test/otbSARDopplerCentroidFreqImageFilterTest.cxx
+++ b/test/otbSARDopplerCentroidFreqImageFilterTest.cxx
@@ -39,8 +39,7 @@ int otbSARDopplerCentroidFreqFilterTest(int argc, char * argv[])
// Image type
const unsigned int Dimension = 2;
- typedef otb::Image <std::complex<float>, 2> ImageSARInType ;
- typedef otb::Image <float, 2> FloatImageType;
+ typedef otb::Image <float, Dimension> FloatImageType;
typedef otb::VectorImage<float> FloatVectorImageType;
// Reader