ENH: Lee filter with correct formula (SAR)

Modules/Filtering/ImageNoise/include/otbLeeImageFilter.h

@@ -29,13 +29,13 @@ namespace otb
  * \brief Anti-speckle image filter
  *
  * This class implements Lee's filter for despeckleing of SAR
- * images. The estimated reflectivity \f$R\f$ is computed as follows:
+ * images. The estimated reflectivity R is computed as follows:
+
+R=I*W+E[I]*(1-W), where
+W=1-Cu*Cu/(Ci*Ci)
+Cu = 1/sqrt(nb of look)
+Ci = sqrt(VAR[i])/E[I]
 
-\f[R = E[I] + b(I-E[I]) \f] with
-\f$ b = C^2_r / ( C^2_r + C^2_v )\f$ and \f$C_v =
-\frac{1}{\sqrt{L}}\f$,  where
-\f$L\f$ the image number of looks and
-\f$C_r = \frac{\sqrt{Var(I)}}{E[I]} \f$ and  \f$Var(I) = E[I^2] - E[I]^2\f$.
  *
  *
  * \ingroup OTBImageNoise

Modules/Filtering/ImageNoise/include/otbLeeImageFilter.txx

@@ -102,7 +102,6 @@ void LeeImageFilter<TInputImage, TOutputImage>::ThreadedGenerateData(
   itk::ConstNeighborhoodIterator<InputImageType> bit;
   itk::ImageRegionIterator<OutputImageType>      it;
 
-  // Allocate output
   typename OutputImageType::Pointer     output = this->GetOutput();
   typename InputImageType::ConstPointer input  = this->GetInput();
 
@@ -120,11 +119,10 @@ void LeeImageFilter<TInputImage, TOutputImage>::ThreadedGenerateData(
   InputRealType sum;
   InputRealType sum2;
 
-  double Cr2, Cv2, E_I, I, Var_I, dPixel;
+  double Ci2, Cu2, w, E_I, I, Var_I, dPixel;
 
   //Compute the ratio using the number of looks
-  Cv2 = 1. / (vcl_sqrt(m_NbLooks));
-  Cv2 *= Cv2;
+  Cu2 = 1.0/m_NbLooks;
 
   // Process each of the boundary faces.  These are N-d regions which border
   // the edge of the buffer.
@@ -134,39 +132,60 @@ void LeeImageFilter<TInputImage, TOutputImage>::ThreadedGenerateData(
     const 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;
+      
       //Parcours du voisinage
       for (i = 0; i < neighborhoodSize; ++i)
         {
         dPixel = static_cast<double>(bit.GetPixel(i));
         sum += dPixel;
-        sum2 += dPixel * dPixel;
         }
       E_I   = sum / static_cast<double>(neighborhoodSize);
-      Var_I = sum2 / static_cast<double>(neighborhoodSize) - E_I * E_I;
+      
+      for (i = 0; i < neighborhoodSize; ++i)
+        {
+        dPixel = static_cast<double>(bit.GetPixel(i));
+        sum2 += (dPixel-E_I) * (dPixel-E_I);
+        }
+      Var_I = sum2 / static_cast<double>(neighborhoodSize -1);
+      
       I = static_cast<double>(bit.GetCenterPixel());
+      
+      Ci2    = Var_I / (E_I * E_I);
 
       const double epsilon = 0.0000000001;
       if (vcl_abs(E_I) < epsilon)
-        {
+      {
         dPixel = itk::NumericTraits<OutputPixelType>::Zero;
-        }
+      }
+      else if (vcl_abs(Var_I) < epsilon)
+      {
+		dPixel = E_I;
+      }
+      else if (Ci2 < Cu2)
+      {
+		dPixel = E_I;
+      }
       else
-        {
-        Cr2    = Var_I / (E_I * E_I);
-        dPixel = E_I + ((I - E_I) * (Cr2)) / (Cr2 + Cv2);
-
-        }
-      // get the mean value
+      {
+		w = 1 - Cu2 / Ci2;
+		dPixel = I*w + E_I*(1-w);
+      }
+      
+      // set the weighted value
       it.Set(static_cast<OutputPixelType>(dPixel));
 
       ++bit;
       ++it;
+      
       progress.CompletedPixel();
       }
     }

Modules/Filtering/ImageNoise/test/CMakeLists.txt

@@ -29,13 +29,14 @@ otb_add_test(NAME bfTvFiltreLee1CanalPoupees COMMAND otbImageNoiseTestDriver
   ${TEMP}/bfFiltreLee_05_05_04.tif
   05 05 4.0)
 
+
 otb_add_test(NAME bfTvFiltreLee COMMAND otbImageNoiseTestDriver
-  --compare-image ${EPSILON_7}  ${BASELINE}/bfFiltreLee_05_05_12.tif
-  ${TEMP}/bfFiltreLee_05_05_12.tif
+  --compare-image ${EPSILON_7}  ${BASELINE}/GomaLee.tif 
+  ${TEMP}/bfFiltreLee_02_02_12.tif
   otbLeeFilter
   ${INPUTDATA}/GomaAvant.png    #poupees.hdr
-  ${TEMP}/bfFiltreLee_05_05_12.tif
-  05 05 12.0)
+  ${TEMP}/bfFiltreLee_02_02_12.tif
+  02 02 12.0)
 
 otb_add_test(NAME bfTuFrostFilterNew COMMAND otbImageNoiseTestDriver
   otbFrostFilterNew)