ENH : New application SARDeramp (part 2)

app/CMakeLists.txt

@@ -83,3 +83,8 @@ OTB_CREATE_APPLICATION(NAME SARGridStatistics
                        SOURCES otbSARGridStatistics.cxx
                        LINK_LIBRARIES ${${otb-module}_LIBRARIES}
 )
+
+OTB_CREATE_APPLICATION(NAME SARDeramp
+                       SOURCES otbSARDeramp.cxx
+                       LINK_LIBRARIES ${${otb-module}_LIBRARIES}
+)

app/otbSARDeramp.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 "otbSARDerampImageFilter.h"
+
+#include "otbWrapperOutputImageParameter.h"
+#include "otbWrapperTypes.h"
+
+#include <iostream>
+#include <string>
+#include <fstream>
+
+namespace otb
+{
+
+  namespace Wrapper
+  {
+
+    class SARDeramp : public Application
+    {
+    public:
+      typedef SARDeramp Self;
+      typedef itk::SmartPointer<Self> Pointer; 
+
+      itkNewMacro(Self);
+      itkTypeMacro(SARDeramp, otb::Wrapper::Application);
+
+      // Filter
+      typedef otb::SARDerampImageFilter<ComplexFloatImageType> DerampFilterType;
+
+    private:
+      void DoInit() override
+      {
+	SetName("SARDeramp");
+	SetDescription("Deramping/Reramping of S1 IW Burst.");
+
+	SetDocName("SAR Deramp");
+	SetDocLongDescription("This application does the deramping or reramping of S1 Iw burst.");
+
+	//Optional descriptors
+	SetDocLimitations("Only Sentinel 1 products are supported for now.");
+	SetDocAuthors("OTB-Team");
+	SetDocSeeAlso(" ");
+	AddDocTag(Tags::SAR);
+
+	//Parameter declarations
+	AddParameter(ParameterType_ComplexInputImage,  "in",   "Input burst");
+	SetParameterDescription("in", "Input burst (from S1 Iw product).");
+
+	AddParameter(ParameterType_ComplexOutputImage, "out", "Output burst");
+	SetParameterDescription("out","Output burst (deramped or reramped).");
+
+	AddParameter(ParameterType_Bool, "reramp", "Activate Reramping mode");
+	SetParameterDescription("reramp", "If true, then reramping else deramping.");
+
+	AddRAMParameter();
+
+	SetDocExampleParameterValue("in","s1b-s4-slc-vv-20160929t014610-20160929t014634-002277-003d71-002.tiff");
+	SetDocExampleParameterValue("out","s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002_Deramp.tiff");
+      }
+
+
+      void DoUpdateParameters() override
+      {
+	// Nothing to do here : all parameters are independent
+      }
+
+      void DoExecute() override
+      { 
+	// Check if reramp is enabled
+	bool derampMode = true;
+	if (IsParameterEnabled("reramp"))
+	  {
+	    derampMode = false;
+	  }
+
+	if (derampMode)
+	  {
+	    otbAppLogINFO(<<"Deramp Mode");
+	  }
+	else
+	  {
+	    otbAppLogINFO(<<"Reramp Mode");
+	  }
+
+	// Instanciate the Interferogram filter
+	DerampFilterType::Pointer filterDeramp = DerampFilterType::New();
+	m_Ref.push_back(filterDeramp.GetPointer());
+	filterDeramp->SetDerampMode(derampMode);
+	
+	// Pipeline
+	filterDeramp->SetInput(GetParameterComplexFloatImage("in"));
+	SetParameterComplexOutputImage("out", filterDeramp->GetOutput());
+
+      }
+      // Vector for filters 
+      std::vector<itk::ProcessObject::Pointer> m_Ref;
+
+    };
+
+  }
+
+}
+
+
+ 
+OTB_APPLICATION_EXPORT(otb::Wrapper::SARDeramp)

include/otbSARDerampImageFilter.h

@@ -31,6 +31,24 @@
 #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 SARDerampImageFilter 
@@ -59,7 +77,7 @@ public:
   itkTypeMacro(SARDerampImageFilter,ImageToImageFilter);
 
   /** Typedef to image input/output type VectorImage (Complex Image)   */
-  typedef TImage                                  ImageInType;
+  typedef TImage                                  ImageType;
   /** Typedef to describe the inout image pointer type. */
   typedef typename ImageType::Pointer             ImagePointer;
   typedef typename ImageType::ConstPointer        ImageConstPointer;
@@ -83,8 +101,53 @@ public:
   // Typedef for iterators
   typedef itk::ImageScanlineConstIterator< ImageType > InputIterator;
   typedef itk::ImageScanlineIterator< ImageType > OutputIterator;
+
+#if defined(USE_BOOST_TIME)
+   typedef boost::posix_time::ptime            TimeType;
+   typedef boost::posix_time::precise_duration DurationType;
+#else
+   typedef ossimplugins::time::ModifiedJulianDate TimeType;
+   typedef ossimplugins::time::Duration           DurationType;
+#endif
  
-  // Setter
+   // Struture declaration
+   struct FMRateRecordType
+   {
+  TimeType      azimuthFMRateTime;
+  double        coef0FMRate;
+  double        coef1FMRate;
+  double        coef2FMRate;
+  double        tau0FMRate;
+  friend std::ostream & operator<<(std::ostream & os, const FMRateRecordType & v)
+  {
+  return os << "{ azimuthFMRateTime: " << v.azimuthFMRateTime
+	    <<        ", coefficient 0: "  << v.coef0FMRate
+	    <<        ", coefficient 1: "        << v.coef1FMRate
+	    <<        ",coefficient 2: "         << v.coef2FMRate
+	    <<        ",slant range time (tau 0): "         << v.tau0FMRate
+	    <<        "}";
+}
+};
+
+   struct DCFRecordType
+   {
+  TimeType      azimuthDCFTime;
+  double        coef0DCF;
+  double        coef1DCF;
+  double        coef2DCF;
+  double        tau0DCF;
+  friend std::ostream & operator<<(std::ostream & os, const DCFRecordType & v)
+  {
+  return os << "{ azimuthDCFTime: " << v.azimuthDCFTime
+	    <<        ", coefficient 0: "  << v.coef0DCF
+	    <<        ", coefficient 1: "        << v.coef1DCF
+	    <<        ",coefficient 2: "         << v.coef2DCF
+	    <<        ",slant range time (tau 0): "         << v.tau0DCF
+	    <<        "}";
+}
+};
+  
+   // Setter
   itkSetMacro(DerampMode, bool);
 
   // Getter
@@ -122,6 +185,12 @@ protected:
   void BeforeThreadedGenerateData() ITK_OVERRIDE;
   
 
+  void getAllCoefs(ImageKeywordlist  const& kwl, std::vector<FMRateRecordType> & FMRateRecords);
+  void getAllCoefs(ImageKeywordlist  const& kwl, std::vector<DCFRecordType> & DCFRecords);
+
+double applyFMRateCoefs(unsigned int index_sample);
+double applyDCFCoefs(unsigned int index_sample);
+   
  private:
   SARDerampImageFilter(const Self&); // purposely not implemented
   void operator=(const Self &); // purposely not 
@@ -134,36 +203,40 @@ protected:
   bool m_DerampMode;
 
   // First azimuth/range time
-  float m_FirstAziTime, m_FirstRangeTime;
+  TimeType m_FirstAziTime;
+  double m_FirstRangeTime;
 
   // Mid burst zero Doppler azimuth time [s] 
-  float m_MidAziTime;
+  TimeType m_MidAziTime;
   
   // Spacecraft velocity computed at mid-burst time (m_MidAziTime) [m/s]
-  float m_VSatAtMidAziTime;
+  double m_VSatAtMidAziTime;
 
   // Doppler Centroid rate introduced by the scanning og the antenna [Hz/s]
-  float m_Ks;
+  double m_Ks;
   
   // Azimuth time interval [s]
-  float m_AziTimeInt;
+  double m_AziTimeInt;
 
   // Range sampling rate [Hz]
-  float m_RangeSamplingRate;
+  double m_RangeSamplingRate;
 
   // Parameters coefficients for Doppler FM Rate
-  float m_FM_C0, m_FM_C1, m_FM_C2, m_FM_Tau0;
+  double 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;
+  double m_DCF_C0, m_DCF_C1, m_DCF_C2, m_DCF_Tau0;
 
   // Reference time 0 [s]
-  float m_RefTime0;
+  double m_RefTime0;
+
+  // Mid Line Burst
+  double m_LineAtMidBurst;
 
   // First estimation
   bool m_FirstEstimation;
 
-  const double C;
+  const double C = 299792458;
 };
 
 } // End namespace otb

include/otbSARDerampImageFilter.txx

@@ -34,6 +34,20 @@
 #include <algorithm>
 #include <omp.h>
 
+#if defined(USE_BOOST_TIME)
+using boost::posix_time::microseconds;
+using boost::posix_time::seconds;
+#else
+using ossimplugins::time::microseconds;
+using ossimplugins::time::seconds;
+#endif
+
+#ifndef M_1_PI
+#define M_1_PI 0.31830988618379067154
+#endif
+
+
+
 namespace otb
 {
   /** 
@@ -43,9 +57,9 @@ namespace otb
   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)
+       m_DCF_C0(0), m_DCF_C1(0), m_DCF_C2(0), m_DCF_Tau0(0), m_RefTime0(0), m_FirstEstimation(true)
   {
-C = 299792458;
+
   }
     
   /** 
@@ -54,6 +68,7 @@ C = 299792458;
   template <class TImage> 
   SARDerampImageFilter< TImage >::~SARDerampImageFilter()
   {
+   
   }
 
   /**
@@ -70,6 +85,66 @@ C = 299792458;
        << std::endl;
   }
 
+  template<class TImage>
+  void
+  SARDerampImageFilter< TImage >
+  ::getAllCoefs(ImageKeywordlist  const& kwl, std::vector<FMRateRecordType> & FMRateRecords)
+  {
+    char fmRatePrefix_[1024];
+    std::size_t nbLists(0);
+    
+    nbLists = std::stoi(kwl.GetMetadataByKey("azimuthFmRate.azi_fm_rate_coef_nb_list"));
+
+    std::string FM_PREFIX = "azimuthFmRate.azi_fm_rate_coef_list";
+    
+    for (std::size_t listId=0; listId!= nbLists ; ++listId) 
+      {
+	const int pos = sprintf(fmRatePrefix_, "%s%d.", FM_PREFIX.c_str(), (listId+1));
+	assert(pos > 0 && pos < sizeof(fmRatePrefix_));
+	const std::string FMPrefix(fmRatePrefix_, pos);
+
+	FMRateRecordType fmRateRecord;
+	fmRateRecord.azimuthFMRateTime = ossimplugins::time::toModifiedJulianDate(kwl.GetMetadataByKey(FMPrefix+ 
+										      "azi_fm_rate_coef_time"));
+	fmRateRecord.coef0FMRate = std::stod(kwl.GetMetadataByKey(FMPrefix + "1.azi_fm_rate_coef"));
+	fmRateRecord.coef1FMRate = std::stod(kwl.GetMetadataByKey(FMPrefix + "2.azi_fm_rate_coef"));
+	fmRateRecord.coef2FMRate = std::stod(kwl.GetMetadataByKey(FMPrefix + "3.azi_fm_rate_coef"));
+	fmRateRecord.tau0FMRate = std::stod(kwl.GetMetadataByKey(FMPrefix + "slant_range_time"));
+
+	FMRateRecords.push_back(fmRateRecord);
+      }
+  }
+
+  template<class TImage>
+  void
+  SARDerampImageFilter< TImage >
+  ::getAllCoefs(ImageKeywordlist  const& kwl, std::vector<DCFRecordType> & DCFRecords)
+  {
+     char dcfPrefix_[1024];
+    std::size_t nbLists(0);
+    
+    nbLists = std::stoi(kwl.GetMetadataByKey("dopplerCentroid.dop_coef_nb_list"));
+
+    std::string DCF_PREFIX = "dopplerCentroid.dop_coef_list";
+    
+    for (std::size_t listId=0; listId!= nbLists ; ++listId) 
+      {
+	const int pos = sprintf(dcfPrefix_, "%s%d.", DCF_PREFIX.c_str(), (listId+1));
+	assert(pos > 0 && pos < sizeof(dcfPrefix_));
+	const std::string DCFPrefix(dcfPrefix_, pos);
+
+	DCFRecordType dcfRecord;
+	dcfRecord.azimuthDCFTime = ossimplugins::time::toModifiedJulianDate(kwl.GetMetadataByKey(DCFPrefix + 
+									      "dop_coef_time"));
+	dcfRecord.coef0DCF = std::stod(kwl.GetMetadataByKey(DCFPrefix + "1.dop_coef"));
+	dcfRecord.coef1DCF = std::stod(kwl.GetMetadataByKey(DCFPrefix + "2.dop_coef"));
+	dcfRecord.coef2DCF = std::stod(kwl.GetMetadataByKey(DCFPrefix + "3.dop_coef"));
+	dcfRecord.tau0DCF = std::stod(kwl.GetMetadataByKey(DCFPrefix + "slant_range_time"));
+
+	 DCFRecords.push_back(dcfRecord);
+      }
+  }
+
   /**
    * Method selectFMRateCoef
    */
@@ -78,7 +153,49 @@ C = 299792458;
   SARDerampImageFilter< TImage >
   ::selectFMRateCoef()
   {
+    // Get input
+    ImagePointer inputPtr = const_cast< ImageType * >( this->GetInput() );
+
+     // Retrieve all polynomials
+    std::vector<FMRateRecordType> FMRateRecords;
+    this->getAllCoefs(inputPtr->GetImageKeywordlist(), FMRateRecords);
+
+    DurationType diffAziTimeMin;
+    unsigned int polySelected = 0;
+    
+    // Select one polynomial (with m_MidAziTime)
+    for (unsigned int i = 0; i < FMRateRecords.size(); i++)
+      {
+	// Difference between midAziTime and aziTime of the current polynomial
+	DurationType diffAziTime;
+	if (m_MidAziTime > FMRateRecords[i].azimuthFMRateTime)
+	  {
+	    diffAziTime = DurationType(m_MidAziTime - FMRateRecords[i].azimuthFMRateTime);
+	  }
+	else
+	  {
+	    diffAziTime = DurationType(FMRateRecords[i].azimuthFMRateTime - m_MidAziTime);
+	  }
+	
+	if (diffAziTime < diffAziTimeMin || i == 0)
+	  {
+	    diffAziTimeMin = diffAziTime;
+	    polySelected = i;
+	  }
+      }
+
+    // Assign m_FM_C0, m_FM_C1, m_FM_C2, m_FM_Tau0
+    m_FM_C0 = FMRateRecords[polySelected].coef0FMRate;
+    m_FM_C1 = FMRateRecords[polySelected].coef1FMRate;
+    m_FM_C2 = FMRateRecords[polySelected].coef2FMRate;
+    m_FM_Tau0 = FMRateRecords[polySelected].tau0FMRate;
+
+    std::cout << "m_FM_C0 = " << m_FM_C0 << std::endl;
+    std::cout << "m_FM_C1 = " << m_FM_C1 << std::endl;
+    std::cout << "m_FM_C2 = " << m_FM_C2 << std::endl;
+    std::cout << "m_FM_Tau0 = " << m_FM_Tau0 << std::endl;
     
+    return true;
   }
 
  /**
@@ -89,7 +206,77 @@ C = 299792458;
   SARDerampImageFilter< TImage >
   ::selectDCFCoef()
   {
+     // Get input
+    ImagePointer inputPtr = const_cast< ImageType * >( this->GetInput() );
+
+     // Retrieve all polynomials
+    std::vector<DCFRecordType> DCFRecords;
+    this->getAllCoefs(inputPtr->GetImageKeywordlist(), DCFRecords);
+
+    DurationType diffAziTimeMin;
+    unsigned int polySelected = 0;
     
+    // Select one polynomial (with m_MidAziTime)
+    for (unsigned int i = 0; i < DCFRecords.size(); i++)
+      {
+	// Difference between midAziTime and aziTime of the current polynomial
+	DurationType diffAziTime;
+	if (m_MidAziTime > DCFRecords[i].azimuthDCFTime)
+	  {
+	    diffAziTime = DurationType(m_MidAziTime - DCFRecords[i].azimuthDCFTime);
+	  }
+	else
+	  {
+	    diffAziTime = DurationType(DCFRecords[i].azimuthDCFTime - m_MidAziTime);
+	  }
+	
+	if (diffAziTime < diffAziTimeMin || i == 0)
+	  {
+	    diffAziTimeMin = diffAziTime;
+	    polySelected = i;
+	  }
+      }
+
+    // Assign m_FM_C0, m_FM_C1, m_FM_C2, m_FM_Tau0
+    m_DCF_C0 = DCFRecords[polySelected].coef0DCF;
+    m_DCF_C1 = DCFRecords[polySelected].coef1DCF;
+    m_DCF_C2 = DCFRecords[polySelected].coef2DCF;
+    m_DCF_Tau0 = DCFRecords[polySelected].tau0DCF;
+
+    std::cout << "m_DCF_C0 = " << m_DCF_C0 << std::endl;
+    std::cout << "m_DCF_C1 = " << m_DCF_C1 << std::endl;
+    std::cout << "m_DCF_C2 = " << m_DCF_C2 << std::endl;
+    std::cout << "m_DCF_Tau0 = " << m_DCF_Tau0 << std::endl;
+    
+    return true;
+  }
+
+  /**
+   * Apply FM Rate coefficients Method
+   */
+  template<class TImage>
+  double
+  SARDerampImageFilter<TImage>
+  ::applyFMRateCoefs(unsigned int index_sample)
+  {
+    double slant_range_time = m_FirstRangeTime + (index_sample / m_RangeSamplingRate);
+    
+    return (m_FM_C0 + m_FM_C1*(slant_range_time - m_FM_Tau0) + m_FM_C2*(slant_range_time - m_FM_Tau0)*
+	    (slant_range_time - m_FM_Tau0));
+  }
+
+  /**
+   * Apply doppler centroid Frequency coefficients Method
+   */
+  template<class TImage>
+  double
+  SARDerampImageFilter<TImage>
+  ::applyDCFCoefs(unsigned int index_sample)
+  {
+    double slant_range_time = m_FirstRangeTime + (index_sample / m_RangeSamplingRate);
+    
+    return (m_DCF_C0 + m_DCF_C1*(slant_range_time - m_DCF_Tau0) + m_DCF_C2*(slant_range_time - m_DCF_Tau0)*
+	    (slant_range_time - m_DCF_Tau0));
   }
 
   /**
@@ -111,30 +298,30 @@ C = 299792458;
 	ImageKeywordlist inputKWL = inputPtr->GetImageKeywordlist();
 
 	// Check version of header/kwl (at least 3)
-	int headerVersion = std::stoi(inputKLW.GetMetadataByKey("header.version"));
+	int headerVersion = std::stoi(inputKWL.GetMetadataByKey("header.version"));
 
 	if (headerVersion < 3)
 	  {
-	    itkExceptionMacro(<<"Header version is < 3. Please Upgrade your geom file");
+	    itkExceptionMacro(<<"Header version is inferior to 3. Please Upgrade your geom file");
 	  }
 
 	// Get some metadata
-	float aziSteeringRate = std::stof(inputKLW.GetMetadataByKey("support_data.azimuth_steering_rate"));
+	double aziSteeringRate = std::stod(inputKWL.GetMetadataByKey("support_data.azimuth_steering_rate"));
 	// Conversion to radians per seconds
-	aziSteeringRate *= (M_PI/180); 
+	aziSteeringRate *= (M_1_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_FirstAziTime = ossimplugins::time::toModifiedJulianDate((inputKWL.GetMetadataByKey("support_data.first_line_time")));
+	m_FirstRangeTime = std::stod(inputKWL.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"));
+	m_AziTimeInt = std::stod(inputKWL.GetMetadataByKey("support_data.line_time_interval"));
+	m_RangeSamplingRate = std::stod(inputKWL.GetMetadataByKey("support_data.range_sampling_rate"));
 
-	float radarFrequency = std::stof(inputKLW.GetMetadataByKey("support_data.radar_frequency"));
+	double radarFrequency = std::stod(inputKWL.GetMetadataByKey("support_data.radar_frequency"));
 	
-	int nbLineBurst = std::stof(inputKLW.GetMetadataByKey("support_data.geom.bursts.number_lines_per_burst"));
+	int nbLineBurst = std::stod(inputKWL.GetMetadataByKey("support_data.geom.bursts.number_lines_per_burst"));
 	// Estimation m_Ks
-	double lineAtMidBurst = nbLineBurst/2.;
-	m_MidAziTime = m_FirstAziTime + m_AziTimeInt * lineAtMidBurst;
+	m_LineAtMidBurst = nbLineBurst/2.;
+	m_MidAziTime = m_FirstAziTime + seconds(m_AziTimeInt * m_LineAtMidBurst);
 	
 	// Try to create a SarSensorModelAdapter
 	SarSensorModelAdapter::Pointer sarSensorModel = SarSensorModelAdapter::New();
@@ -145,7 +332,7 @@ C = 299792458;
 
 	Point3DType satpos, satvel;
  
-	bool lineToSatPosAndVelOk = sarSensorModel->LineToSatPositionAndVelocity(lineAtMidBurst, satpos, satvel);
+	bool lineToSatPosAndVelOk = sarSensorModel->LineToSatPositionAndVelocity(m_LineAtMidBurst, satpos, satvel);
 
 	if (!lineToSatPosAndVelOk)
 	  itkExceptionMacro(<<"Failed to estimate satellite position and velocity.");
@@ -154,8 +341,17 @@ C = 299792458;
 
 	m_Ks = (2*m_VSatAtMidAziTime/C) * radarFrequency * aziSteeringRate;
 	
+	std::cout << "m_FirstAziTime = " << m_FirstAziTime << std::endl;
+	std::cout << "m_MidAziTime = " << m_MidAziTime << std::endl;
+
 	// Polynomial selection (FM Rate and Doppler Centroid Frequency)
-	
+	this->selectFMRateCoef();
+	this->selectDCFCoef();
+
+	// Estimate Reference time at first sample
+	m_RefTime0 = - (this->applyDCFCoefs(0) / this->applyFMRateCoefs(0));
+
+	std::cout << "m_RefTime0 = " << m_RefTime0 << std::endl;
       }
   }
    
@@ -166,168 +362,65 @@ C = 299792458;
   template<class TImage>
   void
   SARDerampImageFilter<TImage>
-  ::ThreadedGenerateData(const ImageOutRegionType & outputRegionForThread,
+  ::ThreadedGenerateData(const ImageRegionType & 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);
+    ImageRegionType inputRegionForThread = 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;
+    // Iterator on input 
+    InputIterator  InIt(this->GetInput(), inputRegionForThread);
+    InIt.GoToBegin();
 
     // For each line
-    while (!OutIt.IsAtEnd() && !InMasterCartMeanIt.IsAtEnd())
+    while (!OutIt.IsAtEnd() && !InIt.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;
-
+	InIt.GoToBeginOfLine();
 
-	// 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];
+	// Output index
+	ImageIndexType index = OutIt.GetIndex();
 
-	// 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;
+	// Zero Doppler azimuth Time
+	double aziTime = (index[1]-m_LineAtMidBurst)*m_AziTimeInt;
+		
+	// For each colunm
+	while (!OutIt.IsAtEndOfLine() && !InIt.IsAtEndOfLine()) 
+	  {
+	    // Reference Time
+	    double refTime =  - (this->applyDCFCoefs(index[0]) / this->applyFMRateCoefs(index[0])) - m_RefTime0;
 
-	////////// Estimate satellite positions for the current line //////////
-	bool checkSlave = m_SarSensorModelAdapterForSlave->LineToSatPositionAndVelocity(static_cast<double>(Le), satPosSlave, 
-								      satVelSlave);
+	    // Kt
+	    double Kt = (this->applyFMRateCoefs(index[0]) * m_Ks) / (this->applyFMRateCoefs(index[0]) - m_Ks);
 
-	bool checkMaster = m_SarSensorModelAdapterForMaster->LineToSatPositionAndVelocity(static_cast<double>(ind_Line), 
-								      satPosMaster, satVelMaster);
+	    // Phi
+	    double Phi = - M_1_PI * Kt * (aziTime -refTime)*(aziTime -refTime);
 
-	// For each colunm
-	while (!OutIt.IsAtEndOfLine() && !InMasterCartMeanIt.IsAtEndOfLine()) 
-	  {
-	    // Check slave Index
-	     if (checkSlave)
+	    // Deramping or Reramping 
+	    if (!m_DerampMode)
 	      {
-		// 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;
-		      }
-		  }
-
+		Phi *= -1;
 	      }
-	    else
-	      {
-		// All components set to 0
-		pixelOut[0] = 0;
-		pixelOut[1] = 0;
-		if (m_MasterCopy) 
-		  {
-		    pixelOut[2] = 0;
-		    pixelOut[3] = 0;
-		    pixelOut[4] = 0;
-		  }
-		}
 
+	    ImagePixelType outPixel; // Complex 
+	    outPixel.real((InIt.Get().real() * std::cos(Phi)) - (InIt.Get().imag() * std::sin(Phi)));
+	    outPixel.imag((InIt.Get().real() * std::sin(Phi)) + (InIt.Get().imag() * std::cos(Phi)));
+	    
 	    //////////// Assign Output ////////////
-	    OutIt.Set(pixelOut);
+	    OutIt.Set(outPixel);
 
 	    // Increment iterators
 	    ++OutIt;
-	    ++InMasterCartMeanIt;
+	    ++InIt;
 	  } // End colunms (ouput)
 
 	// Next Line
 	OutIt.NextLine();
-	InMasterCartMeanIt.NextLine();
+	InIt.NextLine();
       } // End lines (ouput)
   }