Skip to content
Snippets Groups Projects
Commit bf2ada1a authored by Julien Michel's avatar Julien Michel
Browse files

BUG: Fixing TerraSarX calibration code

parent 5a6a23c5
Branches
Tags
No related merge requests found
Showing
with 405 additions and 484 deletions
......@@ -79,13 +79,13 @@ public:
typedef std::vector<DoubleVectorType> DoubleVectorVectorType;
/** Accessors */
void SetCalFactor(double pCalFactor)
void SetCalibrationFactor(double pCalibrationFactor)
{
this->GetFunctor().SetCalFactor( pCalFactor );
this->GetFunctor().SetCalibrationFactor( pCalibrationFactor );
this->Modified();
}
double GetCalFactor(){
return this->GetFunctor().GetCalFactor();
double GetCalibrationFactor(){
return this->GetFunctor().GetCalibrationFactor();
}
protected:
......
......@@ -19,8 +19,6 @@
#define __otbTerraSarBrightnessImageFilter_txx
#include "otbTerraSarBrightnessImageFilter.h"
//#include "otbImageMetadataInterfaceFactory.h"
//#include "otbImageMetadataInterfaceBase.h"
#include "otbTerraSarImageMetadataInterface.h"
namespace otb
......@@ -38,11 +36,11 @@ TerraSarBrightnessImageFilter<TInputImage,TOutputImage>
bool mdIsAvailable = lImageMetadata->CanRead(this->GetInput()->GetMetaDataDictionary());
// If the user doesn't set it AND the metadata is available, set calFactor using image metadata
if (this->GetCalFactor() == itk::NumericTraits<double>::min())
if (this->GetCalibrationFactor() == itk::NumericTraits<double>::min())
{
if (mdIsAvailable)
{
this->SetCalFactor( lImageMetadata->GetCalibrationFactor(this->GetInput()->GetMetaDataDictionary()) );
this->SetCalibrationFactor( lImageMetadata->GetCalibrationFactor(this->GetInput()->GetMetaDataDictionary()) );
}
else
{
......
......@@ -26,9 +26,6 @@
#include "otbUnaryFunctorWithIndexImageFilter.h"
#include "otbTerraSarFunctors.h"
#include "itkMetaDataDictionary.h"
//#include "itkUnaryFunctorImageFilter.h"
//#include "itkConstNeighborhoodIterator.h"
#include "otbMath.h"
namespace otb
{
......@@ -65,14 +62,18 @@ public:
typedef typename OutputImageType::InternalPixelType OutputInternalPixelType;
typedef typename itk::NumericTraits<InputInternalPixelType>::ValueType InputValueType;
typedef typename itk::NumericTraits<OutputInternalPixelType>::ValueType OutputValueType;
typedef typename Functor::TerraSarCalibrationImageFunctor< InputValueType, OutputValueType> FunctorType;
typedef typename Functor::TerraSarCalibrationImageFunctor< InputValueType, OutputValueType> FunctorType;
/** typedef to access metadata */
typedef itk::MetaDataDictionary MetaDataDictionaryType;
typedef typename InputImageType::SizeType SizeType;
/** "typedef" for standard classes. */
typedef TerraSarCalibrationImageFilter Self;
typedef UnaryFunctorWithIndexImageFilter< InputImageType, OutputImageType, FunctorType > Superclass;
typedef TerraSarCalibrationImageFilter Self;
typedef UnaryFunctorWithIndexImageFilter< InputImageType, OutputImageType, FunctorType > Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** object factory method. */
itkNewMacro(Self);
......@@ -80,60 +81,39 @@ public:
// Use a with neighborhood to have access to the pixel coordinates
itkTypeMacro(TerraSarCalibrationImageFilter, UnaryFunctorWithIndexImageFilter);
typedef itk::MetaDataDictionary MetaDataDictionaryType;
typedef typename FunctorType::DoubleVectorType DoubleVectorType;
typedef typename FunctorType::DoubleVectorVectorType DoubleVectorVectorType;
typedef typename InputImageType::SizeType SizeType;
// typedef typename FunctorType::LIntVectorType LIntVectorType;
/** Accessors */
/** Calibration Factor */
void SetCalFactor( double val );
double GetCalFactor() const;
/** Noise minimal range validity */
void SetNoiseRangeValidityMin( double val );
double GetNoiseRangeValidityMin() const;
/** Noise maximal range validity */
void SetNoiseRangeValidityMax( double val );
double GetNoiseRangeValidityMax() const;
/** Noise reference range validity */
void SetNoiseRangeValidityRef( double val );
double GetNoiseRangeValidityRef() const;
void SetCalibrationFactor( double val );
double GetCalibrationFactor() const;
/** Sensor local incident angle in degree */
void SetLocalIncidentAngle( double val );
double GetLocalIncidentAngle() const;
/** Sinus of the sensor local incident angle in degree */
double GetSinLocalIncidentAngle() const;
/** Vector of vector that contain noise polinomial coefficient */
void SetNoisePolynomialCoefficientsList( DoubleVectorVectorType vect );
DoubleVectorVectorType GetNoisePolynomialCoefficientsList() const;
/** Image size (setter is protected)*/
SizeType GetImageSize() const;
/** Fast Calibration Method. If set to trus, will consider only the first noise coefficient else,
* will use all of them and applied it according to its acquisition UTC time and the coordinates
* of the pixel in the image. */
void SetUseFastCalibrationMethod( bool b );
bool GetUseFastCalibrationMethod() const;
/** TimeUTC for each noise coefficient acquisition (in Julian day). */
void SetTimeUTC( DoubleVectorType vect );
DoubleVectorType GetTimeUTC() const;
/** Pulse Repetition Frequency */
void SetPRF( double val );
double GetPRF() const;
/** Inverse Pulse Repetition Frequency */
double GetInvPRF() const;
/** Image Size setter. This is provided in case of calbrating an
* extract of a larger scene. It does not interfer with output image
* regions */
void SetImageSize( SizeType size );
/** Get the image size */
const SizeType& GetImageSize() const;
protected:
/** Constructor */
TerraSarCalibrationImageFilter();
/** Destructor */
virtual ~TerraSarCalibrationImageFilter() {};
/** Image Size setter*/
void SetImageSize( SizeType size );
/** Initialize the functor vector */
void BeforeThreadedGenerateData();
......@@ -141,11 +121,9 @@ protected:
void PrintSelf(std::ostream& os, itk::Indent indent) const;
private:
bool m_ParamLoaded;
// TODO: add private copy and reference construtor
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
......
......@@ -20,8 +20,6 @@
#include <algorithm>
#include "otbTerraSarCalibrationImageFilter.h"
//#include "otbImageMetadataInterfaceFactory.h"
//#include "otbImageMetadataInterfaceBase.h"
#include "otbTerraSarImageMetadataInterface.h"
namespace otb
......@@ -32,119 +30,75 @@ namespace otb
template <class TInputImage, class TOutputImage>
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::TerraSarCalibrationImageFilter()
{
m_ParamLoaded = false;
}
{}
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::BeforeThreadedGenerateData()
{
Superclass::BeforeThreadedGenerateData();
if(m_ParamLoaded)
return;
this->SetImageSize( this->GetInput()->GetLargestPossibleRegion().GetSize() );
// Set the image size
if (this->GetImageSize()[0] == 0 && this->GetImageSize()[1] == 0)
{
this->SetImageSize( this->GetInput()->GetLargestPossibleRegion().GetSize() );
}
// Load metadata
// Build the metadata interface
TerraSarImageMetadataInterface::Pointer lImageMetadata = otb::TerraSarImageMetadataInterface::New();
// Check availability
bool mdIsAvailable = lImageMetadata->CanRead(this->GetInput()->GetMetaDataDictionary());
// If the user doesn't set the data AND the metadata is available, set the data
// If the user did not set the data AND the metadata is available, set the data
// CalFactor
if (this->GetCalFactor() == itk::NumericTraits<double>::min())
// CalibrationFactor
if (this->GetCalibrationFactor() == itk::NumericTraits<double>::min())
{
if (mdIsAvailable)
{
this->SetCalFactor( lImageMetadata->GetCalibrationFactor(this->GetInput()->GetMetaDataDictionary()) );
this->SetCalibrationFactor( lImageMetadata->GetCalibrationFactor(this->GetInput()->GetMetaDataDictionary()) );
}
else
{
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supproted");
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supported");
}
}
// NoiseRangeValidityMin
if (this->GetNoiseRangeValidityMin() == itk::NumericTraits<double>::min())
{
if (mdIsAvailable)
{
DoubleVectorType vecTmp = lImageMetadata->GetNoiseValidityRangeMinList(this->GetInput()->GetMetaDataDictionary());
std::sort(vecTmp.begin(), vecTmp.end());
this->SetNoiseRangeValidityMin(vecTmp[0]);
}
else
{
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supproted");
}
}
// NoiseRangeValidityMax
if (this->GetNoiseRangeValidityMax() == itk::NumericTraits<double>::max())
{
if (mdIsAvailable)
{
DoubleVectorType vecTmp = lImageMetadata->GetNoiseValidityRangeMaxList(this->GetInput()->GetMetaDataDictionary());
std::sort(vecTmp.begin(), vecTmp.end());
this->SetNoiseRangeValidityMax(vecTmp[vecTmp.size()-1]);
}
else
{
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supproted");
}
}
// NoiseRangeValidityRef
if (this->GetNoiseRangeValidityRef() == itk::NumericTraits<double>::Zero)
{
if (mdIsAvailable)
// Noise records
if(this->GetFunctor().GetNumberOfNoiseRecords() == 0)
{
double sum = 0;
// Get vector
DoubleVectorType noiseRefList = lImageMetadata->GetNoiseReferencePointList(this->GetInput()->GetMetaDataDictionary());
// Mean computation
for (int i=0; i<noiseRefList.size(); i++)
if (mdIsAvailable)
{
sum += noiseRefList[i];
}
// Retrieve the number of noise records
unsigned int nbNoiseRecords = lImageMetadata->GetNumberOfNoiseRecords(this->GetInput()->GetMetaDataDictionary());
this->SetNoiseRangeValidityRef( sum / noiseRefList.size() );
}
else
{
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supproted");
}
}
// NoisePolynomialCoefficient
if ( ((this->GetNoisePolynomialCoefficientsList())[0][0] == itk::NumericTraits<double>::Zero)
&& (this->GetNoisePolynomialCoefficientsList().size() == 1)
&& (this->GetNoisePolynomialCoefficientsList()[0].size() == 1) )
{
if (mdIsAvailable)
{
this->SetNoisePolynomialCoefficientsList(lImageMetadata->GetNoisePolynomialCoefficientsList(this->GetInput()->GetMetaDataDictionary()));
}
else
{
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supproted");
}
}
// Time UTC
if ((this->GetTimeUTC()[0] == itk::NumericTraits<double>::Zero) && (this->GetTimeUTC()[1] == 1.) && (this->GetTimeUTC().size() == 2))
{
if (mdIsAvailable)
{
this->SetTimeUTC(lImageMetadata->GetNoiseTimeUTCList(this->GetInput()->GetMetaDataDictionary()));
}
// Retrieve the noise records
ossimplugins::Noise * noiseRecords = lImageMetadata->GetNoise(this->GetInput()->GetMetaDataDictionary());
// Retrieve corresponding times
std::vector<double> noiseTimes = lImageMetadata->GetNoiseTimeUTCList(this->GetInput()->GetMetaDataDictionary());
if(noiseTimes.empty() || nbNoiseRecords)
{
itkExceptionMacro(<<"No noise records found.");
}
// Set the acquisition start time
this->GetFunctor().SetAcquisitionStartTime(noiseTimes[0]);
for(unsigned int i = 0; i < nbNoiseRecords;++i)
{
this->GetFunctor().AddNoiseRecord(noiseTimes[i],noiseRecords->get_imageNoise()[i]);
}
// Free memory
delete noiseRecords;
}
else
{
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supproted");
{
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supported");
}
}
}
// Radar frequency (PRF)
if (this->GetPRF() == 1.)
......@@ -155,7 +109,7 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
}
else
{
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supproted");
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supported");
}
}
......@@ -169,84 +123,29 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
}
else
{
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supproted");
itkExceptionMacro(<<"Invalid input image. Only TerraSar images are supported");
}
}
}
/******************************* ACCESSORS ********************************/
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::SetCalFactor( double val )
{
this->GetFunctor().SetCalFactor( val );
this->Modified();
}
template <class TInputImage, class TOutputImage>
double
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::GetCalFactor() const
{
return this->GetFunctor().GetCalFactor();
}
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::SetNoiseRangeValidityMin( double val )
{
this->GetFunctor().SetNoiseRangeValidityMin( val );
this->Modified();
}
template <class TInputImage, class TOutputImage>
double
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::GetNoiseRangeValidityMin() const
{
return this->GetFunctor().GetNoiseRangeValidityMin();
}
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::SetNoiseRangeValidityMax( double val )
::SetCalibrationFactor( double val )
{
this->GetFunctor().SetNoiseRangeValidityMax( val );
this->GetFunctor().SetCalibrationFactor( val );
this->Modified();
}
template <class TInputImage, class TOutputImage>
double
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::GetNoiseRangeValidityMax() const
{
return this->GetFunctor().GetNoiseRangeValidityMax();
}
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::SetNoiseRangeValidityRef( double val )
{
this->GetFunctor().SetNoiseRangeValidityRef( val );
this->Modified();
}
template <class TInputImage, class TOutputImage>
double
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::GetNoiseRangeValidityRef() const
::GetCalibrationFactor() const
{
return this->GetFunctor().GetNoiseRangeValidityRef();
return this->GetFunctor().GetCalibrationFactor();
}
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
......@@ -264,32 +163,6 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
return this->GetFunctor().GetLocalIncidentAngle();
}
template <class TInputImage, class TOutputImage>
double
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::GetSinLocalIncidentAngle() const
{
return this->GetFunctor().GetSinLocalIncidentAngle();
}
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::SetNoisePolynomialCoefficientsList( DoubleVectorVectorType vect )
{
this->GetFunctor().SetNoisePolynomialCoefficientsList( vect );
this->Modified();
}
template <class TInputImage, class TOutputImage>
typename TerraSarCalibrationImageFilter<TInputImage,TOutputImage>::DoubleVectorVectorType
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::GetNoisePolynomialCoefficientsList() const
{
return this->GetFunctor().GetNoisePolynomialCoefficientsList();
}
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
......@@ -301,7 +174,7 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
template <class TInputImage, class TOutputImage>
typename TerraSarCalibrationImageFilter<TInputImage,TOutputImage>::SizeType
const typename TerraSarCalibrationImageFilter<TInputImage,TOutputImage>::SizeType&
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::GetImageSize() const
{
......@@ -325,26 +198,6 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
return this->GetFunctor().GetUseFastCalibrationMethod();
}
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::SetTimeUTC( DoubleVectorType vect )
{
this->GetFunctor().SetTimeUTC( vect );
this->Modified();
}
template <class TInputImage, class TOutputImage>
typename TerraSarCalibrationImageFilter<TInputImage,TOutputImage>::DoubleVectorType
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::GetTimeUTC() const
{
return this->GetFunctor().GetTimeUTC();
}
template <class TInputImage, class TOutputImage>
void
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
......@@ -363,15 +216,6 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
return this->GetFunctor().GetPRF();
}
template <class TInputImage, class TOutputImage>
double
TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
::GetInvPRF() const
{
return this->GetFunctor().GetInvPRF();
}
/**PrintSelf method */
template <class TInputImage, class TOutputImage>
void
......@@ -380,37 +224,8 @@ TerraSarCalibrationImageFilter<TInputImage,TOutputImage>
{
Superclass::PrintSelf(os, indent);
os << "Calibration Factor : " << this->GetCalFactor() << std::endl;
os << "Calibration Factor : " << this->GetCalibrationFactor() << std::endl;
os << "Sensor local incident angle : " << this->GetLocalIncidentAngle() << std::endl;
os << "Noise minimal range validity: " << this->GetNoiseRangeValidityMin() << std::endl;
os << "Noise maximal range validity: " << this->GetNoiseRangeValidityMax() << std::endl;
os << "Noise reference range : " << this->GetNoiseRangeValidityRef() << std::endl;
if(this->GetUseFastCalibrationMethod())
{
os << "Noise polinomial coefficient: [ ";
for (unsigned int i=0; i<this->GetNoisePolynomialCoefficientsList()[0].size(); ++i)
{
os << this->GetNoisePolynomialCoefficientsList()[0][i] << " ";
}
os << "]" << std::endl;
}
else
{
os << "Pulse Repetition Frequency : " << this->GetPRF() << std::endl;
os << "Noise acquisitions : " << this->GetNoisePolynomialCoefficientsList().size() << std::endl;
for (unsigned int i=0; i<this->GetNoisePolynomialCoefficientsList().size(); ++i)
{
os << "Noise acquisition "<< i << ":" << std::endl;
os << "Noise TimeUTC : " << this->GetTimeUTC()[i] << std::endl;
os << "Noise polinomial coefficient: [ ";
for (unsigned int j=0; j<this->GetNoisePolynomialCoefficientsList()[j].size(); ++j)
{
os << this->GetNoisePolynomialCoefficientsList()[i][j] << " ";
}
os << "]" << std::endl;
}
}
}
......
......@@ -22,17 +22,17 @@
#ifndef __otbTerraSarFunctors_h
#define __otbTerraSarFunctors_h
#include "itkUnaryFunctorImageFilter.h"
#include "itkMetaDataDictionary.h"
#include "otbMath.h"
#include "itkSize.h"
#include "itkIndex.h"
#include <complex>
#include <otb/ImageNoise.h>
namespace otb
{
namespace Functor
{
/**
* \class TerraSarBrightnessImageFunctor
* \brief Compute the radar brightness from an modulus image.
......@@ -52,8 +52,8 @@ public:
typedef itk::Size<2> SizeType;
/** Accessors */
void SetCalFactor( double val ) { m_CalFactor = val; };
double GetCalFactor() { return m_CalFactor; };
void SetCalibrationFactor( double val ) { m_CalibrationFactor = val; };
double GetCalibrationFactor() { return m_CalibrationFactor; };
/** We assume that the input pixel is a scalar -> modulus image */
inline TOutput operator() (const TInput & inPix);
......@@ -62,7 +62,7 @@ public:
private:
/** Calibration Factor */
double m_CalFactor;
double m_CalibrationFactor;
};
......@@ -79,85 +79,176 @@ template<class TInput, class TOutput>
class TerraSarCalibrationImageFunctor
{
public:
/** Constructor */
TerraSarCalibrationImageFunctor();
/** Destructor */
virtual ~TerraSarCalibrationImageFunctor() {};
typedef std::vector<double> DoubleVectorType;
typedef std::vector<DoubleVectorType> DoubleVectorVectorType;
typedef std::vector<long int> LIntVectorType;
typedef itk::Size<2> SizeType;
typedef itk::Index<2> IndexType;
/** Typedef to define the noise records map */
typedef ossimplugins::ImageNoise ImageNoiseType;
typedef std::pair<double,ImageNoiseType> NoiseRecordType;
typedef std::vector<NoiseRecordType > NoiseRecordVectorType;
/** Typedef for image size and index */
typedef itk::Size<2> SizeType;
typedef itk::Index<2> IndexType;
/** Typedef for the brightness functor */
typedef TerraSarBrightnessImageFunctor<double, double> BrightnessFunctorType;
/** Accessors */
void SetCalFactor( double val ) { m_CalFactor = val; m_Brightness.SetCalFactor(val); };
double GetCalFactor() const { return m_CalFactor; };
void SetNoiseRangeValidityMin( double val ) { m_NoiseRangeValidityMin = val; };
double GetNoiseRangeValidityMin() const { return m_NoiseRangeValidityMin; };
void SetNoiseRangeValidityMax( double val ) { m_NoiseRangeValidityMax = val; };
double GetNoiseRangeValidityMax() const { return m_NoiseRangeValidityMax; };
void SetNoiseRangeValidityRef( double val ) { m_NoiseRangeValidityRef = val; };
double GetNoiseRangeValidityRef() const { return m_NoiseRangeValidityRef; };
/**
* Add a new noise record to the map
*/
void AddNoiseRecord(double utcAcquisitionTime, const ImageNoiseType& record)
{
// Create the pair
NoiseRecordType newRecord(utcAcquisitionTime,record);
// Add it to the list
m_NoiseRecords.push_back(newRecord);
}
/** returns the number of noise records */
unsigned int GetNumberOfNoiseRecords() const
{
return m_NoiseRecords.size();
}
/**
* Clear all noise records
*/
void ClearNoiseRecords()
{
m_NoiseRecords.clear();
}
/** Set the calibration factor */
void SetCalibrationFactor( double val )
{
m_CalibrationFactor = val;
m_Brightness.SetCalibrationFactor(val);
};
/** Get the calibration factor */
double GetCalibrationFactor() const
{
return m_CalibrationFactor;
};
/** Set the local incidence angle */
void SetLocalIncidentAngle( double val )
{
m_LocalIncidentAngle = val;
m_SinLocalIncidentAngle = vcl_sin(m_LocalIncidentAngle*CONST_PI_180);
};
double GetLocalIncidentAngle() const { return m_LocalIncidentAngle; };
double GetSinLocalIncidentAngle() const { return m_SinLocalIncidentAngle; };
void SetNoisePolynomialCoefficientsList( DoubleVectorVectorType vect ) { m_NoisePolynomialCoefficientsList = vect; };
DoubleVectorVectorType GetNoisePolynomialCoefficientsList() const { return m_NoisePolynomialCoefficientsList; };
void SetImageSize( SizeType size ) { m_ImageSize = size; };
SizeType GetImageSize() const { return m_ImageSize; };
void SetUseFastCalibrationMethod( bool b ) { m_UseFastCalibrationMethod = b; };
bool GetUseFastCalibrationMethod() const { return m_UseFastCalibrationMethod; };
void SetTimeUTC( DoubleVectorType vect ) { m_TimeUTC = vect; };
DoubleVectorType GetTimeUTC() const { return m_TimeUTC; };
void SetPRF( double val ) { m_PRF = val; m_InvPRF = 1./m_PRF; };
double GetPRF() const { return m_PRF; };
double GetInvPRF() const { return m_InvPRF; };
BrightnessFunctorType GetBrightness() { return m_Brightness; };
double ComputeCurrentNoise( unsigned int colId );
DoubleVectorType ComputeCurrentCoeffs( unsigned int lineId );
/** We assume that the input pixel is a scalar -> modulus image */
/** Get the local incidence angle */
double GetLocalIncidentAngle() const
{
return m_LocalIncidentAngle;
};
/** Set the image size */
void SetImageSize( SizeType size )
{
m_ImageSize = size;
};
/** Get the image size */
const SizeType & GetImageSize() const
{
return m_ImageSize;
};
/** Set the UseFastCalibrationMethod flag */
void SetUseFastCalibrationMethod( bool b )
{
m_UseFastCalibrationMethod = b;
};
/** Get the UseFastCalibrationMethod flag */
bool GetUseFastCalibrationMethod() const
{
return m_UseFastCalibrationMethod;
};
/** Set The Pulse Repetition Frequency (PRF) */
void SetPRF( double val )
{
m_PRF = val;
m_InvPRF = 1./m_PRF;
};
/** Get the PRF */
double GetPRF() const
{
return m_PRF;
};
/** Set the acquisition start time */
void SetAcquisitionStartTime(double startTime)
{
m_AcquisitionStartTime = startTime;
}
/** Get the acqusition start time */
double GetAcquisitionStartTime(void) const
{
return m_AcquisitionStartTime;
}
/** Perform the calibration for one pixel (scalar -> modulus image) */
inline TOutput operator() (const TInput & inPix, IndexType index);
/** We assume that the input pixel is a complex -> complex image */
/** Perform the calibration for one pixel (complex -> complex image) */
inline std::complex<TOutput> operator() (const std::complex<TInput> & inPix, IndexType index);
private:
/** Return the closest noise record in the noise record vector */
inline const ImageNoiseType & FindClosestNoiseRecord(double utcTime) const;
/** Return the current azimuth position (current acquisition time) */
inline double ComputeAzimuthPosition(const IndexType& index) const;
/** Return the current range position */
inline double ComputeRangePosition(const ImageNoiseType& noise, const IndexType& index) const;
/** Return the current NEBN value */
inline double ComputeNoiseEquivalentBetaNaught(const ImageNoiseType& noise, double range) const;
/** Calibration Factor */
double m_CalFactor;
/** Noise minimal range validity */
double m_NoiseRangeValidityMin;
/** Noise maxinimal range validity */
double m_NoiseRangeValidityMax;
/** Noise reference range */
double m_NoiseRangeValidityRef;
double m_CalibrationFactor;
/** Radar Brightness functor */
BrightnessFunctorType m_Brightness;
/** Noise record vector */
NoiseRecordVectorType m_NoiseRecords;
/** Sensor local incident angle in degree */
double m_LocalIncidentAngle;
/** sin of the LocalIncidentAngle */
double m_SinLocalIncidentAngle;
/** Vector of vector that contain noise polinomial coefficient */
DoubleVectorVectorType m_NoisePolynomialCoefficientsList;
/** Image Size */
SizeType m_ImageSize;
/** Fast Calibration Method. If set to trus, will consider only the first noise coefficient else,
/** Fast Calibration Method. If set to true, will consider only the first noise coefficient else,
* will use all of them and applied it according to its acquisition UTC time and the coordinates
* of the pixel in the image. */
bool m_UseFastCalibrationMethod;
/** TimeUTC for each noise coefficient acquisition (in Julian day). */
DoubleVectorType m_TimeUTC;
/** Acquisition start time */
double m_AcquisitionStartTime;
/** Pulse Repetition Frequency */
double m_PRF;
/** Inverse Pulse Repetition Frequency */
double m_InvPRF;
/** Radar Brightness functor */
BrightnessFunctorType m_Brightness;
};
}// end namespace functor
} // end namespace otb
......
......@@ -33,7 +33,7 @@ template <class TInput, class TOutput>
TerraSarBrightnessImageFunctor<TInput, TOutput>
::TerraSarBrightnessImageFunctor()
{
m_CalFactor = itk::NumericTraits<double>::min();
m_CalibrationFactor = itk::NumericTraits<double>::min();
}
......@@ -42,9 +42,13 @@ TOutput
TerraSarBrightnessImageFunctor<TInput, TOutput>
::operator() (const TInput & inPix)
{
// Formula: Beta^0 = Ks * |DN|²
// First, square the input pixel
double squareInPix = vcl_pow( static_cast<double>(inPix), 2.);
// Beta naught computation
double beta = m_CalFactor*squareInPix;
// Then apply the calibration factor
double beta = m_CalibrationFactor*squareInPix;
return static_cast<TOutput>(beta);
}
......@@ -75,80 +79,108 @@ template <class TInput, class TOutput>
TerraSarCalibrationImageFunctor<TInput, TOutput>
::TerraSarCalibrationImageFunctor()
{
m_CalFactor = itk::NumericTraits<double>::min();
m_NoiseRangeValidityMin = itk::NumericTraits<double>::min();
m_NoiseRangeValidityMax = itk::NumericTraits<double>::max();
m_NoiseRangeValidityRef = itk::NumericTraits<double>::Zero;
m_LocalIncidentAngle = itk::NumericTraits<double>::Zero;
m_SinLocalIncidentAngle = itk::NumericTraits<double>::Zero;
m_NoisePolynomialCoefficientsList = DoubleVectorVectorType( 1, DoubleVectorType(1, itk::NumericTraits<double>::Zero) );
m_ImageSize.Fill(0);
m_UseFastCalibrationMethod = true;
m_TimeUTC = DoubleVectorType(2, itk::NumericTraits<double>::Zero);
m_TimeUTC[1] = 1.;
m_PRF = 1.;
// Initialise values
m_CalibrationFactor = itk::NumericTraits<double>::min();
m_LocalIncidentAngle = itk::NumericTraits<double>::Zero;
m_SinLocalIncidentAngle = itk::NumericTraits<double>::Zero;
m_ImageSize.Fill(0);
m_UseFastCalibrationMethod = false;
m_PRF = 1.;
m_InvPRF = 1.;
m_AcquisitionStartTime = itk::NumericTraits<double>::Zero;
}
template <class TInput, class TOutput>
double
const typename TerraSarCalibrationImageFunctor<TInput, TOutput>
::ImageNoiseType &
TerraSarCalibrationImageFunctor<TInput, TOutput>
::ComputeCurrentNoise( unsigned int colId )
::FindClosestNoiseRecord( double utcTime ) const
{
double curRange = 0.;
double width_2 = static_cast<double>(m_ImageSize[0])/2.;
// Use +1 because image start index is 0
if( colId < static_cast<unsigned int>(width_2) )
// Iterators
NoiseRecordVectorType::const_iterator it, closest;
// Initial search values
it = m_NoiseRecords.begin();
double delta = vcl_abs(utcTime - it->first);
closest = it;
// First increment
++it;
// Iterate on the noise records to find the closest
while(it != m_NoiseRecords.end())
{
curRange = m_NoiseRangeValidityMin + ( m_NoiseRangeValidityRef-m_NoiseRangeValidityMin )/width_2 * static_cast<double>(colId+1);
double newDelta = vcl_abs(utcTime-it->first);
// Check if this one is closer
if(newDelta < delta)
{
// If so, update search values
delta = newDelta;
closest = it;
}
// Increment
++it;
}
else
{
curRange = m_NoiseRangeValidityRef + ( m_NoiseRangeValidityMax-m_NoiseRangeValidityRef )/width_2 * (static_cast<double>(colId+1) - width_2 );
}
return curRange;
return closest->second;
}
template <class TInput, class TOutput>
double
TerraSarCalibrationImageFunctor<TInput, TOutput>
::ComputeAzimuthPosition(const IndexType & index) const
{
double currentTime = m_AcquisitionStartTime + (m_ImageSize[1]-index[1]-1) * m_InvPRF;
return currentTime;
}
template <class TInput, class TOutput>
typename TerraSarCalibrationImageFunctor<TInput, TOutput>::DoubleVectorType
double
TerraSarCalibrationImageFunctor<TInput, TOutput>
::ComputeCurrentCoeffs( unsigned int lineId )
::ComputeRangePosition(const ImageNoiseType & noise, const IndexType & index) const
{
DoubleVectorType curCoeffs;
if(m_UseFastCalibrationMethod)
{
curCoeffs = m_NoisePolynomialCoefficientsList[0];
}
else
// First compute the range step for the given noise record
double rangeStep = (noise.get_validityRangeMax()-noise.get_validityRangeMin())
/static_cast<double>(m_ImageSize[0]);
// Compute the range position using the rangeStep
double position = index[0] * rangeStep + noise.get_validityRangeMin();
return position;
}
template <class TInput, class TOutput>
double
TerraSarCalibrationImageFunctor<TInput, TOutput>
::ComputeNoiseEquivalentBetaNaught(const ImageNoiseType & noise, double range) const
{
// Formula: NEBN = Ks * SUM( coef_i * (tau - tau_ref)^i)
// Retrieve the polynomial degree
unsigned int polynomialDegree = noise.get_polynomialDegree();
// Compute tau - tau_ref
double deltaTau = range - noise.get_referencePoint();
// Get polynomial coefficients
std::vector<double> coefficients = noise.get_polynomialCoefficients();
// Initialize nebn value
double nebn = 0.;
// For each degree
for(unsigned int degree = 0; degree <= polynomialDegree; ++degree)
{
// (m_ImageSize[1]-1)-(lineId) because the first acquisition line is the last image one.
// (m_ImageSize[1]-1) because image starts to 0.
//double interval = static_cast<double>(m_ImageSize[1]) / static_cast<double>(m_NoisePolynomialCoefficientsList.size());
// compute utc time of the line
double currTimeUTC = m_TimeUTC[0] + static_cast<double>((m_ImageSize[1]-1)-lineId)*m_InvPRF;
unsigned int id = 0;
bool go = true;
// deduct the corresponding noise acquisition index
while( id<m_TimeUTC.size() && go)
{
if( currTimeUTC < m_TimeUTC[id] )
go = false;
id++;
}
id--;
double timeCoef = 1. / (m_TimeUTC[id]- m_TimeUTC[id-1]) * (currTimeUTC-m_TimeUTC[id-1]);
for(unsigned int j=0; j<m_NoisePolynomialCoefficientsList.size(); j++)
{
curCoeffs.push_back( m_NoisePolynomialCoefficientsList[id-1][j] + (m_NoisePolynomialCoefficientsList[id][j] - m_NoisePolynomialCoefficientsList[id-1][j]) * timeCoef );
}
// Cumulate polynomial
nebn += coefficients[degree] * vcl_pow(deltaTau,static_cast<double>(degree));
}
return curCoeffs;
}
// Do not forget to multiply by the calibration factor
nebn*=m_CalibrationFactor;
// Return the nebn value
return nebn;
}
template <class TInput, class TOutput>
......@@ -156,49 +188,60 @@ TOutput
TerraSarCalibrationImageFunctor<TInput, TOutput>
::operator()(const TInput & inPix, IndexType index)
{
double diffCurRange = ComputeCurrentNoise( static_cast<unsigned int>(index[0]) ) - this->GetNoiseRangeValidityRef();
// Formula: sigma = (Ks.|DN|²-NEBN) * sin Theta_local
DoubleVectorType curCoeff = this->ComputeCurrentCoeffs( static_cast<unsigned int>(index[1]) );
// First, compute the brightness using the brightness functor
double beta0 = m_Brightness( static_cast<double>(inPix) );
double outRadBr = this->GetBrightness()( static_cast<double>(inPix) );
// Then compute the current acquisition time
double currentAzimuth = this->ComputeAzimuthPosition(index);
double NEBN = 0.;
for(unsigned int i=0; i<curCoeff.size(); i++)
// Retrieve the associated noise record
ImageNoiseType currentNoiseRecord;
// If fast calibration is on, use always the first noise record
if(m_UseFastCalibrationMethod)
{
NEBN += curCoeff[i]*vcl_pow( diffCurRange, static_cast<double>(i));
currentNoiseRecord = m_NoiseRecords[0].second;
}
else
{
// Else, find the closest noise record
currentNoiseRecord = this->FindClosestNoiseRecord(currentAzimuth);
}
// Compute the range position for this noise record
double currentRange = this->ComputeRangePosition(currentNoiseRecord,index);
double sigma = ( outRadBr - this->GetCalFactor()*NEBN ) * this->GetSinLocalIncidentAngle();
// Compute the NEBN
double NEBN = this->ComputeNoiseEquivalentBetaNaught(currentNoiseRecord,currentRange);
// Last, apply formula
double sigma = ( beta0 - NEBN ) * m_SinLocalIncidentAngle;
// Cast and return
return static_cast<TOutput>(sigma);
}
template <class TInput, class TOutput>
std::complex<TOutput>
TerraSarCalibrationImageFunctor<TInput, TOutput>
::operator()(const std::complex<TInput> & inPix, IndexType index)
{
double diffCurRange = this->ComputeCurrentNoise( static_cast<unsigned int>(index[0]) ) - this->GetNoiseRangeValidityRef();
DoubleVectorType curCoeff = this->ComputeCurrentCoeffs( static_cast<unsigned int>(index[1]) );
// First, extract modulus and phase
double modulus = vcl_abs(inPix);
double outRadBr = static_cast<double>(this->GetBrightness()( modulus ));
double phase = std::arg(inPix);
double NEBN = 0.;
for(unsigned int i=0; i<curCoeff.size(); i++)
{
NEBN += curCoeff[i]*vcl_pow( diffCurRange, static_cast<double>(i));
}
double sigma = ( outRadBr - this->GetCalFactor()*NEBN ) * this->GetSinLocalIncidentAngle();
double phase = std::arg(inPix);
// Then, calibrate the modulus
double sigma = this->operator()(modulus,index);
// Last, put back the phase
std::complex<TOutput> out(sigma*vcl_cos(phase), sigma*vcl_sin(phase));
// And return the result
return out;
}
}// namespace Functor
} // namespace otb
......
......@@ -66,7 +66,7 @@ int otbTerraSarBrightnessImageFilterTest(int argc, char * argv[])
}
else
{
filter->SetCalFactor( 10 );
filter->SetCalibrationFactor( 10 );
writer->SetInput(filter->GetOutput());
}
......
......@@ -32,7 +32,7 @@ int otbTerraSarBrightnessImageFunctor(int argc, char * argv[])
typedef otb::Functor::TerraSarBrightnessImageFunctor<ScalarType, ScalarType> FunctorType;
FunctorType funct;
funct.SetCalFactor(calFact);
funct.SetCalibrationFactor(calFact);
ScalarType inPix = static_cast<ScalarType>(atof(argv[1]));
ComplexType inCplxPix(static_cast<ScalarType>(atof(argv[2])), static_cast<ScalarType>(atof(argv[3])));
......
......@@ -36,8 +36,6 @@ int otbTerraSarCalibrationImageFilterTest(int argc, char * argv[])
typedef otb::ImageFileReader<ImageType> ReaderType;
typedef otb::ImageFileWriter<ImageType> WriterType;
typedef otb::TerraSarCalibrationImageFilter<ImageType, ImageType> FilterType;
typedef FilterType::DoubleVectorType DoubleVectorType;
typedef FilterType::DoubleVectorVectorType DoubleVectorVectorType;
typedef itk::ExtractImageFilter<ImageType, ImageType> ExtractorType;
ReaderType::Pointer reader = ReaderType::New();
......@@ -51,68 +49,24 @@ int otbTerraSarCalibrationImageFilterTest(int argc, char * argv[])
reader->UpdateOutputInformation();
filter->SetInput(reader->GetOutput());
if( !useMetadata )
{
DoubleVectorType coefs;
coefs.push_back(1.);
coefs.push_back(0.5);
coefs.push_back(1.);
coefs.push_back(0.1);
DoubleVectorVectorType coefVect(1, coefs);
coefs.clear();
coefs.push_back(10.);
coefs.push_back(5);
coefs.push_back(10);
coefs.push_back(1);
coefVect.push_back( coefs );
writer->SetInput(filter->GetOutput());
coefs.clear();
coefs.push_back(100);
coefs.push_back(50);
coefs.push_back(100);
coefs.push_back(10);
coefVect.push_back( coefs );
filter->SetNoisePolynomialCoefficientsList(coefVect);
// Generate an extract from the large input
ImageType::RegionType region;
ImageType::IndexType id;
id[0] = 50; id[1] = 100;
ImageType::SizeType size;
size[0] = 150; size[1] = 100;
region.SetIndex(id);
region.SetSize(size);
extractor->SetExtractionRegion(region);
filter->SetCalFactor( 10 );
filter->SetNoiseRangeValidityMin( 0 );
filter->SetNoiseRangeValidityMax( 1 );
filter->SetNoiseRangeValidityRef( 0.5 );
filter->SetLocalIncidentAngle( 15 );
std::vector<double> timeUtc;
timeUtc.push_back(1.);
timeUtc.push_back(2.);
timeUtc.push_back(3.);
filter->SetTimeUTC(timeUtc);
filter->SetPRF(50);
extractor->SetInput(filter->GetOutput());
writer->SetInput(extractor->GetOutput());
writer->SetInput(filter->GetOutput());
}
else
{
// Generate an extract from the large input
ImageType::RegionType region;
ImageType::IndexType id;
id[0] = 50; id[1] = 100;
ImageType::SizeType size;
size[0] = 150; size[1] = 100;
region.SetIndex(id);
region.SetSize(size);
extractor->SetExtractionRegion(region);
extractor->SetInput(filter->GetOutput());
writer->SetInput(extractor->GetOutput());
}
filter->SetUseFastCalibrationMethod( false );
writer->Update();
return EXIT_SUCCESS;
}
......@@ -1202,6 +1202,7 @@ bool ossimplugins::ossimTerraSarProductDoc::initIncidenceAngles(
return result;
}
return false;
}
bool ossimplugins::ossimTerraSarProductDoc::initNoise(
const ossimXmlDocument* xdoc, ossimplugins::Noise* noise) const
......
......@@ -67,27 +67,63 @@ public:
{
_timeUTC = value;
}
const ossimString & get_timeUTC(void) const
{
return _timeUTC;
}
void set_validityRangeMin(double value)
{
_validityRangeMin = value;
}
double get_validityRangeMin(void) const
{
return _validityRangeMin;
}
void set_validityRangeMax(double value)
{
_validityRangeMax = value;
}
double get_validityRangeMax(void) const
{
return _validityRangeMax;
}
void set_referencePoint(double value)
{
_referencePoint = value;
}
double get_referencePoint(void) const
{
return _referencePoint;
}
void set_polynomialDegree(unsigned int value)
{
_polynomialDegree = value;
}
unsigned int get_polynomialDegree() const
{
return _polynomialDegree;
}
void set_polynomialCoefficients(const std::vector<double>& value)
{
_polynomialCoefficients = value;
}
const std::vector<double> & get_polynomialCoefficients(void) const
{
return _polynomialCoefficients;
}
protected:
/**
......
......@@ -76,6 +76,11 @@ public:
{
_tabImageNoise = image_noise;
}
const std::vector<ImageNoise> & get_imageNoise() const
{
return _tabImageNoise;
}
protected:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment