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

ENH: the norm of the input is now also a parameter of SpectralAngleDetails:ComputeSpectralAngle

parent 8c7f4cd3
...@@ -30,6 +30,9 @@ namespace otb ...@@ -30,6 +30,9 @@ namespace otb
/** \class BinarySpectralAngleFunctor /** \class BinarySpectralAngleFunctor
* \brief This functor computes the spectral angle between two pixels. * \brief This functor computes the spectral angle between two pixels.
* *
* If the pixels have different sizes, only the first components of the
* largest pixel will be considered.
*
* It can be used as a functor in a BinaryFunctorImageFilter to * It can be used as a functor in a BinaryFunctorImageFilter to
* compute the pixel-by-pixel spectral angles. * compute the pixel-by-pixel spectral angles.
* *
...@@ -46,9 +49,18 @@ public: ...@@ -46,9 +49,18 @@ public:
virtual ~BinarySpectralAngleFunctor() = default; virtual ~BinarySpectralAngleFunctor() = default;
// Binary operator // Binary operator
inline TOutputValue operator()(const TInput1& a, const TInput2& b) const inline TOutputValue operator()(const TInput1& in1, const TInput2& in2) const
{ {
return SpectralAngleDetails::ComputeSpectralAngle<TInput1, TInput2, TOutputValue>(a, b, b.GetNorm()); // Compute norms.
auto in1Norm = 0;
auto in2Norm = 0;
for (unsigned int i = 0; i < std::min(in1.Size(), in2.Size()); ++i)
{
in1Norm += in1[i] * in1[i];
in2Norm += in2[i] * in2[i];
}
return SpectralAngleDetails::ComputeSpectralAngle<TInput1, TInput2, TOutputValue>(in1, in1Norm, in2, in2Norm);
} }
}; };
......
...@@ -33,18 +33,16 @@ namespace SpectralAngleDetails ...@@ -33,18 +33,16 @@ namespace SpectralAngleDetails
{ {
template <class TInput, class TReference, class TOutput> template <class TInput, class TReference, class TOutput>
TOutput ComputeSpectralAngle(TInput const & input, TReference const & reference, TOutput ComputeSpectralAngle(TInput const & input, typename TInput ::ValueType const & inputNorm,
typename TReference::ValueType refNorm) TReference const & reference, typename TReference::ValueType refNorm)
{ {
// Compute norm and scalar product. // Compute scalar product.
double scalarProduct = 0.0; double scalarProduct = 0.0;
double squaredNormInput = 0.0;
for (unsigned int i = 0; i < std::min(input.Size(), reference.Size()); ++i) for (unsigned int i = 0; i < std::min(input.Size(), reference.Size()); ++i)
{ {
scalarProduct += input[i] * reference[i]; scalarProduct += input[i] * reference[i];
squaredNormInput += input[i] * input[i];
} }
auto normProd = std::sqrt(squaredNormInput) * refNorm; auto normProd = inputNorm * refNorm;
if ((normProd == 0.0) || (scalarProduct / normProd > 1)) if ((normProd == 0.0) || (scalarProduct / normProd > 1))
{ {
return static_cast<TOutput>(0.0); return static_cast<TOutput>(0.0);
...@@ -77,7 +75,7 @@ public: ...@@ -77,7 +75,7 @@ public:
// Binary operator // Binary operator
inline TOutputValue operator()(TInput const & inPix) const inline TOutputValue operator()(TInput const & inPix) const
{ {
return SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, TOutputValue>(inPix, m_ReferencePixel, m_RefNorm); return SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, TOutputValue>(inPix, inPix.GetNorm(), m_ReferencePixel, m_RefNorm);
} }
void SetReferencePixel(TInput const & ref) void SetReferencePixel(TInput const & ref)
...@@ -109,10 +107,12 @@ public: ...@@ -109,10 +107,12 @@ public:
TOutput res; TOutput res;
res.SetSize(m_ReferencePixels.size()); res.SetSize(m_ReferencePixels.size());
auto inputNorm = inPix.GetNorm();
for (unsigned int i = 0; i< m_ReferencePixels.size(); i++) for (unsigned int i = 0; i< m_ReferencePixels.size(); i++)
{ {
res[i] = SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, typename TOutput::ValueType> res[i] = SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, typename TOutput::ValueType>
(inPix, m_ReferencePixels[i], m_ReferenceNorm[i]); (inPix, inputNorm, m_ReferencePixels[i], m_ReferenceNorm[i]);
} }
return res; return res;
......
...@@ -56,7 +56,7 @@ public: ...@@ -56,7 +56,7 @@ public:
// Binary operator // Binary operator
inline TOutputValue operator()(TInput const & inPix) const inline TOutputValue operator()(TInput const & inPix) const
{ {
return std::sqrt(SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, TOutputValue>(inPix, m_ReferencePixel, m_RefNorm)); return std::sqrt(SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, TOutputValue>(inPix, inPix.GetNorm(), m_ReferencePixel, m_RefNorm));
} }
void SetReferencePixel(TInput const & ref) void SetReferencePixel(TInput const & ref)
......
...@@ -82,9 +82,8 @@ public: ...@@ -82,9 +82,8 @@ public:
pix[3] = this->Value(CommonBandNames::NIR, inPix); pix[3] = this->Value(CommonBandNames::NIR, inPix);
return std::sqrt(SpectralAngleDetails::ComputeSpectralAngle<PixelType, PixelType, TOutput> return std::sqrt(SpectralAngleDetails::ComputeSpectralAngle<PixelType, PixelType, TOutput>
(pix, (pix, pix.GetNorm(),
m_ReferencePixel, m_ReferencePixel, m_RefNorm));
m_RefNorm));
} }
/**@{*/ /**@{*/
......
...@@ -25,7 +25,7 @@ ...@@ -25,7 +25,7 @@
#include "otbParser.h" #include "otbParser.h"
#include "otbMacro.h" #include "otbMacro.h"
#include "otbSpectralAngleFunctor.h" #include "otbBinarySpectralAngleFunctor.h"
#include <vnl/algo/vnl_lsqr.h> #include <vnl/algo/vnl_lsqr.h>
#include <vnl/vnl_sparse_matrix_linear_system.h> #include <vnl/vnl_sparse_matrix_linear_system.h>
...@@ -115,8 +115,9 @@ public: ...@@ -115,8 +115,9 @@ public:
m_IntensityP2 = m_IntensityP2 / (static_cast<double>(m_NbOfBands)); m_IntensityP2 = m_IntensityP2 / (static_cast<double>(m_NbOfBands));
m_Distance = std::sqrt(m_Distance); m_Distance = std::sqrt(m_Distance);
m_SpectralAngle = SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, double>(p1, p2, p2.GetNorm()); BinarySpectralAngleFunctor<TInput, TInput, double> spectralAngleFunctor;
m_SpectralAngle = spectralAngleFunctor(p1, p2);
value = m_Parser->Eval(); value = m_Parser->Eval();
......
Supports Markdown
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