ENH : New application SARDeramp (part 1)

include/otbSARDerampImageFilter.h

+/*
+ * Copyright (C) 2005-2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbSARDerampImageFilter_h
+#define otbSARDerampImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkSmartPointer.h"
+#include "itkPoint.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageScanlineIterator.h"
+
+#include "otbImageKeywordlist.h"
+#include "otbSarSensorModelAdapter.h"
+
+namespace otb
+{
+/** \class SARDerampImageFilter 
+ * \brief Deramps or Reramps a burst from Sentinel-1 IW product. 
+ * 
+ * This filter deramps or reramps an input burst.
+ *
+ * \ingroup DiapOTBModule
+ */
+
+  template <typename TImage> 
+  class ITK_EXPORT SARDerampImageFilter :
+    public itk::ImageToImageFilter<TImage,TImage>
+{
+public:
+
+  // Standard class typedefs
+  typedef SARDerampImageFilter                    Self;
+  typedef itk::ImageToImageFilter<TImage,TImage>         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(SARDerampImageFilter,ImageToImageFilter);
+
+  /** Typedef to image input/output type VectorImage (Complex Image)   */
+  typedef TImage                                  ImageInType;
+  /** Typedef to describe the inout image pointer type. */
+  typedef typename ImageType::Pointer             ImagePointer;
+  typedef typename ImageType::ConstPointer        ImageConstPointer;
+  /** Typedef to describe the inout image region type. */
+  typedef typename ImageType::RegionType          ImageRegionType;
+  /** Typedef to describe the type of pixel and point for inout image. */
+  typedef typename ImageType::PixelType           ImagePixelType;
+  typedef typename ImageType::PointType           ImagePointType;
+  /** Typedef to describe the image index, size types and spacing for inout image. */
+  typedef typename ImageType::IndexType           ImageIndexType;
+  typedef typename ImageType::IndexValueType      ImageIndexValueType;
+  typedef typename ImageType::SizeType            ImageSizeType;
+  typedef typename ImageType::SizeValueType       ImageSizeValueType;
+  typedef typename ImageType::SpacingType         ImageSpacingType;
+  typedef typename ImageType::SpacingValueType    ImageSpacingValueType;
+
+  // Define Point2DType and Point3DType
+  using Point2DType = itk::Point<double,2>;
+  using Point3DType = itk::Point<double,3>;
+
+  // Typedef for iterators
+  typedef itk::ImageScanlineConstIterator< ImageType > InputIterator;
+  typedef itk::ImageScanlineIterator< ImageType > OutputIterator;
+ 
+  // Setter
+  itkSetMacro(DerampMode, bool);
+
+  // Getter
+  itkGetMacro(DerampMode, bool);
+
+protected:
+  // Constructor
+  SARDerampImageFilter();
+
+  // Destructor
+  virtual ~SARDerampImageFilter() ITK_OVERRIDE;
+
+  // Print
+  void PrintSelf(std::ostream & os, itk::Indent indent) const ITK_OVERRIDE;
+  
+  /** 
+   * SARDerampImageFilter can be implemented as a multithreaded filter.
+   * Therefore, this implementation provides a ThreadedGenerateData() routine
+   * which is called for each processing thread. The main output image data is
+   * allocated automatically by the superclass prior to calling
+   * ThreadedGenerateData().  ThreadedGenerateData can only write to the
+   * portion of the output image specified by the parameter
+   * "outputRegionForThread"
+   *
+   * \sa ImageToImageFilter::ThreadedGenerateData(),
+   *     ImageToImageFilter::GenerateData() */
+  virtual void ThreadedGenerateData(const ImageRegionType& outputRegionForThread, 
+				    itk::ThreadIdType threadId) ITK_OVERRIDE;
+
+  /** 
+   * SARDerampImageFilter reuses calculations into ThreadedGeneratedData. The aim is to estimate once and for all
+   * some calculations and store results into argument of this class.
+   *
+   * \sa ImageToImageFilter::BeforeThreadedGenerateData()*/
+  void BeforeThreadedGenerateData() ITK_OVERRIDE;
+  
+
+ private:
+  SARDerampImageFilter(const Self&); // purposely not implemented
+  void operator=(const Self &); // purposely not 
+
+  // Polynomial Selection
+  bool selectFMRateCoef();
+  bool selectDCFCoef();
+
+  // Deramp (if true) or Reramp (if false) mode
+  bool m_DerampMode;
+
+  // First azimuth/range time
+  float m_FirstAziTime, m_FirstRangeTime;
+
+  // Mid burst zero Doppler azimuth time [s] 
+  float m_MidAziTime;
+  
+  // Spacecraft velocity computed at mid-burst time (m_MidAziTime) [m/s]
+  float m_VSatAtMidAziTime;
+
+  // Doppler Centroid rate introduced by the scanning og the antenna [Hz/s]
+  float m_Ks;
+  
+  // Azimuth time interval [s]
+  float m_AziTimeInt;
+
+  // Range sampling rate [Hz]
+  float m_RangeSamplingRate;
+
+  // Parameters coefficients for Doppler FM Rate
+  float m_FM_C0, m_FM_C1, m_FM_C2, m_FM_Tau0;
+
+  // Parameters coefficients for Doppler Centroid Frequency
+  float m_DCF_C0, m_DCF_C1, m_DCF_C2, m_DCF_Tau0;
+
+  // Reference time 0 [s]
+  float m_RefTime0;
+
+  // First estimation
+  bool m_FirstEstimation;
+
+  const double C;
+};
+
+} // End namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSARDerampImageFilter.txx"
+#endif
+
+
+
+#endif

