ENH: Frost filter with correct formula (SAR)

Modules/Filtering/ImageNoise/include/otbFrostImageFilter.h

@@ -29,18 +29,11 @@ namespace otb
  * \brief Anti-speckle image filter
  *
  * Uses a negative exponential convolution kernel.
- * The output of the filter for pixel p is:
- *      \f$ \hat I_{s}=\sum_{p\in\eta_{p}} m_{p}I_{p} \f$
- *
- * where :   \f$ m_{p}=\frac{KC_{s}^{2}\exp(-KC_{s}^{2}d_{s, p})}{\sum_{p\in\eta_{p}} KC_{s}^{2}\exp(-KC_{s}^{2}d_{s, p})} \f$
- *    and  \f$ d_{s, p}=\sqrt{(i-i_{p})^2+(j-j_{p})^2} \f$
- *
- * \f$ K \f$     : the decrease coefficient
- * \f$ (i, j)\f$ : the coordinates of the pixel inside the region
- * defined by \f$ \eta_{s} \f$
- * \f$ (i_{p}, j_{p})\f$ : the coordinates of the pixels belonging to \f$ \eta_{p} \subset \eta_{s} \f$
- * \f$ C_{s}\f$ : the variation coefficient computed over \f$ \eta_{p}\f$
- *
+ * The kernel is defined as follows:
+ * exp(-A*D), where 
+ * D is the distance from the current pixel to the center pixel
+ * A = k*Ci*Ci  with Ci = VAR[I]/ (E[I]*E[I])
+ * The final result is normalized by the sum of the kernel coefficients.
  *
  * \ingroup OTBImageNoise
  */
@@ -92,9 +85,9 @@ public:
   /** Get the radius used to define the neighborhood for the filter calculation. */
   itkGetConstReferenceMacro(Radius, SizeType);
 
-  /** Set The numbers of view used for the filter calculation. */
+  /** Set the damping factor. */
   itkSetMacro(Deramp, double);
-  /** Get The numbers of view used for the filter calculation. */
+  /** Get the damping factor. */
   itkGetConstReferenceMacro(Deramp, double);
 
   /** To be allowed to use the pipeline method FrostImageFilter needs

Modules/Filtering/ImageNoise/include/otbFrostImageFilter.txx

@@ -38,7 +38,7 @@ template <class TInputImage, class TOutputImage>
 FrostImageFilter<TInputImage, TOutputImage>::FrostImageFilter()
 {
   m_Radius.Fill(1);
-  m_Deramp = 0.1;
+  m_Deramp = 2;
 }
 
 template <class TInputImage, class TOutputImage>
@@ -135,54 +135,66 @@ void FrostImageFilter<TInputImage, TOutputImage>::ThreadedGenerateData(
     unsigned int neighborhoodSize = bit.Size();
     it = itk::ImageRegionIterator<OutputImageType>(output, *fit);
     bit.OverrideBoundaryCondition(&nbc);
+    
     bit.GoToBegin();
+    it.GoToBegin();
 
     while (!bit.IsAtEnd())
       {
       sum  = itk::NumericTraits<InputRealType>::Zero;
       sum2 = itk::NumericTraits<InputRealType>::Zero;
+      
       for (i = 0; i < neighborhoodSize; ++i)
         {
         dPixel = static_cast<double>(bit.GetPixel(i));
         sum += dPixel;
-        sum2 += dPixel * dPixel;
         }
-      Mean   = sum  / double(neighborhoodSize);
-      Variance  = sum2 / double(neighborhoodSize) - Mean * Mean;
-
-      if (Mean == 0)
-        {
-        Alpha = 0;
-        }
-      else
-        {
-        Alpha = m_Deramp * Variance / (Mean * Mean);
-        }
-
-      NormFilter  = 0.0;
-      FrostFilter = 0.0;
-
-      const int rad_x = m_Radius[0];
-      const int rad_y = m_Radius[1];
-
-      for (int x = -rad_x; x <= rad_x; ++x)
+      Mean   = sum / static_cast<double>(neighborhoodSize);
+      
+      for (i = 0; i < neighborhoodSize; ++i)
         {
-        for (int y = -rad_y; y <= rad_y; ++y)
-          {
-          double Dist = vcl_sqrt(static_cast<double>(x * x + y * y));
-          off[0] = x;
-          off[1] = y;
-
-          dPixel = static_cast<double>(bit.GetPixel(off));
-
-          CoefFilter = Alpha * vcl_exp(-Alpha * Dist);
-          NormFilter += CoefFilter;
-          FrostFilter += (CoefFilter * dPixel);
-          }
+        dPixel = static_cast<double>(bit.GetPixel(i));
+        sum2 += (dPixel-Mean) * (dPixel-Mean);
         }
+      Variance  = sum2 / double(neighborhoodSize-1);
 
-      if (NormFilter == 0.) dPixel = 0.;
-      else dPixel = FrostFilter / NormFilter;
+      const double epsilon = 0.0000000001;
+      if (vcl_abs(Mean) < epsilon)
+      {
+        dPixel = itk::NumericTraits<OutputPixelType>::Zero;
+      }
+      else if (vcl_abs(Variance) < epsilon)
+      {
+		dPixel = Mean;
+      }
+	  else
+	  {
+		  Alpha = m_Deramp * Variance / (Mean * Mean);
+
+		  NormFilter  = 0.0;
+		  FrostFilter = 0.0;
+
+		  const int rad_x = m_Radius[0];
+		  const int rad_y = m_Radius[1];
+
+		  for (int x = -rad_x; x <= rad_x; ++x)
+			{
+			for (int y = -rad_y; y <= rad_y; ++y)
+			  {
+			  double Dist = vcl_sqrt(static_cast<double>(x * x + y * y));
+			  off[0] = x;
+			  off[1] = y;
+
+			  dPixel = static_cast<double>(bit.GetPixel(off));
+
+			  CoefFilter = vcl_exp(-Alpha * Dist);
+			  NormFilter += CoefFilter;
+			  FrostFilter += (CoefFilter * dPixel);
+			  }
+			}
+
+		  dPixel = FrostFilter / NormFilter;
+	  }
 
       it.Set(static_cast<OutputPixelType>(dPixel));