From 80175522b8008a3e802234b17d2c1f7c64ab17f4 Mon Sep 17 00:00:00 2001
From: Christophe Palmann <christophe.palmann@c-s.fr>
Date: Tue, 30 Jun 2015 14:11:11 +0200
Subject: [PATCH] ENH: Frost filter with correct formula (SAR)
---
.../ImageNoise/include/otbFrostImageFilter.h | 21 ++---
.../include/otbFrostImageFilter.txx | 82 +++++++++++--------
2 files changed, 54 insertions(+), 49 deletions(-)
diff --git a/Modules/Filtering/ImageNoise/include/otbFrostImageFilter.h b/Modules/Filtering/ImageNoise/include/otbFrostImageFilter.h
index 487224a17e..15de095c98 100644
--- a/Modules/Filtering/ImageNoise/include/otbFrostImageFilter.h
+++ b/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
diff --git a/Modules/Filtering/ImageNoise/include/otbFrostImageFilter.txx b/Modules/Filtering/ImageNoise/include/otbFrostImageFilter.txx
index 17197ea98a..bca00f9c7d 100644
--- a/Modules/Filtering/ImageNoise/include/otbFrostImageFilter.txx
+++ b/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));
--
GitLab