ENH: otbReciprocalHAlphaImageFilter.h (SAR)

Modules/Filtering/Polarimetry/include/otbReciprocalHAlphaImageFilter.h

@@ -95,6 +95,7 @@ public:
     // Entropy estimation
     double totalEigenValues(0.0);
     double p[3];
+    double plog[3];
     double entropy;
     double alpha;
     double anisotropy;
@@ -120,46 +121,32 @@ public:
           }
       }
 
-    totalEigenValues = sortedRealEigenValues[0] + sortedRealEigenValues[1] + sortedRealEigenValues[2];
-    if (totalEigenValues <m_Epsilon)
+    totalEigenValues = 0.0;
+    for (unsigned int k = 0; k < 3; ++k)
       {
-        totalEigenValues = m_Epsilon;
+        sortedRealEigenValues[k] = std::max(sortedRealEigenValues[k], 0.);
+        totalEigenValues += sortedRealEigenValues[k];
       }
 
+      
     for (unsigned int k = 0; k < 3; ++k)
       {
-        p[k] = std::max(sortedRealEigenValues[k], 0.) / totalEigenValues;
+        p[k] = sortedRealEigenValues[k] / totalEigenValues;
+        
+        if (p[k]<m_Epsilon) //n=log(n)-->0 when n-->0
+			plog[k]=0.0;
+		else
+			plog[k]=-p[k]*log(p[k])/log(3);
       }
 
-    if ( (p[0] < m_Epsilon) || (p[1] < m_Epsilon) || (p[2] < m_Epsilon) )
-      {
-        entropy =0.0;
-      }
-    else
-      {
-        entropy = p[0]*log(p[0]) + p[1]*log(p[1]) + p[2]*log(p[2]);
-        entropy /= -log(3.);
-      }
+	entropy = 0.0;
+	for (unsigned int k = 0; k < 3; ++k)
+			entropy += plog[k];
 
     // alpha estimation
     double val0, val1, val2;
     double a0, a1, a2;
 
-    for(unsigned int k = 0; k < 3; ++k)
-      {
-         p[k] = sortedRealEigenValues[k] / totalEigenValues;
-
-         if (p[k] < 0.)
-           {
-             p[k] = 0.;
-           }
-
-         if (p[k] > 1.)
-           {
-             p[k] = 1.;
-           }
-      }
-
     val0 = std::abs(sortedGreaterEigenVector[0]);
     a0=acos(vcl_abs(val0)) * CONST_180_PI;
 
@@ -171,14 +158,10 @@ public:
 
     alpha=p[0]*a0 + p[1]*a1 + p[2]*a2;
 
-    if (alpha>90)
-      {
-        alpha=0.;
-      }
-
     // Anisotropy estimation
     anisotropy=(sortedRealEigenValues[1] - sortedRealEigenValues[2])/(sortedRealEigenValues[1] + sortedRealEigenValues[2] + m_Epsilon);
 
+
     result[0] = static_cast<OutputValueType>(entropy);
     result[1] = static_cast<OutputValueType>(alpha);
     result[2] = static_cast<OutputValueType>(anisotropy);