Commit c93be2c6 authored by Cédric Traizet's avatar Cédric Traizet
Browse files

REFAC: factorize spectral angle code in the library

parent 544cb182
Pipeline #4097 passed with stages
in 10 minutes and 10 seconds
......@@ -23,6 +23,7 @@
#include <algorithm>
#include "otbMath.h"
#include "otbSpectralAngleFunctor.h"
namespace otb
{
......@@ -40,40 +41,14 @@ template <class TInput1, class TInput2, class TOutputValue>
class BinarySpectralAngleFunctor
{
public:
BinarySpectralAngleFunctor()
{
}
BinarySpectralAngleFunctor() = default;
virtual ~BinarySpectralAngleFunctor()
{
}
virtual ~BinarySpectralAngleFunctor() = default;
// Binary operator
inline TOutputValue operator()(const TInput1& a, const TInput2& b) const
{
const double Epsilon = 1E-10;
double dist = 0.0;
double scalarProd = 0.0;
double norma = 0.0;
double normb = 0.0;
double sqrtNormProd = 0.0;
for (unsigned int i = 0; i < std::min(a.Size(), b.Size()); ++i)
{
scalarProd += a[i] * b[i];
norma += a[i] * a[i];
normb += b[i] * b[i];
}
sqrtNormProd = std::sqrt(norma * normb);
if (std::abs(sqrtNormProd) < Epsilon || scalarProd / sqrtNormProd > 1)
{
dist = 0.0;
}
else
{
dist = std::acos(scalarProd / sqrtNormProd);
}
return static_cast<TOutputValue>(dist);
return SpectralAngleDetails::ComputeSpectralAngle<TInput1, TInput2, TOutputValue>(a, b, b.GetNorm());
}
};
......
......@@ -26,13 +26,42 @@
namespace otb
{
namespace Functor
{
namespace SpectralAngleDetails
{
template <class TInput, class TReference, class TOutput>
TOutput ComputeSpectralAngle(TInput const & input, TReference const & reference,
typename TReference::ValueType refNorm)
{
// Compute norm and scalar product.
double scalarProduct = 0.0;
double squaredNormInput = 0.0;
for (unsigned int i = 0; i < std::min(input.Size(), reference.Size()); ++i)
{
scalarProduct += input[i] * reference[i];
squaredNormInput += input[i] * input[i];
}
auto normProd = std::sqrt(squaredNormInput) * refNorm;
if ((normProd == 0.0) || (scalarProduct / normProd > 1))
{
return static_cast<TOutput>(0.0);
}
else
{
return static_cast<TOutput>(std::acos(scalarProduct / normProd));
}
}
} // end namespace SpectralAngleDetails
/** \class SpectralAngleFunctor
* \brief This functor computes the spectral angle according to a reference pixel.
*
* \ingroup OTBImageManipulation
*/
namespace Functor
{
template <class TInput, class TOutputValue>
class SpectralAngleFunctor
{
......@@ -43,66 +72,34 @@ public:
m_ReferencePixel.Fill(1);
}
virtual ~SpectralAngleFunctor()
{
}
virtual ~SpectralAngleFunctor() = default;
// Binary operator
inline TOutputValue operator()(const TInput& inPix) const
inline TOutputValue operator()(TInput const & inPix) const
{
return this->Evaluate(inPix);
return SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, TOutputValue>(inPix, m_ReferencePixel, m_RefNorm);
}
void SetReferencePixel(TInput ref)
void SetReferencePixel(TInput const & ref)
{
m_ReferencePixel = ref;
m_RefNorm = 0.0;
for (unsigned int i = 0; i < ref.Size(); ++i)
{
m_RefNorm += ref[i] * ref[i];
}
m_RefNorm = std::sqrt(static_cast<double>(m_RefNorm));
m_RefNorm = ref.GetNorm();
}
TInput GetReferencePixel() const
{
return m_ReferencePixel;
}
protected:
// This method can be reimplemented in subclasses to actually
// compute the index value
virtual TOutputValue Evaluate(const TInput& inPix) const
{
TOutputValue out;
double dist = 0.0;
double scalarProd = 0.0;
double normProd = 0.0;
double normProd1 = 0.0;
double sqrtNormProd = 0.0;
for (unsigned int i = 0; i < std::min(inPix.Size(), m_ReferencePixel.Size()); ++i)
{
scalarProd += inPix[i] * m_ReferencePixel[i];
normProd1 += inPix[i] * inPix[i];
}
normProd = normProd1 * m_RefNorm * m_RefNorm;
sqrtNormProd = std::sqrt(normProd);
if ((sqrtNormProd == 0.0) || (scalarProd / sqrtNormProd > 1))
{
dist = 0.0;
}
else
{
dist = std::acos(scalarProd / sqrtNormProd);
}
out = static_cast<TOutputValue>(dist);
return out;
}
private :
TInput m_ReferencePixel;
double m_RefNorm;
};
} // end namespace functor
} // end namespace otb
......
......@@ -35,25 +35,44 @@ namespace Functor
*
* \ingroup OTBImageManipulation
*/
template <class TInputVectorPixel, class TOutputPixel>
class SqrtSpectralAngleFunctor : public SpectralAngleFunctor<TInputVectorPixel, TOutputPixel>
/** \class SpectralAngleFunctor
* \brief This functor computes the spectral angle according to a reference pixel.
*
* \ingroup OTBImageManipulation
*/
template <class TInput, class TOutputValue>
class SqrtSpectralAngleFunctor
{
public:
typedef SqrtSpectralAngleFunctor Self;
typedef SpectralAngleFunctor<TInputVectorPixel, TOutputPixel> Superclass;
SqrtSpectralAngleFunctor()
{
m_ReferencePixel.SetSize(4);
m_ReferencePixel.Fill(1);
}
~SqrtSpectralAngleFunctor() override
virtual ~SqrtSpectralAngleFunctor() = default;
// Binary operator
inline TOutputValue operator()(TInput const & inPix) const
{
return std::sqrt(SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, TOutputValue>(inPix, m_ReferencePixel, m_RefNorm));
}
protected:
TOutputPixel Evaluate(const TInputVectorPixel& inPix) const override
void SetReferencePixel(TInput const & ref)
{
return static_cast<TOutputPixel>(std::sqrt(Superclass::Evaluate(inPix)));
m_ReferencePixel = ref;
m_RefNorm = ref.GetNorm();
}
TInput GetReferencePixel() const
{
return m_ReferencePixel;
}
private :
TInput m_ReferencePixel;
double m_RefNorm;
};
} // end namespace Functor
......
......@@ -23,6 +23,7 @@
#include "otbSqrtSpectralAngleFunctor.h"
#include "itkUnaryFunctorImageFilter.h"
#include "otbRadiometricIndex.h"
namespace otb
{
......@@ -38,100 +39,103 @@ namespace Functor
*
* \ingroup OTBIndices
*/
template <class TInputVectorPixel, class TOutputPixel>
class WaterSqrtSpectralAngleFunctor : public SqrtSpectralAngleFunctor<TInputVectorPixel, TOutputPixel>
template <class TInput, class TOutput>
class WaterSqrtSpectralAngleFunctor : public RadiometricIndex<TInput, TOutput>
{
public:
typedef WaterSqrtSpectralAngleFunctor Self;
typedef SqrtSpectralAngleFunctor<TInputVectorPixel, TOutputPixel> Superclass;
typedef TInputVectorPixel InputVectorPixelType;
WaterSqrtSpectralAngleFunctor()
{
using typename RadiometricIndex<TInput, TOutput>::PixelType;
// Set the channels indices
m_BlueIndex = 0;
m_GreenIndex = 1;
m_RedIndex = 2;
m_NIRIndex = 3;
// Set reference water value
InputVectorPixelType reference;
reference.SetSize(4);
reference[0] = 136.0;
reference[1] = 132.0;
reference[2] = 47.0;
reference[3] = 24.0;
this->SetReferenceWaterPixel(reference);
WaterSqrtSpectralAngleFunctor()
: RadiometricIndex<TInput, TOutput>({CommonBandNames::BLUE, CommonBandNames::GREEN, CommonBandNames::RED, CommonBandNames::NIR} )
{
// Set default reference water value
m_ReferencePixel.SetSize(4);
m_ReferencePixel[0] = 136.0;
m_ReferencePixel[1] = 132.0;
m_ReferencePixel[2] = 47.0;
m_ReferencePixel[3] = 24.0;
m_RefNorm = m_ReferencePixel.GetNorm();
}
~WaterSqrtSpectralAngleFunctor() override
virtual ~WaterSqrtSpectralAngleFunctor() = default;
/** Set Reference Pixel
* \param ref : The new reference pixel, the band indices will be used.
* */
void SetReferenceWaterPixel(PixelType ref)
{
assert(m_ReferencePixel.Size == 4);
m_ReferencePixel[0] = this->Value(CommonBandNames::BLUE, ref);
m_ReferencePixel[1] = this->Value(CommonBandNames::GREEN, ref);
m_ReferencePixel[2] = this->Value(CommonBandNames::RED, ref);
m_ReferencePixel[3] = this->Value(CommonBandNames::NIR, ref);
m_RefNorm = m_ReferencePixel.GetNorm();
}
/** Set Reference Pixel */
void SetReferenceWaterPixel(InputVectorPixelType ref)
// Binary operator
inline TOutput operator()(PixelType const & inPix) const override
{
if (ref.GetSize() != 4)
{
}
InputVectorPixelType reference;
reference.SetSize(4);
reference[m_BlueIndex] = ref[0];
reference[m_GreenIndex] = ref[1];
reference[m_RedIndex] = ref[2];
reference[m_NIRIndex] = ref[3];
this->SetReferencePixel(reference);
PixelType pix(4);
pix[0] = this->Value(CommonBandNames::BLUE, inPix);
pix[1] = this->Value(CommonBandNames::GREEN, inPix);
pix[2] = this->Value(CommonBandNames::RED, inPix);
pix[3] = this->Value(CommonBandNames::NIR, inPix);
return std::sqrt(SpectralAngleDetails::ComputeSpectralAngle<PixelType, PixelType, TOutput>
(pix,
m_ReferencePixel,
m_RefNorm));
}
/** Getters and setters */
/**@{*/
/** Legacy getters and setters (use SetBandIndex/GetBandIndex instead)
* \deprecated
* */
void SetBlueChannel(unsigned int channel)
{
m_BlueIndex = channel;
this->SetBandIndex(CommonBandNames::BLUE, channel + 1);
}
unsigned int GetBlueChannel() const
{
return m_BlueIndex;
return this->GetBandIndex(CommonBandNames::BLUE) - 1;
}
void SetGreenChannel(unsigned int channel)
{
m_GreenIndex = channel;
this->SetBandIndex(CommonBandNames::GREEN, channel + 1);
}
unsigned int GetGreenChannel() const
{
return m_GreenIndex;
return this->GetBandIndex(CommonBandNames::GREEN) - 1;
}
void SetRedChannel(unsigned int channel)
{
m_RedIndex = channel;
this->SetBandIndex(CommonBandNames::RED, channel + 1);
}
unsigned int GetRedChannel() const
{
return m_RedIndex;
return this->GetBandIndex(CommonBandNames::RED) - 1;
}
void SetNIRChannel(unsigned int channel)
{
m_NIRIndex = channel;
this->SetBandIndex(CommonBandNames::NIR, channel + 1);
}
unsigned int GetNIRChannel() const
{
return m_NIRIndex;
return this->GetBandIndex(CommonBandNames::NIR) - 1;
}
/**@}*/
protected:
inline TOutputPixel Evaluate(const TInputVectorPixel& inPix) const override
{
return static_cast<TOutputPixel>(Superclass::Evaluate(inPix));
}
private:
/** Channels */
int m_BlueIndex;
int m_GreenIndex;
int m_RedIndex;
int m_NIRIndex;
PixelType m_ReferencePixel;
double m_RefNorm;
};
} // End namespace Functor
/** \class WaterSqrtSpectralAngleImageFilter
* \deprecated
* \brief Compute a radiometric water indice
*
* This filter calculates a pixel wise water indice by calculating
......@@ -147,11 +151,11 @@ protected:
* probability of water.
*
* \sa SpectralAngleDistanceImageFilter
*
*
* \ingroup OTBIndices
*/
template <class TInputVectorImage, class TOutputImage,
class TFunction = Functor::WaterSqrtSpectralAngleFunctor<typename TInputVectorImage::PixelType, typename TOutputImage::PixelType>>
class TFunction = Functor::WaterSqrtSpectralAngleFunctor<double, typename TOutputImage::PixelType>>
class ITK_EXPORT WaterSqrtSpectralAngleImageFilter : public itk::UnaryFunctorImageFilter<TInputVectorImage, TOutputImage, TFunction>
{
public:
......
......@@ -39,6 +39,7 @@ int otbWaterSqrtSpectralAngleImageFilter(int itkNotUsed(argc), char* argv[])
// Instantiating objects
WaterSqrtSpectralAngleImageFilterType::Pointer filter = WaterSqrtSpectralAngleImageFilterType::New();
ReaderType::Pointer reader = ReaderType::New();
WriterType::Pointer writer = WriterType::New();
......@@ -60,7 +61,7 @@ int otbWaterSqrtSpectralAngleImageFilter(int itkNotUsed(argc), char* argv[])
filter->SetInput(reader->GetOutput());
writer->SetInput(filter->GetOutput());
writer->Update();
return EXIT_SUCCESS;
}
......@@ -25,7 +25,7 @@
#include "otbParser.h"
#include "otbMacro.h"
#include "otbSpectralAngleFunctor.h"
#include <vnl/algo/vnl_lsqr.h>
#include <vnl/vnl_sparse_matrix_linear_system.h>
......@@ -116,30 +116,7 @@ public:
m_Distance = std::sqrt(m_Distance);
// compute spectralAngle
double scalarProd = 0.0;
double normProd = 0.0;
double normProd1 = 0.0;
double normProd2 = 0.0;
for (unsigned int i = 0; i < p1.Size(); ++i)
{
scalarProd += p1[i] * p2[i];
normProd1 += p1[i] * p1[i];
normProd2 += p2[i] * p2[i];
}
normProd = normProd1 * normProd2;
if (normProd == 0.0)
{
m_SpectralAngle = 0.0;
}
else
{
m_SpectralAngle = std::acos(scalarProd / std::sqrt(normProd));
}
//
m_SpectralAngle = SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, double>(p1, p2, p2.GetNorm());
value = m_Parser->Eval();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment