otbSARESD.cxx

/*
 * Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"

#include "otbMultiToMonoChannelExtractROI.h"
#include "otbImageList.h"
#include "otbImageListToVectorImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "otbSARStreamingMeanPhaseAndAzimutShiftImageFilter.h"
#include "otbSARESDOffsetsImageFilter.h"
#include "otbSARDopplerCentroidFreqImageFilter.h"

#include <iostream>
#include <string>
#include <fstream>

namespace otb
{
namespace Wrapper
{

class SARESD : public Application
{
public:
  typedef SARESD Self;
  typedef itk::SmartPointer<Self> Pointer; 

  itkNewMacro(Self);
  itkTypeMacro(SARESD, otb::Wrapper::Application);

  /** Filters typedef */
  typedef otb::SARStreamingMeanPhaseAndAzimutShiftImageFilter<FloatVectorImageType>   MeanFilter;;
  typedef otb::ImageList<FloatImageType>                                              ImageListType;
  typedef otb::ImageListToVectorImageFilter<ImageListType,
                                       FloatVectorImageType >                         ListConcatenerFilterType;
  typedef otb::MultiToMonoChannelExtractROI<FloatVectorImageType::InternalPixelType,
					    FloatImageType::PixelType>                ExtractROIFilterType;
  typedef itk::SubtractImageFilter <FloatImageType, FloatImageType, FloatImageType>   SubFilterType;

  typedef ObjectList<ExtractROIFilterType>                                            ExtractROIFilterListType;

  typedef otb::SARDopplerCentroidFreqImageFilter<FloatImageType>                      DCFFilterType;
  typedef otb::SARESDOffsetsImageFilter<FloatImageType>                               ESDFilterType;
  

private:
  void DoInit() override
  {
    SetName("SARESD");
    SetDescription("ESD processing to correct phase between two bursts.");

    SetDocName("SAR ESD Offset");
    SetDocLongDescription("This application applies esd processing on two consecutives bursts"
      " dimension.");

    //Optional descriptors
    SetDocLimitations("Only Sentinel 1 products are supported for now.");
    SetDocAuthors("OTB-Team");
    SetDocSeeAlso(" ");
    AddDocTag(Tags::SAR);

    //Parameter declarations
    AddParameter(ParameterType_InputImage,  "ininterfup",   "First interferogram (Vector Image)");
    SetParameterDescription("ininterfup", "First interferogram.");

    AddParameter(ParameterType_InputImage,  "ininterflow",   "Second interferogram (Vector Image)");
    SetParameterDescription("ininterflow", "Second interferogram.");

    AddParameter(ParameterType_ComplexInputImage,  "insar",   "Input SAR  image (for metadata only");
    SetParameterDescription("insar", "SAR Image to extract SAR geometry.");

    AddParameter(ParameterType_Int,  "burstindex", "Index of first interferogram");
    SetParameterDescription("burstindex", "Index of first interferogram.");
    SetDefaultParameterInt("burstindex", 0);

    AddParameter(ParameterType_Float, "threshold", "Threshold for correlation rate");
    SetParameterDescription("threshold","Threshold for correlation rate.");
    SetDefaultParameterFloat("threshold", 0.3);
    SetMinimumParameterFloatValue("threshold", 0.0001);
    MandatoryOff("threshold");

    AddParameter(ParameterType_Int, "mlazi", "Averaging on azimut");
    SetParameterDescription("mlazi","Averaging on azimut");
    SetDefaultParameterInt("mlazi", 1);
    SetMinimumParameterIntValue("mlazi", 1);
    MandatoryOff("mlazi");
      
    AddParameter(ParameterType_OutputImage, "out", "Output Tmp");
    SetParameterDescription("out","Output Tmp.");

    AddParameter(ParameterType_Float,"azishift","Azimut shift to correct phase jump into two consecutive interferograms");
    SetParameterDescription("azishift", "Azimut shift to correct phase jump.");
    SetParameterRole("azishift", Role_Output);
    SetDefaultParameterFloat("azishift", 0);
 
    AddRAMParameter();

    SetDocExampleParameterValue("ininterfup","./interferogram_burst0.tiff");
    SetDocExampleParameterValue("ininterflow", "./interferogram_burst1.tiff");
    SetDocExampleParameterValue("insar","s1a-iw1-slc-vv-20150821t152433-20150821t152458-007363-00a1e0-004.tiff");
    SetDocExampleParameterValue("burstindex", "0");
    SetDocExampleParameterValue("out","phaseDiff.tif");
  }

  void DoUpdateParameters() override
  {
    // Nothing // TODO: o do here : all parameters are independent
  }

void DoExecute() override
{  
  ////////// Instanciate all filters //////////
  ExtractROIFilterType::Pointer extractorUp = ExtractROIFilterType::New();
  ExtractROIFilterType::Pointer extractorLow = ExtractROIFilterType::New();

  ExtractROIFilterType::Pointer extractorCohUp = ExtractROIFilterType::New();
  ExtractROIFilterType::Pointer extractorCohLow = ExtractROIFilterType::New();

  ExtractROIFilterType::Pointer extractorIsDataUp = ExtractROIFilterType::New();
  ExtractROIFilterType::Pointer extractorIsDataLow = ExtractROIFilterType::New();

  ExtractROIFilterType::Pointer extractorDcfUp = ExtractROIFilterType::New();
  ExtractROIFilterType::Pointer extractorDcfLow = ExtractROIFilterType::New();

  MeanFilter::Pointer meanPhase = MeanFilter::New();
  meanPhase->SetPhaseMean(true);
  
  MeanFilter::Pointer meanAziShift = MeanFilter::New();
  meanAziShift->SetPhaseMean(false);
    
  SubFilterType::Pointer subPhase = SubFilterType::New();

  ListConcatenerFilterType::Pointer m_Concatener = ListConcatenerFilterType::New();
  ImageListType::Pointer m_ImageList = ImageListType::New();

  ListConcatenerFilterType::Pointer m_Concatener1 = ListConcatenerFilterType::New();
  ImageListType::Pointer m_ImageList1 = ImageListType::New();
  
  DCFFilterType::Pointer dcfUp = DCFFilterType::New();
  DCFFilterType::Pointer dcfLow = DCFFilterType::New();

  ESDFilterType::Pointer esdFilter = ESDFilterType::New();
  
  ////////// Get and Check inputs //////////
  // Get Number of burst index
  float threshold = GetParameterFloat("threshold");
  int ml_azi = GetParameterInt("mlazi");
  int burstindex = GetParameterInt("burstindex");
  otbAppLogINFO(<< "burst index : " << burstindex);
  otbAppLogINFO(<< "ML factor for azimut dimension : " << ml_azi);
  otbAppLogINFO(<< "Threshold for phase difference : " << threshold);

  // Get input interferograms
  FloatVectorImageType::Pointer InterfUpPtr = GetParameterFloatVectorImage("ininterfup");
  FloatVectorImageType::Pointer InterfLowPtr = GetParameterFloatVectorImage("ininterflow");

  // Check number of components (at least three : modulus, phase and coherency)
  int nbComponentsUp = InterfUpPtr->GetNumberOfComponentsPerPixel();
  int nbComponentsLow = InterfLowPtr->GetNumberOfComponentsPerPixel();
  if (nbComponentsUp < 3 || nbComponentsLow < 3 || nbComponentsUp != nbComponentsLow)
    {
      otbAppLogFATAL(<< "Input interferograms have to contain at least three components");
    }  

  // Check invalid Pixel Key
  const bool inputWithInvalidPixels1 = InterfUpPtr->GetImageKeywordlist().HasKey("support_data.invalid_pixels")
    && InterfUpPtr->GetImageKeywordlist().GetMetadataByKey("support_data.invalid_pixels") == "yes";
  const bool inputWithInvalidPixels2 = InterfLowPtr->GetImageKeywordlist().HasKey("support_data.invalid_pixels")
    && InterfLowPtr->GetImageKeywordlist().GetMetadataByKey("support_data.invalid_pixels") == "yes";
  
  if (inputWithInvalidPixels1!= inputWithInvalidPixels2)
    {
      // Throw an execption
      otbAppLogFATAL(<< "Incoherency between input images (for support_data.invalid_pixels key).");
    }
  

  ////////// Determine overlap area between the two bursts //////////
  std::pair<unsigned long,unsigned long> linesUp;
  std::pair<unsigned long,unsigned long> linesLow;
  std::pair<unsigned long,unsigned long> samplesUp;
  std::pair<unsigned long,unsigned long> samplesLow;
  meanPhase->getOverlapLinesAndSamples(GetParameterComplexFloatImage("insar")->GetImageKeywordlist(),
					    linesUp, linesLow, samplesUp, samplesLow, burstindex, 
					    inputWithInvalidPixels1);

  
  // Check origin (msut have the whole burst/interferogram)
  InterfUpPtr->UpdateOutputInformation();
  InterfLowPtr->UpdateOutputInformation();
  unsigned long originOffset_samplesUp = static_cast<long>(InterfUpPtr->GetOrigin()[0]-0.5);
  unsigned long originOffset_linesUp = static_cast<long>(InterfUpPtr->GetOrigin()[1]-0.5);
  unsigned long originOffset_samplesLow = static_cast<long>(InterfLowPtr->GetOrigin()[0]-0.5);
  unsigned long originOffset_linesLow = static_cast<long>(InterfLowPtr->GetOrigin()[1]-0.5);

  if (originOffset_samplesLow != 0 || originOffset_linesLow != 0)
    {
      otbAppLogFATAL(<< "Input interferogram (Low) have not the expected origin (0,0)");
    }

  if (originOffset_samplesUp != 0 || originOffset_linesUp != 0)
    {
      otbAppLogFATAL(<< "Input interferogram (Up) have not the expected origin (0,0)");
    }

  // Retrieve start and size for each burst
  unsigned long minSamples_Up = std::min(samplesUp.second, 
					 static_cast<unsigned long>(InterfUpPtr->GetLargestPossibleRegion().GetSize()[0]-1));
  unsigned long minLines_Up = std::min(linesUp.second, 
				       static_cast<unsigned long>(InterfUpPtr->GetLargestPossibleRegion().GetSize()[1]-1));

  unsigned long minSamples_Low = std::min(samplesLow.second, 
					 static_cast<unsigned long>(InterfLowPtr->GetLargestPossibleRegion().GetSize()[0]-1));
  unsigned long minLines_Low = std::min(linesLow.second, 
				       static_cast<unsigned long>(InterfLowPtr->GetLargestPossibleRegion().GetSize()[1]-1));

  unsigned long startL_Up = linesUp.first;
  unsigned long sizeL_Up = minLines_Up - linesUp.first + 1;
  unsigned long startS_Up = samplesUp.first;
  unsigned long sizeS_Up = minSamples_Up - samplesUp.first + 1;

  unsigned long startL_Low = linesLow.first;
  unsigned long sizeL_Low = minLines_Low - linesLow.first + 1;
  unsigned long startS_Low = samplesLow.first;
  unsigned long sizeS_Low = minSamples_Low - samplesLow.first + 1;

  std::cout << "startL_Up = " << startL_Up << std::endl;
  std::cout << "sizeL_Up = " << sizeL_Up << std::endl;
  std::cout << "startS_Up = " << startS_Up << std::endl;
  std::cout << "sizeS_Up = " << sizeS_Up << std::endl;

  std::cout << "startL_Low = " << startL_Low << std::endl;
  std::cout << "sizeL_Low = " << sizeL_Low << std::endl;
  std::cout << "startS_Low = " << startS_Low << std::endl;
  std::cout << "sizeS_Low = " << sizeS_Low << std::endl;
  
  std::cout << "samplesLow.first = " << samplesLow.first << std::endl;
  std::cout << "samplesUp.first = " << samplesUp.first << std::endl;
  std::cout << "samplesLow.second = " << samplesLow.second << std::endl;
  std::cout << "samplesUp.second = " << samplesUp.second << std::endl;
  std::cout << "minSamples_Low = " << minSamples_Low << std::endl;
  std::cout << "minSamples_Up = " << minSamples_Up << std::endl;

  ////////// Extract Phase for overlap area //////////
  // Set the channel to extract : Phase = 2
  extractorUp->SetInput(InterfUpPtr);
  extractorUp->SetChannel(2);      
  extractorUp->SetStartX(startS_Up);
  extractorUp->SetStartY(startL_Up);
  extractorUp->SetSizeX(sizeS_Up);
  extractorUp->SetSizeY(sizeL_Up);
  extractorUp->UpdateOutputInformation();

  extractorLow->SetInput(InterfLowPtr);
  extractorLow->SetChannel(2);    
  extractorLow->SetStartX(startS_Low);
  extractorLow->SetStartY(startL_Low);
  extractorLow->SetSizeX(sizeS_Low);
  extractorLow->SetSizeY(sizeL_Low);
  extractorLow->UpdateOutputInformation();

  //////////// Subtract phase ////////////
  subPhase->SetCoordinateTolerance(10e08);
  subPhase->SetInput1(extractorUp->GetOutput());
  subPhase->SetInput2(extractorLow->GetOutput());

  //////////// Extract Coherency (Coherency = 3) ////////////
  extractorCohUp->SetInput(InterfUpPtr);
  extractorCohUp->SetChannel(3);      
  extractorCohUp->SetStartX(startS_Up);
  extractorCohUp->SetStartY(startL_Up);
  extractorCohUp->SetSizeX(sizeS_Up);
  extractorCohUp->SetSizeY(sizeL_Up);
  extractorCohUp->UpdateOutputInformation();

  extractorCohLow->SetInput(InterfLowPtr);
  extractorCohLow->SetChannel(3);    
  extractorCohLow->SetStartX(startS_Low);
  extractorCohLow->SetStartY(startL_Low);
  extractorCohLow->SetSizeX(sizeS_Low);
  extractorCohLow->SetSizeY(sizeL_Low);
  extractorCohLow->UpdateOutputInformation();

  //////////// Extract isData (isData Channel = 4) ////////////
  extractorIsDataUp->SetInput(InterfUpPtr);
  extractorIsDataUp->SetChannel(4);      
  extractorIsDataUp->SetStartX(startS_Up);
  extractorIsDataUp->SetStartY(startL_Up);
  extractorIsDataUp->SetSizeX(sizeS_Up);
  extractorIsDataUp->SetSizeY(sizeL_Up);
  extractorIsDataUp->UpdateOutputInformation();

  extractorIsDataLow->SetInput(InterfLowPtr);
  extractorIsDataLow->SetChannel(4);    
  extractorIsDataLow->SetStartX(startS_Low);
  extractorIsDataLow->SetStartY(startL_Low);
  extractorIsDataLow->SetSizeX(sizeS_Low);
  extractorIsDataLow->SetSizeY(sizeL_Low);
  extractorIsDataLow->UpdateOutputInformation();
  
  //////////// Concatenation of phase difference with coherency ////////////
  m_ImageList->PushBack( subPhase->GetOutput() );
  m_ImageList->PushBack( extractorCohUp->GetOutput() );
  m_ImageList->PushBack( extractorCohLow->GetOutput() );
  m_ImageList->PushBack( extractorIsDataUp->GetOutput() );
  m_ImageList->PushBack( extractorIsDataLow->GetOutput() );
  
  m_Concatener->SetInput( m_ImageList );
  
  ////////// Mean of phase difference //////////
  meanPhase->SetInput(m_Concatener->GetOutput());
  meanPhase->SetThreshold(threshold);
  meanPhase->Update();

  std::cout << "meanPhase->GetMean() = " << meanPhase->GetMean() << std::endl;

  // Doppler Centroid Filter
  dcfUp->SetInput(extractorUp->GetOutput());
  dcfLow->SetInput(extractorLow->GetOutput());

  // Extract overlap area into Doppler Centroid outputs
  // extractorDcfUp->SetInput(dcfUp->GetOutput());
  // extractorDcfUp->SetChannel(1);      
  // extractorDcfUp->SetStartX(startS_Up);
  // extractorDcfUp->SetStartY(startL_Up);
  // extractorDcfUp->SetSizeX(sizeS_Up);
  // extractorDcfUp->SetSizeY(sizeL_Up);
  // extractorDcfUp->UpdateOutputInformation();

  // extractorDcfLow->SetInput(dcfLow->GetOutput());
  // extractorDcfLow->SetChannel(1);    
  // extractorDcfLow->SetStartX(startS_Low);
  // extractorDcfLow->SetStartY(startL_Low);
  // extractorDcfLow->SetSizeX(sizeS_Low);
  // extractorDcfLow->SetSizeY(sizeL_Low);
  // extractorDcfLow->UpdateOutputInformation();

  // //////////// ESD Filter ////////////
  esdFilter->SetPhaseDiffMean(meanPhase->GetMean());
  double aziInterval = std::stod(InterfUpPtr->GetImageKeywordlist().GetMetadataByKey("support_data.line_time_interval"));
  esdFilter->SetAziTimeInt(aziInterval);
  //esdFilter->SetInput(subPhase->GetOutput());
  esdFilter->SetInput1(dcfUp->GetOutput());
  esdFilter->SetInput2(dcfLow->GetOutput());
  esdFilter->SetCoordinateTolerance(10e08);

  //////////// Concatenation with coherency ////////////
  m_ImageList1->PushBack( esdFilter->GetOutput() );
  m_ImageList1->PushBack( extractorCohUp->GetOutput() );
  m_ImageList1->PushBack( extractorCohLow->GetOutput() );
  
  m_Concatener1->SetInput( m_ImageList1 );

  ////////// Mean of azimut shift ///////////
  meanAziShift->SetInput(m_Concatener1->GetOutput());
  meanAziShift->Update();

  std::cout << "meanAziShift->GetMean() = " << meanAziShift->GetMean() << std::endl;

  ////////// Outputs : Image of azimut shift (for overlap area) and mean of azimut shift ////////////
  //SetParameterOutputImage("out", dcfLow->GetOutput());
  SetParameterOutputImage("out", subPhase->GetOutput());
  //SetParameterOutputImage("out", esdFilter->GetOutput());
  SetParameterFloat("azishift", meanAziShift->GetMean());
  RegisterPipeline();
}
   // Vector for filters 
  std::vector<itk::ProcessObject::Pointer> m_Ref;
};


}

}

OTB_APPLICATION_EXPORT(otb::Wrapper::SARESD)