Commit 0290d845 authored by Cyrille Valladeau's avatar Cyrille Valladeau
Browse files

ENH : change TSX brightness and calib functors : 1 functor for scalar and complex

parent 0bd9d86c
......@@ -24,7 +24,7 @@
#include "otbUnaryFunctorWithIndexImageFilter.h"
#include "otbRadarFunctors.h"
#include "otbTerraSarFunctors.h"
#include "itkMetaDataDictionary.h"
//#include "itkConstNeighborhoodIterator.h"
#include "otbMath.h"
......
......@@ -19,8 +19,8 @@
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __otbRadarFunctors_h
#define __otbRadarFunctors_h
#ifndef __otbTerraSarFunctors_h
#define __otbTerraSarFunctors_h
#include "itkUnaryFunctorImageFilter.h"
......@@ -31,6 +31,8 @@ namespace otb
{
namespace Functor
{
/**
* \class TerraSarRadarBrightnessImageFunctor
* \brief Compute the radar brightness from an modulus image.
......@@ -42,29 +44,21 @@ template <class TInput, class TOutput>
class TerraSarRadarBrightnessImageFunctor
{
public:
TerraSarRadarBrightnessImageFunctor()
{
m_CalFactor = 1.;
};
virtual ~TerraSarRadarBrightnessImageFunctor() {};
typedef std::vector<double> DoubleVectorType;
typedef std::vector<DoubleVectorType> DoubleVectorVectorType;
typedef itk::Size<2> SizeType;
/** Accessors */
void SetCalFactor( double val ) { m_CalFactor = val; };
double GetCalFactor() { return m_CalFactor; };
inline TOutput operator() (const TInput & inPix)
{
double squareInPix = vcl_pow( static_cast<double>(inPix), 2.);
// Beta naught computation
double beta = m_CalFactor*squareInPix;
return static_cast<TOutput>(beta);
}
TerraSarRadarBrightnessImageFunctor();
virtual ~TerraSarRadarBrightnessImageFunctor() {};
typedef std::vector<double> DoubleVectorType;
typedef std::vector<DoubleVectorType> DoubleVectorVectorType;
typedef itk::Size<2> SizeType;
/** Accessors */
void SetCalFactor( double val ) { m_CalFactor = val; };
double GetCalFactor() { return m_CalFactor; };
/** We assume that the input pixel is a scalar -> modulus image */
inline TOutput operator() (const TInput & inPix);
/** We assume that the input pixel is a complex -> complex image */
inline std::complex<TOutput> operator() (const std::complex<TInput> & inPix);
private:
/** Calibration Factor */
......@@ -73,48 +67,6 @@ private:
/**
* \class TerraSarRadarBrightnessComplexImageFunctor
* \brief Compute the radar brightness from an complexe image.
*
* \ingroup Functor
* \ingroup Radiometry
*/
template <class TInput, class TOutput>
class TerraSarRadarBrightnessComplexImageFunctor
{
public:
TerraSarRadarBrightnessComplexImageFunctor() {};
virtual ~TerraSarRadarBrightnessComplexImageFunctor() {};
typedef TerraSarRadarBrightnessImageFunctor<double, double> BetaNaughtFunctorType;
/** Accessors */
void SetCalFactor( double val ) { m_BetaNaughtFunctor.SetCalFactor(val); };
double GetCalFactor() { return m_BetaNaughtFunctor.GetCalFactor(); };
/* We assume that the input pixel is a complex */
inline TOutput operator() (const TInput & inPix)
{
// Beta naught computation, will be the Modulus of the result
double beta = m_BetaNaughtFunctor(static_cast<double>(std::abs(inPix)));
// Phase
double phase = std::arg(inPix);
// We retrieve the complex value from the modulus and the phase.
std::complex<double> res = std::complex<double>(beta*vcl_cos(phase), beta*vcl_sin(phase) );
return static_cast<TOutput>(res);
}
private:
/** Calibration Factor */
BetaNaughtFunctorType m_BetaNaughtFunctor;
};
/**
* \class TerraSarCalibrationImageFunctor
......@@ -130,16 +82,15 @@ public:
TerraSarCalibrationImageFunctor();
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 typename TInputIt::PixelType InputPixelType;
typedef TerraSarRadarBrightnessImageFunctor<TInput, TOutput> BrightnessFunctorType;
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 TerraSarRadarBrightnessImageFunctor<double, double> BrightnessFunctorType;
/** Accessors */
void SetCalFactor( double val ) { m_CalFactor = val; m_RadarBrightness.SetCalFactor(val); };
void SetCalFactor( double val ) { m_CalFactor = val; m_RadarBrightness.SetCalFactor(val); };
double GetCalFactor() const { return m_CalFactor; };
void SetNoiseRangeValidityMin( double val ) { m_NoiseRangeValidityMin = val; };
double GetNoiseRangeValidityMin() const { return m_NoiseRangeValidityMin; };
......@@ -165,13 +116,14 @@ public:
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 GetRadarBrightness() const { return m_RadarBrightness; };
//BrightnessFunctorType GetRadarBrightness() { return m_RadarBrightness; };
BrightnessFunctorType GetRadarBrightness() { return m_RadarBrightness; };
double ComputeCurrentNoise( unsigned int colId );
DoubleVectorType ComputeCurrentCoeffs( unsigned int lineId );
inline TOutput operator() (const TInput & inPix, IndexType index);
/** We assume that the input pixel is a scalar -> modulus image */
inline TOutput operator() (const TInput & inPix, IndexType index);
/** We assume that the input pixel is a complex -> complex image */
inline std::complex<TOutput> operator() (const std::complex<TInput> & inPix, IndexType index);
private:
/** Calibration Factor */
......@@ -200,106 +152,17 @@ private:
double m_PRF;
/** Inverse Pulse Repetition Frequency */
double m_InvPRF;
/** Radar Brightness functor */
/** Radar Brightness functor */
BrightnessFunctorType m_RadarBrightness;
};
/** TODO : Use inheritance **/
/**
* \class TerraSarCalibrationComplexImageFunctor
* \brief Compute sigma naught coefficient from a modulus image.
*
* \ingroup Functor
* \ingroup Radiometry
*/
template <class TInput, class TOutput>
class TerraSarCalibrationComplexImageFunctor
{
public:
TerraSarCalibrationComplexImageFunctor();
virtual ~TerraSarCalibrationComplexImageFunctor() {};
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 typename TInputIt::PixelType InputPixelType;
typedef TerraSarRadarBrightnessImageFunctor<double, double> BrightnessFunctorType;
/** Accessors */
void SetCalFactor( double val ) { m_CalFactor = val; m_RadarBrightness.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; };
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( LIntVectorType vect ) { m_TimeUTC = vect; };
LIntVectorType 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 GetRadarBrightness() const { return m_RadarBrightness; };
//BrightnessFunctorType GetRadarBrightness() { return m_RadarBrightness; };
double ComputeCurrentNoise( unsigned int colId );
DoubleVectorType ComputeCurrentCoeffs( unsigned int lineId );
inline TOutput operator() (const TInput & inPix, IndexType index);
private:
/** 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;
/** 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,
* 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 second). */
LIntVectorType m_TimeUTC;
/** Pulse Repetition Frequency */
double m_PRF;
/** Inverse Pulse Repetition Frequency */
double m_InvPRF;
/** Radar Brightness functor */
BrightnessFunctorType m_RadarBrightness;
};
}// end namespace functor
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbRadarFunctors.txx"
#include "otbTerraSarFunctors.txx"
#endif
#endif
......@@ -15,20 +15,61 @@
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __otbRadarFunctors_txx
#define __otbRadarFunctors_txx
#ifndef __otbTerraSarFunctors_txx
#define __otbTerraSarFunctors_txx
#include "otbRadarFunctors.h"
#include "otbTerraSarFunctors.h"
namespace otb
{
/**************************************************************************
***************** TerraSarCalibrationImageFunctor *******************
**************************************************************************/
namespace Functor
{
/** Constructor */
/******************************************************************/
/*********** TerraSarRadarBrightnessImageFunctor *****************/
/******************************************************************/
template <class TInput, class TOutput>
TerraSarRadarBrightnessImageFunctor<TInput, TOutput>
::TerraSarRadarBrightnessImageFunctor()
{
m_CalFactor = 1.;
}
template <class TInput, class TOutput>
TOutput
TerraSarRadarBrightnessImageFunctor<TInput, TOutput>
::operator() (const TInput & inPix)
{
double squareInPix = vcl_pow( static_cast<double>(inPix), 2.);
// Beta naught computation
double beta = m_CalFactor*squareInPix;
return static_cast<TOutput>(beta);
}
template <class TInput, class TOutput>
std::complex<TOutput>
TerraSarRadarBrightnessImageFunctor<TInput, TOutput>
::operator() (const std::complex<TInput> & inPix)
{
// Beta naught computation, will be the Modulus of the result
double beta = operator()(static_cast<double>(std::abs(inPix)));
// Phase
double phase = std::arg(inPix);
// We retrieve the complex value from the modulus and the phase.
std::complex<TOutput> res = std::complex<TOutput>(beta*vcl_cos(phase), beta*vcl_sin(phase) );
return static_cast<TOutput>(res);
}
/******************************************************************/
/*********** TerraSarCalibrationImageFunctor ****************/
/******************************************************************/
template <class TInput, class TOutput>
TerraSarCalibrationImageFunctor<TInput, TOutput>
::TerraSarCalibrationImageFunctor()
......@@ -38,6 +79,7 @@ TerraSarCalibrationImageFunctor<TInput, TOutput>
m_NoiseRangeValidityMax = 0.;
m_NoiseRangeValidityRef = 0.;
m_LocalIncidentAngle = 0.;
m_SinLocalIncidentAngle = 0.;
m_NoisePolynomialCoefficientsList.clear();
m_ImageSize.Fill(0);
m_UseFastCalibrationMethod = true;
......@@ -50,187 +92,110 @@ template <class TInput, class TOutput>
double
TerraSarCalibrationImageFunctor<TInput, TOutput>
::ComputeCurrentNoise( unsigned int colId )
{
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) )
{
curRange = m_NoiseRangeValidityMin + ( m_NoiseRangeValidityRef-m_NoiseRangeValidityMin )/width_2 * static_cast<double>(colId+1);
}
else
{
curRange = m_NoiseRangeValidityRef + ( m_NoiseRangeValidityMax-m_NoiseRangeValidityRef )/width_2 * (static_cast<double>(colId+1) - width_2 );
}
return curRange;
}
{
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) )
{
curRange = m_NoiseRangeValidityMin + ( m_NoiseRangeValidityRef-m_NoiseRangeValidityMin )/width_2 * static_cast<double>(colId+1);
}
else
{
curRange = m_NoiseRangeValidityRef + ( m_NoiseRangeValidityMax-m_NoiseRangeValidityRef )/width_2 * (static_cast<double>(colId+1) - width_2 );
}
return curRange;
}
template <class TInput, class TOutput>
typename TerraSarCalibrationImageFunctor<TInput, TOutput>::DoubleVectorType
TerraSarCalibrationImageFunctor<TInput, TOutput>
::ComputeCurrentCoeffs( unsigned int lineId )
{
DoubleVectorType curCoeffs;
if(m_UseFastCalibrationMethod)
{
curCoeffs = m_NoisePolynomialCoefficientsList[0];
}
else
{
// m_ImageSize[1]-(lineId+1) because the first acquisition line is the last image one.
// line+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]-(lineId-1))*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 );
}
{
DoubleVectorType curCoeffs;
if(m_UseFastCalibrationMethod)
{
curCoeffs = m_NoisePolynomialCoefficientsList[0];
}
else
{
// m_ImageSize[1]-(lineId+1) because the first acquisition line is the last image one.
// line+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]-(lineId-1))*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 );
}
}
return curCoeffs;
}
return curCoeffs;
}
template <class TInput, class TOutput>
TOutput
TerraSarCalibrationImageFunctor<TInput, TOutput>
::operator()(const TInput & inPix, IndexType index)
{
double diffCurRange = ComputeCurrentNoise( static_cast<unsigned int>(index[0]) ) - m_NoiseRangeValidityRef;
DoubleVectorType curCoeff = ComputeCurrentCoeffs( static_cast<unsigned int>(index[1]) );
double diffCurRange = ComputeCurrentNoise( static_cast<unsigned int>(index[0]) ) - this->GetNoiseRangeValidityRef();
DoubleVectorType curCoeff = this->ComputeCurrentCoeffs( static_cast<unsigned int>(index[1]) );
TOutput outRadBr = m_RadarBrightness( inPix );
double outRadBr = this->GetRadarBrightness()( static_cast<double>(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 - m_CalFactor*NEBN ) * m_SinLocalIncidentAngle;
double sigma = ( outRadBr - this->GetCalFactor()*NEBN ) * this->GetSinLocalIncidentAngle();
return static_cast<TOutput>(sigma);
}
/** Constructor */
template <class TInput, class TOutput>
TerraSarCalibrationComplexImageFunctor<TInput, TOutput>
::TerraSarCalibrationComplexImageFunctor()
{
m_CalFactor = 1.;
m_NoiseRangeValidityMin = 0.;
m_NoiseRangeValidityMax = 0.;
m_NoiseRangeValidityRef = 0.;
m_LocalIncidentAngle = 0.;
m_NoisePolynomialCoefficientsList.clear();
m_ImageSize.Fill(0);
m_UseFastCalibrationMethod = true;
m_TimeUTC.clear();
m_PRF = 1.;
}
template <class TInput, class TOutput>
double
TerraSarCalibrationComplexImageFunctor<TInput, TOutput>
::ComputeCurrentNoise( unsigned int colId )
{
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) )
{
curRange = m_NoiseRangeValidityMin + ( m_NoiseRangeValidityRef-m_NoiseRangeValidityMin )/width_2 * static_cast<double>(colId+1);
}
else
{
curRange = m_NoiseRangeValidityRef + ( m_NoiseRangeValidityMax-m_NoiseRangeValidityRef )/width_2 * (static_cast<double>(colId+1) - width_2 );
}
return curRange;
}
template <class TInput, class TOutput>
typename TerraSarCalibrationComplexImageFunctor<TInput, TOutput>::DoubleVectorType
TerraSarCalibrationComplexImageFunctor<TInput, TOutput>
::ComputeCurrentCoeffs( unsigned int lineId )
{
DoubleVectorType curCoeffs;
if(m_UseFastCalibrationMethod)
{
curCoeffs = m_NoisePolynomialCoefficientsList[0];
}
else
{
// m_ImageSize[1]-(lineId+1) because the first acquisition line is the last image one.
// line+1 because image starts to 0.
// compute utc time of the line
double currTimeUTC = m_TimeUTC[0] + static_cast<double>(m_ImageSize[1]-(lineId-1))*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 );
}
}
return curCoeffs;
}
template <class TInput, class TOutput>
TOutput
TerraSarCalibrationComplexImageFunctor<TInput, TOutput>
::operator()(const TInput & inPix, IndexType index)
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]) );
double modulus = std::abs(inPix);
double outRadBr = static_cast<double>(m_RadarBrightness( modulus ));
double outRadBr = static_cast<double>(this->GetRadarBrightness()( modulus ));
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);
TOutput out(sigma*vcl_cos(phase), sigma*vcl_sin(phase));
std::complex<TOutput> out(sigma*vcl_cos(phase), sigma*vcl_sin(phase));
return out;
}
}// namespace Functor
} // namespace otb
......
......@@ -26,7 +26,7 @@
#include "itkUnaryFunctorImageFilter.h"
#include "itkMetaDataDictionary.h"
#include "otbMath.h"
#include "otbRadarFunctors.h"
#include "otbTerraSarFunctors.h"
namespace otb
{
......
......@@ -1056,38 +1056,8 @@ ADD_TEST(raTvTerraSarRadarBrightnessImageFilterTest ${RADIOMETRY_TESTS9}
${TEMP}/raTvTerraSarRadarBrightnessImageFilterTest.tif
)
# ------- TerraSarRadarBrightnessComplexImageFilter ------------------------------
ADD_TEST(raTuTerraSarRadarBrightnessComplexImageFilterNew ${RADIOMETRY_TESTS9}