From 7c62d1a24c18a49a0ba08d49858ef8e228cb41d8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ga=C3=ABlle=20USSEGLIO?= <gaelle.usseglio@cnes.fr>
Date: Mon, 3 Dec 2018 16:32:00 +0000
Subject: [PATCH] ENH : New application SARDeramp (part 2)
---
app/CMakeLists.txt | 5 +
app/otbSARDeramp.cxx | 127 +++++++++
include/otbSARDerampImageFilter.h | 97 ++++++-
include/otbSARDerampImageFilter.txx | 393 +++++++++++++++++-----------
4 files changed, 460 insertions(+), 162 deletions(-)
create mode 100644 app/otbSARDeramp.cxx
diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 25a31f4..101e0fd 100644
--- a/app/CMakeLists.txt
+++ b/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}
+)
diff --git a/app/otbSARDeramp.cxx b/app/otbSARDeramp.cxx
new file mode 100644
index 0000000..3ff808c
--- /dev/null
+++ b/app/otbSARDeramp.cxx
@@ -0,0 +1,127 @@
+/*
+ * 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)
diff --git a/include/otbSARDerampImageFilter.h b/include/otbSARDerampImageFilter.h
index b10fa48..ff583a7 100644
--- a/include/otbSARDerampImageFilter.h
+++ b/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
diff --git a/include/otbSARDerampImageFilter.txx b/include/otbSARDerampImageFilter.txx
index 709089c..fcf4204 100644
--- a/include/otbSARDerampImageFilter.txx
+++ b/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)
}
--
GitLab