From b31695cc580db7d0ccb0d5e85df756e7565487a6 Mon Sep 17 00:00:00 2001
From: Julien Michel <julien.michel@orfeo-toolbox.org>
Date: Fri, 15 Jan 2010 16:25:48 +0100
Subject: [PATCH] ENH: In TSX calib, theta_loc is now estimated via
least-square linear regression on the incidence angle records available
---
.../otbTerraSarCalibrationFunctor.h | 27 +--
.../otbTerraSarCalibrationFunctor.txx | 12 +-
.../otbTerraSarCalibrationImageFilter.h | 59 +++++-
.../otbTerraSarCalibrationImageFilter.txx | 175 ++++++++++++++++--
.../otbTerraSarCalibrationFunctor.cxx | 4 +-
5 files changed, 222 insertions(+), 55 deletions(-)
diff --git a/Code/Radiometry/otbTerraSarCalibrationFunctor.h b/Code/Radiometry/otbTerraSarCalibrationFunctor.h
index d817c98113..6d6859c777 100644
--- a/Code/Radiometry/otbTerraSarCalibrationFunctor.h
+++ b/Code/Radiometry/otbTerraSarCalibrationFunctor.h
@@ -98,20 +98,7 @@ public:
{
return m_CalibrationFactor;
};
-
- /** Set the local incidence angle */
- void SetLocalIncidentAngle( double val )
- {
- m_LocalIncidentAngle = val;
- m_SinLocalIncidentAngle = vcl_sin(m_LocalIncidentAngle*CONST_PI_180);
- };
-
- /** Get the local incidence angle */
- double GetLocalIncidentAngle() const
- {
- return m_LocalIncidentAngle;
- };
-
+
/** Set the image size */
void SetOriginalProductSize( SizeType size )
{
@@ -155,16 +142,16 @@ public:
}
/** Get the default value */
- void GetDefaultValue(void)
+ double GetDefaultValue(void)
{
return m_DefaultValue;
}
/** Perform the calibration for one pixel (scalar -> modulus image) */
- inline TOutput operator() (const TInput & inPix, IndexType index);
+ inline TOutput operator() (const TInput & inPix, const IndexType& index, double angle);
/** Perform the calibration for one pixel (complex -> complex image) */
- inline std::complex<TOutput> operator() (const std::complex<TInput> & inPix, IndexType index);
+ inline std::complex<TOutput> operator() (const std::complex<TInput> & inPix, const IndexType& index, double angle);
private:
/** Return the current range position */
@@ -182,12 +169,6 @@ private:
/** Noise record vector */
ImageNoiseType m_NoiseRecord;
- /** Sensor local incident angle in degree */
- double m_LocalIncidentAngle;
-
- /** sin of the LocalIncidentAngle */
- double m_SinLocalIncidentAngle;
-
/** Image Size */
SizeType m_OriginalProductSize;
diff --git a/Code/Radiometry/otbTerraSarCalibrationFunctor.txx b/Code/Radiometry/otbTerraSarCalibrationFunctor.txx
index 2426d23d61..78bd9ef877 100644
--- a/Code/Radiometry/otbTerraSarCalibrationFunctor.txx
+++ b/Code/Radiometry/otbTerraSarCalibrationFunctor.txx
@@ -31,8 +31,6 @@ TerraSarCalibrationFunctor<TInput, TOutput>
{
// Initialise values
m_CalibrationFactor = itk::NumericTraits<double>::Zero;
- m_LocalIncidentAngle = itk::NumericTraits<double>::Zero;
- m_SinLocalIncidentAngle = itk::NumericTraits<double>::Zero;
m_OriginalProductSize.Fill(0);
m_UseFastCalibration = false;
m_ResultsInDecibels = false;
@@ -97,7 +95,7 @@ TerraSarCalibrationFunctor<TInput, TOutput>
template <class TInput, class TOutput>
TOutput
TerraSarCalibrationFunctor<TInput, TOutput>
-::operator()(const TInput & inPix, IndexType index)
+::operator()(const TInput & inPix, const IndexType& index, double angle)
{
// Formula: sigma = (Ks.|DN|²-NEBN) * sin Theta_local
@@ -110,7 +108,7 @@ TerraSarCalibrationFunctor<TInput, TOutput>
// If fast calibration is off, compute noise
if(m_UseFastCalibration)
{
- sigma = beta0 * m_SinLocalIncidentAngle;
+ sigma = beta0 * vcl_sin(angle);
}
else
{
@@ -121,7 +119,7 @@ TerraSarCalibrationFunctor<TInput, TOutput>
double NEBN = this->ComputeNoiseEquivalentBetaNaught(currentRange);
// Last, apply formula
- sigma = (beta0 - NEBN) * m_SinLocalIncidentAngle;
+ sigma = (beta0 - NEBN) * vcl_sin(angle);
// Handle negative sigma case
if(sigma <=0)
@@ -143,14 +141,14 @@ TerraSarCalibrationFunctor<TInput, TOutput>
template <class TInput, class TOutput>
std::complex<TOutput>
TerraSarCalibrationFunctor<TInput, TOutput>
-::operator()(const std::complex<TInput> & inPix, IndexType index)
+::operator()(const std::complex<TInput> & inPix, const IndexType& index, double angle)
{
// First, extract modulus and phase
double modulus = vcl_sqrt(inPix.real()*inPix.real() + inPix.imag()*inPix.imag());
double phase = vcl_atan2(inPix.imag(),inPix.real());
// Then, calibrate the modulus
- double sigma = this->operator()(modulus,index);
+ double sigma = this->operator()(modulus,index,angle);
// Last, put back the phase
std::complex<TOutput> out(std::polar(sigma,phase));
diff --git a/Code/Radiometry/otbTerraSarCalibrationImageFilter.h b/Code/Radiometry/otbTerraSarCalibrationImageFilter.h
index 674303f721..8860f845bc 100644
--- a/Code/Radiometry/otbTerraSarCalibrationImageFilter.h
+++ b/Code/Radiometry/otbTerraSarCalibrationImageFilter.h
@@ -45,8 +45,12 @@ namespace otb
* NEBN is interpolated for each range position according to the
* polynomial coefficients from the most recent noise record.
*
- * \f$ \theta_{loc} \f$ is the mean incident angle over the scene. For
- * now, this filter does not support more accurate behaviour like
+ * \f$ \theta_{loc} \f$ is estimated from the incidence angles from
+ * the metadata following a least-square regression.
+ * One can also add its own incident angle records using the
+ * appropriate add method.
+ *
+ * For now, this filter does not support more accurate behaviour like
* incident angle interpolation from metadata or the Geocoded
* Incidence Angle mask optional 1B product (GIM).
*
@@ -107,6 +111,11 @@ public:
typedef typename OutputImageType::InternalPixelType OutputInternalPixelType;
typedef typename itk::NumericTraits<InputInternalPixelType>::ValueType InputValueType;
typedef typename itk::NumericTraits<OutputInternalPixelType>::ValueType OutputValueType;
+ typedef typename InputImageType::IndexType IndexType;
+
+ /** Incidence angle samples */
+ typedef std::pair<IndexType,double> IncidenceAngleRecordType;
+ typedef std::vector<IncidenceAngleRecordType> IncidenceAngleRecordVectorType;
/** Calibration functor typedef */
typedef typename Functor::TerraSarCalibrationFunctor<InputValueType,
@@ -135,10 +144,6 @@ public:
/** Calibration Factor */
itkSetMacro(CalibrationFactor,double);
itkGetMacro(CalibrationFactor,double);
-
- /** Sensor local incident angle in degree */
- itkSetMacro(LocalIncidentAngle,double);
- itkGetMacro(LocalIncidentAngle,double);
/** If fast calibration is On, noise is ignored */
itkSetMacro(UseFastCalibration,bool);
@@ -179,6 +184,17 @@ public:
*/
void ClearNoiseRecords();
+ /**
+ * Add an incidence angle record to estimate the local incidence.
+ * Angle has to be given in degrees.
+ */
+ void AddIncidenceAngleRecord(const IndexType& index, double angle);
+
+ /**
+ * Clear the incidence angle records.
+ */
+ void ClearIncidenceAngleRecords();
+
protected:
/** Constructor */
TerraSarCalibrationImageFilter();
@@ -194,6 +210,10 @@ protected:
/** PrintSelf method */
void PrintSelf(std::ostream& os, itk::Indent indent) const;
+ /** Reimplement the Modified() method to handle the
+ * ParametersUpToDate flag */
+ virtual void Modified();
+
private:
TerraSarCalibrationImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
@@ -202,12 +222,19 @@ private:
* date */
static bool CompareNoiseRecords(const NoiseRecordType& record1, const NoiseRecordType& record2);
+ /** Estimate angular plane. The angular is interpolated by a plane
+ * fitted by least squares to the angular corners. More Appropriate
+ * methods (using DEM for instance) could be implemented later */
+ void EstimateAngularPlaneParameters();
+
+ /**
+ * Compute the incidence angle from the index
+ */
+ inline double ComputeIncidenceAngle(const IndexType& index) const;
+
/** Calibration Factor */
double m_CalibrationFactor;
- /** Sensor local incident angle in degree */
- double m_LocalIncidentAngle;
-
/** Pulse Repetition Frequency */
double m_PRF;
@@ -224,10 +251,22 @@ private:
NoiseRecordVectorType m_NoiseRecords;
/** Default value (for negative sigma) */
- bool m_DefaultValue;
+ double m_DefaultValue;
+
+ /** Incidence angle records vector */
+ IncidenceAngleRecordVectorType m_IncidenceAngleRecords;
+
+ /** Incidence angle linear regression parameters */
+ double m_IncidenceAngleAx;
+ double m_IncidenceAngleAy;
+ double m_IncidenceAngleOffset;
+
+ /** True if parameters are up-to-date (avoid recomputing) */
+ bool m_ParametersUpToDate;
};
+
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
diff --git a/Code/Radiometry/otbTerraSarCalibrationImageFilter.txx b/Code/Radiometry/otbTerraSarCalibrationImageFilter.txx
index 424ad42899..5e8db6198a 100644
--- a/Code/Radiometry/otbTerraSarCalibrationImageFilter.txx
+++ b/Code/Radiometry/otbTerraSarCalibrationImageFilter.txx
@@ -28,6 +28,10 @@
#include "itkImageRegionIterator.h"
#include "itkProgressReporter.h"
+#include <vnl/algo/vnl_lsqr.h>
+#include <vnl/vnl_sparse_matrix_linear_system.h>
+#include <vnl/vnl_least_squares_function.h>
+
namespace otb
{
/**
@@ -36,14 +40,17 @@ namespace otb
template <class TInputImage, class TOutputImage>
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::TerraSarCalibrationImageFilter() : m_CalibrationFactor(itk::NumericTraits<double>::Zero),
- m_LocalIncidentAngle(itk::NumericTraits<double>::Zero),
- m_PRF(1.),
- m_OriginalProductSize(),
- m_UseFastCalibration(false),
- m_ResultsInDecibels(false),
- m_NoiseRecords(),
- m_DefaultValue(0.00001)
-
+ m_PRF(1.),
+ m_OriginalProductSize(),
+ m_UseFastCalibration(false),
+ m_ResultsInDecibels(false),
+ m_NoiseRecords(),
+ m_DefaultValue(0.00001),
+ m_IncidenceAngleRecords(),
+ m_IncidenceAngleAx(0.),
+ m_IncidenceAngleAy(0.),
+ m_IncidenceAngleOffset(0.),
+ m_ParametersUpToDate(false)
{}
template <class TInputImage, class TOutputImage>
@@ -54,6 +61,16 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
this->ClearNoiseRecords();
}
+template <class TInputImage, class TOutputImage>
+void
+TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
+::Modified()
+{
+ // Call Superclass implementation
+ Superclass::Modified();
+ m_ParametersUpToDate = false;
+}
+
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
@@ -64,6 +81,9 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
// Add it to the list
m_NoiseRecords.push_back(newRecord);
+
+ // Call modified method
+ this->Modified();
}
template <class TInputImage, class TOutputImage>
@@ -72,6 +92,32 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::ClearNoiseRecords()
{
m_NoiseRecords.clear();
+
+ // Call modified method
+ this->Modified();
+}
+
+template <class TInputImage, class TOutputImage>
+void
+TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
+::AddIncidenceAngleRecord(const IndexType& index, double angle)
+{
+ IncidenceAngleRecordType record(index,angle);
+ m_IncidenceAngleRecords.push_back(record);
+
+ // Call modified method
+ this->Modified();
+}
+
+template <class TInputImage, class TOutputImage>
+void
+TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
+::ClearIncidenceAngleRecords()
+{
+ m_IncidenceAngleRecords.clear();
+
+ // Call modified method
+ this->Modified();
}
template <class TInputImage, class TOutputImage>
@@ -87,6 +133,12 @@ void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::BeforeThreadedGenerateData()
{
+ // Check if we need to update parameters
+ if(m_ParametersUpToDate)
+ {
+ return;
+ }
+
// Set the product original size
if (m_OriginalProductSize[0] == 0 && m_OriginalProductSize[1] == 0)
{
@@ -162,11 +214,29 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
}
// Incident Angle
- if (m_LocalIncidentAngle == itk::NumericTraits<double>::Zero)
+ if (m_IncidenceAngleRecords.empty())
{
if (mdIsAvailable)
{
- m_LocalIncidentAngle = lImageMetadata->GetMeanIncidenceAngles(this->GetInput()->GetMetaDataDictionary());
+ // Retrieve center incidence angle
+ double centerAngle = lImageMetadata->GetCenterIncidenceAngle(this->GetInput()->GetMetaDataDictionary());
+ IndexType centerIndex = lImageMetadata->GetCenterIncidenceAngleIndex(this->GetInput()->GetMetaDataDictionary());
+ this->AddIncidenceAngleRecord(centerIndex,centerAngle);
+
+ // Retrieve corners incidence angle
+ std::vector<double> cangles = lImageMetadata->GetCornersIncidenceAngles(this->GetInput()->GetMetaDataDictionary());
+ std::vector<IndexType> cindex = lImageMetadata->GetCornersIncidenceAnglesIndex(this->GetInput()->GetMetaDataDictionary());
+
+ std::vector<double>::const_iterator angIt = cangles.begin();
+ typename std::vector<IndexType>::const_iterator indIt = cindex.begin();
+
+ // Add each corner angle record as well
+ while(angIt != cangles.end() && indIt != cindex.end())
+ {
+ this->AddIncidenceAngleRecord(*indIt,*angIt);
+ ++angIt;
+ ++indIt;
+ }
}
else
{
@@ -176,6 +246,85 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
// Ensure that noise records are sorted by decreasing acquisition
// date
std::sort(m_NoiseRecords.begin(),m_NoiseRecords.end(), &Self::CompareNoiseRecords);
+
+ // Estimate angular plane parameters
+ this->EstimateAngularPlaneParameters();
+
+ // Mark parameters as up-to-date
+ m_ParametersUpToDate = true;
+}
+
+template <class TInputImage, class TOutputImage>
+void
+TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
+::EstimateAngularPlaneParameters()
+{
+ // Check if there are any incidence angle records
+ if(m_IncidenceAngleRecords.empty())
+ {
+ itkExceptionMacro(<<"No incidence angle records found, can not perform calibration!");
+ }
+
+ // If we do not have enough point to do the least square regression,
+ // set default value as a constant mean angle
+ if(m_IncidenceAngleRecords.size() < 3)
+ {
+ m_IncidenceAngleAx = 0.;
+ m_IncidenceAngleAy = 0.;
+ m_IncidenceAngleOffset = 0.;
+
+ typename IncidenceAngleRecordVectorType::const_iterator
+ aIt = m_IncidenceAngleRecords.begin();
+
+ // Compute the mean incidence angle
+ while(aIt != m_IncidenceAngleRecords.end())
+ {
+ m_IncidenceAngleOffset += aIt->second;
+ ++aIt;
+ }
+ // Set the offset to the mean angle
+ m_IncidenceAngleOffset/=m_IncidenceAngleRecords.size();
+ m_IncidenceAngleOffset*= M_PI / 180.;
+ return;
+ }
+
+ // Perform the plane least square estimation
+ unsigned int nbRecords = m_IncidenceAngleRecords.size();
+ vnl_sparse_matrix<double> a(nbRecords,3);
+ vnl_vector<double> b(nbRecords), bestParams(3);
+
+ // Fill the linear system
+ for(unsigned int i = 0; i < nbRecords; ++i)
+ {
+ std::cout<<"Incidence record: "<<m_IncidenceAngleRecords.at(i).first<<" <-> "<<m_IncidenceAngleRecords.at(i).second<<std::endl;
+ a(i,0) = m_IncidenceAngleRecords.at(i).first[0];
+ a(i,1) = m_IncidenceAngleRecords.at(i).first[1];
+ a(i,2) = 1.;
+ b[i] = m_IncidenceAngleRecords.at(i).second * M_PI / 180.;
+ }
+
+ // Create the linear system
+ vnl_sparse_matrix_linear_system<double> linearSystem(a,b);
+ vnl_lsqr linearSystemSolver(linearSystem);
+
+ // And solve it
+ linearSystemSolver.minimize(bestParams);
+
+ m_IncidenceAngleAx = bestParams[0];
+ m_IncidenceAngleAy = bestParams[1];
+ m_IncidenceAngleOffset = bestParams[2];
+
+ std::cout<<"Parameters :" <<bestParams<<std::endl;
+}
+
+template <class TInputImage, class TOutputImage>
+double
+TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
+::ComputeIncidenceAngle(const IndexType& index) const
+{
+ // Apply the regression
+ double angle = m_IncidenceAngleAx * index[0] + m_IncidenceAngleAy * index[1] + m_IncidenceAngleOffset;
+ return angle;
}
template <class TInputImage, class TOutputImage>
@@ -205,7 +354,6 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
// Set up the functor
CalibrationFunctorType calibrationFunctor;
calibrationFunctor.SetCalibrationFactor(m_CalibrationFactor);
- calibrationFunctor.SetLocalIncidentAngle(m_LocalIncidentAngle);
calibrationFunctor.SetUseFastCalibration(m_UseFastCalibration);
calibrationFunctor.SetResultsInDecibels(m_ResultsInDecibels);
calibrationFunctor.SetOriginalProductSize(m_OriginalProductSize);
@@ -272,7 +420,7 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
}
}
// Apply the calibration functor
- outputIt.Set( calibrationFunctor( inputIt.Get(), inputIt.GetIndex() ) );
+ outputIt.Set( calibrationFunctor( inputIt.Get(), inputIt.GetIndex(), this->ComputeIncidenceAngle(inputIt.GetIndex()) ) );
++inputIt;
++outputIt;
progress.CompletedPixel(); // potential exception thrown here
@@ -290,10 +438,11 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
os << indent << "Calibration Factor: " << m_CalibrationFactor << std::endl;
os << indent << "PRF: "<<m_PRF <<std::endl;
os << indent << "Original product size: "<<m_OriginalProductSize << std::endl;
- os << indent << "Sensor local incidence angle: " << m_LocalIncidentAngle << std::endl;
os << indent << "Fast calibration: " << (m_UseFastCalibration ? "On" : "Off")<<std::endl;
os << indent << "Results in decibels: " << (m_ResultsInDecibels ? "Yes" : "No") << std::endl;
os << indent << "Number of noise records: " << m_NoiseRecords.size() <<std::endl;
+ os << indent << "Number of angle records: " << m_IncidenceAngleRecords.size() << std::endl;
+ os << indent << "Angle regression: "<<m_IncidenceAngleAx<<"* col + "<<m_IncidenceAngleAx<<" * row +"<<m_IncidenceAngleOffset<<std::endl;
}
diff --git a/Testing/Code/Radiometry/otbTerraSarCalibrationFunctor.cxx b/Testing/Code/Radiometry/otbTerraSarCalibrationFunctor.cxx
index a81652d684..a0fadf9d55 100644
--- a/Testing/Code/Radiometry/otbTerraSarCalibrationFunctor.cxx
+++ b/Testing/Code/Radiometry/otbTerraSarCalibrationFunctor.cxx
@@ -34,10 +34,10 @@ int otbTerraSarCalibrationFunctor(int argc, char * argv[])
id[0] = 125; id[0] = 150;
ScalarType inPix = 150.2;
- std::cout << inPix << " -> " << funct.operator()(inPix, id) << std::endl;
+ std::cout << inPix << " -> " << funct.operator()(inPix, id, 0.) << std::endl;
ComplexType inCplxPix(12, 180);
- std::cout << inCplxPix << " -> " << funct.operator()(inCplxPix, id) << std::endl;
+ std::cout << inCplxPix << " -> " << funct.operator()(inCplxPix, id,0.) << std::endl;
return EXIT_SUCCESS;
}
--
GitLab