From dcafc1bdcaafd036f5d5f857b66eb911153ff42b Mon Sep 17 00:00:00 2001
From: Guillaume Borrut <guillaume.borrut@c-s.fr>
Date: Thu, 29 Jan 2009 18:30:51 +0100
Subject: [PATCH] BUG : SpectralAngleFunctor, arcos(x) gives us a nan value if
x>1
---
.../otbSpectralAngleFunctor.h | 29 ++++++++++---------
1 file changed, 15 insertions(+), 14 deletions(-)
diff --git a/Code/FeatureExtraction/otbSpectralAngleFunctor.h b/Code/FeatureExtraction/otbSpectralAngleFunctor.h
index 5b724980f4..a234f6e710 100644
--- a/Code/FeatureExtraction/otbSpectralAngleFunctor.h
+++ b/Code/FeatureExtraction/otbSpectralAngleFunctor.h
@@ -40,28 +40,29 @@ class SpectralAngleFunctor
~SpectralAngleFunctor(){};
inline TOutputValue operator()(const TInput& inPix)
{
- TOutputValue out;
+ 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];
- }
+ {
+ scalarProd += inPix[i]*m_ReferencePixel[i];
+ normProd1 += inPix[i]*inPix[i];
+ }
normProd = normProd1 * m_RefNorm*m_RefNorm;
-
- if ( normProd == 0.0)
- {
- dist = 0.0;
- }
+ sqrtNormProd = vcl_sqrt(normProd);
+ if ( (sqrtNormProd == 0.0) || ( scalarProd / sqrtNormProd > 1) )
+ {
+ dist = 0.0;
+ }
else
- {
- dist = vcl_acos(scalarProd/vcl_sqrt(normProd));
- }
-
+ {
+ dist = vcl_acos(scalarProd/sqrtNormProd);
+ }
+
out = static_cast<TOutputValue>(dist);
return out;
--
GitLab