ENH: modified moments normalization

Code/FeatureExtraction/otbComplexMomentImageFunction.txx

@@ -61,19 +61,13 @@ ComplexMomentImageFunction<TInputImage, TCoordRep>
   ComplexType moments;
   moments.resize(m_Pmax+1);
   
-  std::vector<ScalarComplexType> valXpY, valXqY;
-  valXpY.resize(m_Pmax+1);
-  valXqY.resize(m_Qmax+1);
-    
   // Initialize moments
   for (unsigned int p = 0; p <= m_Pmax; p++)
     {
     moments.at(p).resize(m_Qmax+1);
-    valXpY.at(p) = ScalarComplexType(1.0,0.0);
     for (unsigned int q = 0; q <= m_Qmax; q++)
       {
       moments.at(p).at(q) =  ScalarComplexType(0.0,0.0);
-      valXqY.at(q)        =  ScalarComplexType(1.0,0.0);
       }
     }
 
@@ -105,33 +99,18 @@ ComplexMomentImageFunction<TInputImage, TCoordRep>
     {
     // Retrieve value, and centered-reduced position
     ScalarRealType value = static_cast<ScalarRealType>(it.GetPixel(i));
-    ScalarRealType     x = static_cast<ScalarRealType>(it.GetOffset(i)[0]);
-    ScalarRealType     y = static_cast<ScalarRealType>(it.GetOffset(i)[1]);
+    ScalarRealType     x = static_cast<ScalarRealType>(it.GetOffset(i)[0])/(2*m_NeighborhoodRadius+1);
+    ScalarRealType     y = static_cast<ScalarRealType>(it.GetOffset(i)[1])/(2*m_NeighborhoodRadius+1);
     
     // Build complex value
     ScalarComplexType xpy(x,y),xqy(x,-y);
     
-    unsigned int pTmp = 1;
-    unsigned int qTmp = 1;
-    
-    while (pTmp <= m_Pmax)
-      {
-      valXpY.at(pTmp) = valXpY.at(pTmp-1) * xpy;
-      pTmp ++;
-      }
-    while (qTmp <= m_Qmax)
-      {
-      valXqY.at(qTmp) = valXqY.at(qTmp-1) * xqy;
-      qTmp ++;
-      }
-    
-
     // Update cumulants
     for (unsigned int p = 0; p <= m_Pmax; p++)
       {
       for (unsigned int q= 0; q <= m_Qmax; q++)
         {
-        moments.at(p).at(q) += valXpY.at(p) * valXqY.at(q) * value;   
+        moments.at(p).at(q) += vcl_pow(xpy,p) * vcl_pow(xqy,q) * value;   
         }
       }
     }
@@ -141,7 +120,10 @@ ComplexMomentImageFunction<TInputImage, TCoordRep>
     {
     for (unsigned int q= 0; q <= m_Qmax; q++)
       {
-      moments.at(p).at(q) /= vcl_pow(moments.at(0).at(0), (p+q)/2);   
+      if (p+q != 0)
+        {
+        moments.at(p).at(q) /= moments.at(0).at(0);
+        }
       }
     }
 

Code/FeatureExtraction/otbFlusserImageFunction.txx

@@ -101,8 +101,8 @@ FlusserImageFunction<TInputImage,TCoordRep>
     {
     // Retrieve value, and centered-reduced position
     ScalarRealType value = static_cast<ScalarRealType>(it.GetPixel(i));
-    ScalarRealType x = static_cast<ScalarRealType>(it.GetOffset(i)[0]);
-    ScalarRealType y = static_cast<ScalarRealType>(it.GetOffset(i)[1]);
+    ScalarRealType x = static_cast<ScalarRealType>(it.GetOffset(i)[0])/(2*m_NeighborhoodRadius+1);
+    ScalarRealType y = static_cast<ScalarRealType>(it.GetOffset(i)[1])/(2*m_NeighborhoodRadius+1);
     
     // Build complex value
     ComplexType xpy(x,y),xqy(x,-y);
@@ -120,15 +120,15 @@ FlusserImageFunction<TInputImage,TCoordRep>
     }
   
   // Nomalisation
-  c11 /= vcl_pow(c00, (1+1)/2);
-  c12 /= vcl_pow(c00, (1+2)/2);
-  c21 /= vcl_pow(c00, (2+1)/2);
-  c20 /= vcl_pow(c00, (2+0)/2);
-  c30 /= vcl_pow(c00, (3+0)/2);
-  c22 /= vcl_pow(c00, (2+2)/2);
-  c31 /= vcl_pow(c00, (3+1)/2);
-  c40 /= vcl_pow(c00, (4+0)/2);
-  
+  c11 /= c00;
+  c12 /= c00;
+  c21 /= c00;
+  c20 /= c00;
+  c30 /= c00;
+  c22 /= c00;
+  c31 /= c00;
+  c40 /= c00;
+
   // Compute moments combinations
   moments[0]  = static_cast<ScalarRealType>(c11.real());
   moments[1]  = static_cast<ScalarRealType>((c21*c12).real());

Code/FeatureExtraction/otbHuImageFunction.txx

@@ -100,9 +100,8 @@ HuImageFunction<TInputImage, TCoordRep>
     {
     // Retrieve value, and centered-reduced position
     ScalarRealType value = static_cast<ScalarRealType>(it.GetPixel(i));
-    ScalarRealType x = static_cast<ScalarRealType>(it.GetOffset(i)[0]);
-    ScalarRealType y = static_cast<ScalarRealType>(it.GetOffset(i)[1]);
-    
+    ScalarRealType x = static_cast<ScalarRealType>(it.GetOffset(i)[0])/(2*m_NeighborhoodRadius+1);
+    ScalarRealType y = static_cast<ScalarRealType>(it.GetOffset(i)[1])/(2*m_NeighborhoodRadius+1);
     // Build complex value
     ComplexType xpy(x,y),xqy(x,-y);
     
@@ -118,13 +117,13 @@ HuImageFunction<TInputImage, TCoordRep>
     }
 
   // Nomalisation
-  c11 /= vcl_pow(c00, (1+1)/2);
-  c20 /= vcl_pow(c00, (2+0)/2);
-  c02 /= vcl_pow(c00, (0+2)/2);
-  c30 /= vcl_pow(c00, (3+0)/2);
-  c03 /= vcl_pow(c00, (0+3)/2);
-  c21 /= vcl_pow(c00, (2+1)/2);
-  c12 /= vcl_pow(c00, (1+2)/2);
+  c11 /= c00;
+  c20 /= c00;
+  c02 /= c00;
+  c30 /= c00;
+  c03 /= c00;
+  c21 /= c00;
+  c12 /= c00;
  
   // Compute moments combinations
   moments[0]  = static_cast<ScalarRealType>(c11.real());