converted mean and gauss filtering to FilterT

core/image/FilterT.h

@@ -96,6 +96,63 @@ class FilterT
     * @return void
     **/
     static void gradientStrength ( const NICE::ImageT<SrcType> &src, NICE::ImageT<DstType> &dst );
+    
+    
+    /**
+    * Blurs Image \c src into the Image \c dst using a mean filter.
+    * @param src  source gray image
+    * @param size mean mask of size (2*\c size + 1) x (2*\c size + 1)
+    * @param dst  Optional buffer to be used as target.<br>
+    *             Create a new Image if \c dst == NULL.<br>
+    *             If \c dst != NULL then size must be equal to \c src 's size!
+    * @throw ImageException will be thrown if \c dst != NULL and the size of \c src and \c dst is not equal.
+    * @return Pointer to Image
+    */
+    ImageT<DstType> * filterMean ( const ImageT<SrcType>& src, const uint& size, ImageT<DstType>* dst = NULL );     
+    
+    /**
+      * Blurs Image \c src into the Image \c dst using a mean filter. This is the
+      * implementation of the filter with a runtime independent of the filter size.
+      * It is significantly faster for large filter sizes. An analysis can be done with the program
+      * testApproxGaussFilter.
+      * The second template parameter gives the type used to store sums of gray-values. The third template
+      * parameter specifies the data type of the temporary image used to handle the seperability of the mean filter.
+      * Boundary handling is done with a cropped filter if use_filtersize_independent_implementation is set to true, otherwise the borders of the image
+      * are ignored.
+      * @param src  source gray image
+      * @param size mean mask of size (2*\c size + 1) x (2*\c size + 1)
+      * @param dst  Optional buffer to be used as target.<br>
+      *             Create a new Image if \c dst == NULL.<br>
+      *             If \c dst != NULL then size must be equal to \c src 's size!
+      * @throw ImageException will be thrown if \c dst != NULL and the size of \c src and \c dst is not equal.
+      * @return Pointer to Image
+      */
+    ImageT<DstType> * filterMeanLargeFS ( const ImageT<SrcType>& src, const uint& size, ImageT<DstType>* dst = NULL );    
+      
+    /**
+    * @brief Blurs Image \c src into the Image \c dst by approximating a gauss filter with a specific standard deviation through multiple mean filters
+    * The size of the mean filters is automatically calculated
+    * The second template parameter gives the type used to store sums of gray-values. The third template parameter,
+    * which is only necessary if use_filtersize_independent_implementation
+    * is set to true, specifies the data type of the temporary image used to handle the seperability of the mean filter.
+    * Boundary handling is done with a cropped filter if use_filtersize_independent_implementation is set to true, otherwise the borders of the image
+    * are ignored.
+      * @param src  source gray image
+      * @param sigma standard deviation of the gauss filter
+      * @param dst  Optional buffer to be used as target.<br>
+      *             Create a new Image if \c dst == NULL.<br>
+      *             If \c dst != NULL then size must be equal to \c src 's size!
+    * @param use_filtersize_independent_implementation  if this flag is set to true, we use the implementation of the mean filter with a runtime
+    *  independent of the filter size. This is beneficial for large filter sizes. However, for small filter sizes (<10) the standard implementation is faster.
+    *  An analysis can be done with the program testApproxGaussFilter.
+      * @throw ImageException will be thrown if \c dst != NULL and the size of \c src and \c dst is not equal
+      *                       or an invalid size \c size is specified.
+      * @return Pointer to Image
+      */
+    NICE::ImageT<DstType> * filterGaussSigmaApproximate ( const NICE::ImageT<SrcType> &src, double sigma, NICE::ImageT<DstType> *dst = NULL,
+        bool use_filtersize_independent_implementation = true );  
+  
+ 
 };
 
 // float specializations using IPP

core/image/FilterT.tcc