include/otbSARDerampImageFilter.txx

+/*
+ * 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 otbSARDerampImageFilter_txx
+#define otbSARDerampImageFilter_txx
+
+#include "otbSARDerampImageFilter.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageScanlineIterator.h"
+#include "itkProgressReporter.h"
+#include "itkNumericTraitsPointPixel.h"
+
+#include "otbSarSensorModelAdapter.h"
+
+#include <cmath>
+#include <algorithm>
+#include <omp.h>
+
+namespace otb
+{
+  /** 
+   * Constructor with default initialization
+   */
+  template <class TImage> 
+  SARDerampImageFilter< TImage >::SARDerampImageFilter()
+    :  m_DerampMode(true),  m_FirstAziTime(0), m_FirstRangeTime(0), m_MidAziTime(0), m_VSatAtMidAziTime(0),
+       m_Ks(0), m_AziTimeInt(0),m_RangeSamplingRate(0), m_FM_C0(0), m_FM_C1(0), m_FM_C2(0), m_FM_Tau0(0), 
+       m_DCF_C0(0), m_DCF_C1(0), m_DCF_C2(0), m_DCF_Tau0(0), m_RefTime0(0), m_FirstEstimation(false)
+  {
+C = 299792458;
+  }
+    
+  /** 
+   * Destructor
+   */
+  template <class TImage> 
+  SARDerampImageFilter< TImage >::~SARDerampImageFilter()
+  {
+  }
+
+  /**
+   * Print
+   */
+  template<class TImage>
+  void
+  SARDerampImageFilter< TImage >
+  ::PrintSelf(std::ostream & os, itk::Indent indent) const
+  {
+    Superclass::PrintSelf(os, indent);
+
+    os << indent << "Mode for deramp filter, if true then Deramp Mode else Reramp mode : " << m_DerampMode 
+       << std::endl;
+  }
+
+  /**
+   * Method selectFMRateCoef
+   */
+  template<class TImage>
+  bool 
+  SARDerampImageFilter< TImage >
+  ::selectFMRateCoef()
+  {
+    
+  }
+
+ /**
+   * Method selectDCFCoef
+   */
+  template<class TImage>
+  bool 
+  SARDerampImageFilter< TImage >
+  ::selectDCFCoef()
+  {
+    
+  }
+
+  /**
+   * Method ThreadedGenerateData
+   */
+  template<class TImage>
+  void
+  SARDerampImageFilter<TImage>
+  ::BeforeThreadedGenerateData()
+  {
+    // Estimates general parameters for the current burst, if m_FirstEstimation == true
+    if (m_FirstEstimation)
+      {
+	m_FirstEstimation = false;
+
+	// Get input
+	ImagePointer inputPtr = const_cast< ImageType * >( this->GetInput() );
+	// KeyWordList
+	ImageKeywordlist inputKWL = inputPtr->GetImageKeywordlist();
+
+	// Check version of header/kwl (at least 3)
+	int headerVersion = std::stoi(inputKLW.GetMetadataByKey("header.version"));
+
+	if (headerVersion < 3)
+	  {
+	    itkExceptionMacro(<<"Header version is < 3. Please Upgrade your geom file");
+	  }
+
+	// Get some metadata
+	float aziSteeringRate = std::stof(inputKLW.GetMetadataByKey("support_data.azimuth_steering_rate"));
+	// Conversion to radians per seconds
+	aziSteeringRate *= (M_PI/180); 
+
+	m_FirstAziTime = std::stof(inputKLW.GetMetadataByKey("support_data.first_line_time"));
+	m_FirstRangeTime = std::stof(inputKLW.GetMetadataByKey("support_data.slant_range_to_first_pixel"));
+
+	m_AziTimeInt = std::stof(inputKLW.GetMetadataByKey("support_data.line_time_interval"));
+	m_RangeTimeInt = std::stof(inputKLW.GetMetadataByKey("support_data.range_sampling_rate"));
+
+	float radarFrequency = std::stof(inputKLW.GetMetadataByKey("support_data.radar_frequency"));
+	
+	int nbLineBurst = std::stof(inputKLW.GetMetadataByKey("support_data.geom.bursts.number_lines_per_burst"));
+	// Estimation m_Ks
+	double lineAtMidBurst = nbLineBurst/2.;
+	m_MidAziTime = m_FirstAziTime + m_AziTimeInt * lineAtMidBurst;
+	
+	// Try to create a SarSensorModelAdapter
+	SarSensorModelAdapter::Pointer sarSensorModel = SarSensorModelAdapter::New();
+	bool loadOk = sarSensorModel->LoadState(inputKWL);
+
+	if(!loadOk || !sarSensorModel->IsValidSensorModel())
+	  itkExceptionMacro(<<"Input image does not contain a valid SAR sensor model.");
+
+	Point3DType satpos, satvel;
+ 
+	bool lineToSatPosAndVelOk = sarSensorModel->LineToSatPositionAndVelocity(lineAtMidBurst, satpos, satvel);
+
+	if (!lineToSatPosAndVelOk)
+	  itkExceptionMacro(<<"Failed to estimate satellite position and velocity.");
+  
+	m_VSatAtMidAziTime = std::sqrt(satvel[0]*satvel[0] + satvel[1]*satvel[1] + satvel[2]*satvel[2]);
+
+	m_Ks = (2*m_VSatAtMidAziTime/C) * radarFrequency * aziSteeringRate;
+	
+	// Polynomial selection (FM Rate and Doppler Centroid Frequency)
+	
+      }
+  }
+   
+
+  /**
+   * Method ThreadedGenerateData
+   */
+  template<class TImage>
+  void
+  SARDerampImageFilter<TImage>
+  ::ThreadedGenerateData(const ImageOutRegionType & outputRegionForThread,
+			 itk::ThreadIdType /*threadId*/)
+  {
+    // Compute corresponding input region for master and slave cartesian mean
+    ImageInRegionType inputMasterRegionForThread = outputRegionForThread;
+    
+    // Compute corresponding input region for grid
+    GridRegionType inputGridRegionForThread = OutputRegionToInputGridRegion(outputRegionForThread);
+
+    // Iterator on output
+    OutputIterator OutIt(this->GetOutput(), outputRegionForThread);
+    OutIt.GoToBegin();
+
+    // Iterator on input master cartesian mean
+    InputIterator  InMasterCartMeanIt(this->GetMasterCartesianMeanInput(), inputMasterRegionForThread);
+    InMasterCartMeanIt.GoToBegin();
+
+    // Allocate output pixel
+    ImageOutPixelType pixelOut;
+    if (m_MasterCopy)
+      {
+	pixelOut.Reserve(5);
+      }
+    else
+      {
+	pixelOut.Reserve(2);
+      }
+
+    double constMul = static_cast<double>(m_Factor*2*M_PI)/m_Lambda;
+    
+    if (m_ApproxDiapason)
+      {
+	constMul = static_cast<double>(m_Factor*256)/m_Lambda;
+      }
+
+    Point3DType worldSlave;
+    Point3DType satPosSlave;
+    Point3DType satVelSlave;
+    Point3DType worldMaster;
+    Point3DType satPosMaster;
+    Point3DType satVelMaster;
+
+    // For each line
+    while (!OutIt.IsAtEnd() && !InMasterCartMeanIt.IsAtEnd())
+      {
+	OutIt.GoToBeginOfLine();
+	InMasterCartMeanIt.GoToBeginOfLine();
+
+	// Index of current line (into output Geometry)
+	int ind_Line = OutIt.GetIndex()[1] + int(this->GetOutput()->GetOrigin()[1]);
+
+	// Get Master Cartesian Mean Per line
+	ImageInIndexType indexMasterPerLine;
+	indexMasterPerLine[0] = 0; // Always 0 since Master/Slave Cartesain Per Line are vectors
+	indexMasterPerLine[1] = ind_Line;
+
+
+	// Get the index of current tile into grid to retrive the shifts (the closest (round) grid point 
+	// at the center of current tile). Output Geo = Master Cart Mean Geo = (Grid geo / GridStep) 
+	int Lgrid =  std::round( ind_Line / m_GridStep[1]); 
+	    
+	Lgrid = std::min (std::max (Lgrid, 0), 
+			  static_cast<int>(this->GetGridInput()->GetLargestPossibleRegion().GetSize()[1])-1); 
+	    
+	GridIndexType gridIndex;
+	gridIndex[0] = 0;	    
+	gridIndex[1] = Lgrid;
+	    
+	double gridShift_Azi = this->GetGridInput()->GetPixel(gridIndex)[1];
+
+	// Apply on slave, the integer shifts
+	int Le = std::round(ind_Line + gridShift_Azi); 
+	
+	// Get Slave Cartesian Means 
+	ImageInIndexType indexSlavePerLine;
+	indexSlavePerLine[0] = 0; // Always 0 since Master/Slave Cartesain Per Line are vectors
+	indexSlavePerLine[1] = Le;
+
+	////////// Estimate satellite positions for the current line //////////
+	bool checkSlave = m_SarSensorModelAdapterForSlave->LineToSatPositionAndVelocity(static_cast<double>(Le), satPosSlave, 
+								      satVelSlave);
+
+	bool checkMaster = m_SarSensorModelAdapterForMaster->LineToSatPositionAndVelocity(static_cast<double>(ind_Line), 
+								      satPosMaster, satVelMaster);
+
+	// For each colunm
+	while (!OutIt.IsAtEndOfLine() && !InMasterCartMeanIt.IsAtEndOfLine()) 
+	  {
+	    // Check slave Index
+	     if (checkSlave)
+	      {
+		// Check if Value into Master Cartesian Mean with IsData Mask
+		if (InMasterCartMeanIt.Get()[3] != 0)
+		  {
+		    float Xcart_Ground = InMasterCartMeanIt.Get()[0];
+		    float Ycart_Ground = InMasterCartMeanIt.Get()[1];
+		    float Zcart_Ground = InMasterCartMeanIt.Get()[2];
+
+		    //////////// Estimate Topographic phase (P = (factor*256/lambda) * (De-Dm))  //////////
+		    double De = sqrt(pow((Xcart_Ground - satPosSlave[0]), 2) + 
+				     pow((Ycart_Ground - satPosSlave[1]), 2) + 
+				     pow((Zcart_Ground - satPosSlave[2]), 2));
+		
+		    double Dm = sqrt(pow((Xcart_Ground - satPosMaster[0]), 2) + 
+				     pow((Ycart_Ground - satPosMaster[1]), 2) + 
+				     pow((Zcart_Ground - satPosMaster[2]), 2));
+		    
+		    
+		    pixelOut[0] = constMul * (De-Dm);
+		    // Mod 2*Pi
+		    //pixelOut[0] =  pixelOut[0]-(2*M_PI)*floor(pixelOut[0]/(2*M_PI));
+
+		    // IsData set to 1
+		    pixelOut[1] = 1;
+		    
+		    if (m_MasterCopy)
+		      {
+			////////////// Copy of Master Cartesian Mean  //////////////
+			pixelOut[2] = Xcart_Ground;
+			pixelOut[3] = Ycart_Ground;
+			pixelOut[4] = Zcart_Ground;
+		      }
+		    
+		  }
+		else
+		  {
+		    // All components set to 0
+		    pixelOut[0] = 0;
+		    pixelOut[1] = 0;
+		    if (m_MasterCopy) 
+		      {
+			pixelOut[2] = 0;
+			pixelOut[3] = 0;
+			pixelOut[4] = 0;
+		      }
+		  }
+
+	      }
+	    else
+	      {
+		// All components set to 0
+		pixelOut[0] = 0;
+		pixelOut[1] = 0;
+		if (m_MasterCopy) 
+		  {
+		    pixelOut[2] = 0;
+		    pixelOut[3] = 0;
+		    pixelOut[4] = 0;
+		  }
+		}
+
+	    //////////// Assign Output ////////////
+	    OutIt.Set(pixelOut);
+
+	    // Increment iterators
+	    ++OutIt;
+	    ++InMasterCartMeanIt;
+	  } // End colunms (ouput)
+
+	// Next Line
+	OutIt.NextLine();
+	InMasterCartMeanIt.NextLine();
+      } // End lines (ouput)
+  }
+
+
+} /*namespace otb*/
+
+#endif