@@ -146,4 +146,175 @@ void FilterT<SrcType, CalcType, DstType>::gradientStrength ( const NICE::ImageT<
   return;
 }
 
+template<class SrcType, class CalcType, class DstType>
+ImageT<DstType> * FilterT<SrcType, CalcType, DstType>::filterMean ( const NICE::ImageT<SrcType>& src, const uint& size, ImageT<DstType>* dst )
+{
+  ImageT<DstType>* result = createResultBuffer ( src, dst );
+
+  int isize = size;
+  int is = isize + 1;
+  int msize = ( 2 * isize + 1 ) * ( 2 * isize + 1 );
+  CalcType sum = 0;
+
+  int xend = src.width() - isize;
+
+  const SrcType* pSrc;
+  DstType* pDst;
+  for ( int y = isize; y < src.height() - isize; ++y ) {
+    sum = 0;
+
+    for ( int j = y - isize; j <= y + isize; ++j ) {
+      pSrc = src.getPixelPointerY ( j );
+      for ( uint i = 0; i <= 2*size; ++i, ++pSrc )
+        sum += *pSrc;
+    }
+
+    pSrc = src.getPixelPointerXY ( is, y );
+    pDst = result->getPixelPointerXY ( isize, y );
+    *pDst = sum / msize;
+    ++pDst;
+
+    for ( int x = is; x < xend; ++x, ++pSrc, ++pDst ) {
+      for ( int i = -isize; i <= isize; ++i ) {
+        sum += * ( pSrc + size + i * src.getStepsize() );
+        sum -= * ( pSrc - is + i * src.getStepsize() );
+      }
+      *pDst = sum / msize;
+    }
+  }
+
+  return result;
+}
+
+template<class SrcType, class CalcType, class DstType>
+ImageT<DstType> * FilterT<SrcType, CalcType, DstType>::filterMeanLargeFS ( const ImageT<SrcType>& src, const uint& size, ImageT<DstType>* dst )
+{
+  ImageT<DstType>* result = createResultBuffer ( src, dst );
+  ImageT<CalcType> tmp ( src.width(), src.height() );
+
+  unsigned int isize = size;
+  unsigned int n = 2 * isize + 1;
+  CalcType sum = 0;
+
+  unsigned int xend = src.width();
+  unsigned int yend = src.height();
+
+  // legend of the following comments
+  // current pixel in dst = x
+  // end of filter mask and current pixel in src = e
+  // the examples use a size=2 filter, is=3, n=5
+  for ( int y = 0; y < src.height(); ++y ) {
+    sum = 0;
+
+    const SrcType* pSrc = src.getPixelPointerY ( y );
+    DstType *pDst = tmp.getPixelPointerY ( y );
+
+    // first position: [eoooooooooooooooooo]
+    // last position:  [-eooooooooooooooooo]
+    for ( unsigned int e = 0; e < size; ++e ) {
+      sum += * ( pSrc );
+      pSrc++;
+    }
+    // first position: [x-eoooooooooooooooo]
+    // last position:  [s-x-eoooooooooooooo]
+    for ( unsigned int e = size; e < n; ++e ) {
+      sum += * ( pSrc );
+      pSrc++;
+      *pDst = sum / ( e + 1 );
+      pDst++;
+    }
+    // first position: [os-x-eooooooooooooo]
+    // last position:  [oooooooooooooos-x-e]
+    for ( unsigned int e = n;e < xend; ++e ) {
+      sum -= * ( pSrc - n );
+      sum += * ( pSrc );
+      pSrc++;
+      *pDst = sum / n;
+      pDst++;
+    }
+    // first position: [ooooooooooooooos-x-]
+    // last position:  [oooooooooooooooos-x]
+    for ( unsigned int x = xend - size; x < xend; ++x ) {
+      sum -= * ( pSrc - ( xend - x + 1 ) );
+      *pDst = sum / ( size + xend - x );
+      pDst++;
+    }
+  }
+
+  // now let us filter along the columns
+  long long linestep = src.rowStepsize() / src.bytedepth();
+
+  for ( unsigned int x = 0; x < xend; ++x ) {
+    sum = 0;
+
+    const CalcType* pSrc = tmp.getPixelPointerXY ( x, 0 );
+    SrcType *pDst = result->getPixelPointerXY ( x, 0 );
+
+    for ( unsigned int e = 0; e < size; ++e ) {
+      sum += * ( pSrc );
+      pSrc += linestep;
+    }
+    for ( unsigned int e = size; e < n; ++e ) {
+      sum += * ( pSrc );
+      pSrc += linestep;
+      *pDst = sum / ( e + 1 );
+      pDst += linestep;
+    }
+    for ( unsigned int e = n;e < yend; ++e ) {
+      sum -= * ( pSrc - n * linestep );
+      sum += * ( pSrc );
+      pSrc += linestep;
+      *pDst = sum / n;
+      pDst += linestep;
+    }
+    for ( unsigned int y = yend - size; y < yend; ++y ) {
+      sum -= * ( pSrc - ( yend - y + 1 ) * linestep );
+      *pDst = sum / ( size + yend - y );
+      pDst += linestep;
+    }
+  }
+
+  return result;
+}    
+
+template<class SrcType, class CalcType, class DstType>
+ImageT<DstType> * FilterT<SrcType, CalcType, DstType>::filterGaussSigmaApproximate ( const NICE::ImageT<SrcType> &src, double sigma, NICE::ImageT<DstType> *dst,
+    bool use_filtersize_independent_implementation )
+{
+  // We use the idea of Wells 1986
+  // Efficient Synthesis of Gaussian Filters by Cascaded Uniform Filters
+  // http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=04767776
+  // However, the idea is mainly from the lecture script of Prof. Schukat-Talamazzini.
+  //
+  // 4x mean filter of size (2M+1) x (2M+1) is nearly the same as a
+  // gauss filter with sigma = sqrt(4/3(M^2+M))
+  // stddev of a mean filter of size 2M+1 is 1/3 M (M+1)
+  // after applying it 4 times (convolution) the stddev is 4/3 M (M+1)
+  // After some school math stuff (solving quadratic equations),
+  // we can derive the formula used below.
+  // copy-n-paste for gnuplot: plot [0:10] (sqrt(1+3*x*x)-1)/2
+  ImageT<DstType>* result = createResultBuffer ( src, dst );
+  NICE::ImageT<CalcType> tmp ( src.width(), src.height() );
+
+  int M = ( int ) ( ( sqrt ( 1 + 3 * sigma * sigma ) - 1 ) / 2.0 + 0.5 );
+  if ( M < 1 ) {
+    fthrow ( Exception, "Unable to approximate an Gaussian filter of this small scale (sigma=" << sigma << ")" );
+  }
+
+  if ( use_filtersize_independent_implementation )
+  {
+    this->filterMeanLargeFS ( src, M, &tmp );
+    this->filterMeanLargeFS ( tmp, M, result );
+    this->filterMeanLargeFS ( *result, M, &tmp );
+    this->filterMeanLargeFS ( tmp, M, result );
+  } else {
+    this->filterMean ( src, M, &tmp );
+    this->filterMean ( tmp, M, result );
+    this->filterMean ( *result, M, &tmp );
+    this->filterMean ( tmp, M, result );
+  }
+
+  return result;
+}
+
 } //namespace

core/vector/Distance.tcc

@@ -17,7 +17,7 @@ template<class T>
 T EuclidianDistance<T>::doCalculate(const VectorT<T>& v1, const VectorT<T>& v2) const {
 
     if(v1.size()!=v2.size())
-        _THROW_EVector("Input vectors must have the same size.");
+        _THROW_EVector("Input vectors must have the same size. v1: ");
 
     T dist = T(0);