Gaussian filter, or Gaussian blur source code

//   gaussianblur.cpp - impelementation of template class
//   of 1D and 2D gaussian blur or filter
//
//   The code is property of LIBROW
//   You can use it on your own
//   When utilizing credit LIBROW site

#ifndef _GAUSSIANBLUR_CPP_
#define _GAUSSIANBLUR_CPP_

//   Include declaration file
#include "gaussianblur.h"
//   Include math library
#include <math.h>

//   1D GAUSSIAN BLUR
//     pSignal - input signal;
//     pResult - output signal, NULL for inplace processing
//     N       - length of the signal
//     W       - window size, odd positive number
template <class T> bool TGaussianBlur<T>::Filter(T *pSignal, T *pResult,
   unsigned int N, unsigned int W) const
{
   //   Check input data cosnsitency
   if (!Consistent(pSignal, N, W))
      return false;
   //   Allocate extension
   CExtension Extension;
   if (!Extension.Allocate(N, W))
      return false;
   //   Create signal descriptor
   const CArray Signal(pSignal, N);
   //   Paste signal into extension
   Extension.Paste(Signal.Buffer);
   //   Extend signal
   Extension.Extend();
   //   Create filter window
   CWindow Window;
   if (!Window.Create(W))
      return false;
   //   Extension iterator
   const T *ExtIter = Extension.Buffer;
   //   Extension stop position
   const T *const ExtStop = Extension.Buffer + Extension.Size.x;
   //   Result iterator
   T *ResIter = pResult ? pResult : pSignal;
   //   Filter - apply to every element
   while (ExtIter < ExtStop)
      *(ResIter++) = Window.Apply(ExtIter++);
   //   Succeeded
   return true;
}

//   2D GAUSSIAN BLUR
//     pImage  - input image
//     pResult - output image, NULL for inplace processing
//     N       - width of the image
//     M       - height of the image
//     W       - window size
template <class T> bool TGaussianBlur<T>::Filter(T *pImage, T *pResult,
   unsigned int N, unsigned int M, unsigned int W) const
{
   //   Check input data consistency
   if (!Consistent(pImage, CSize(N, M), W))
     return false;
   //   Allocate extension
   CExtension Extension;
   if (!Extension.Allocate(CSize(N, M), W))
      return false;
   //   Create image descriptor
   CArray Image(pImage, CSize(N, M));
   //   Create filter window
   CWindow Window;
   if (!Window.Create(W))
      return false;
   //   Stop postion
   const T * ExtStop = Extension.Buffer + Extension.Size.x;
   //   Result iterator
   T *ResIter = pResult ? pResult : pImage;
   //   Image iterator
   const T *ImIter = Image.Buffer;
   //   Image stop position
   const T * ImStop = Image.Buffer + Image.Size.Area();
   //   Filter line by line
   while (ImIter < ImStop)
   {
      //   Paste image line into extension
      Extension.Paste(ImIter);
      //   Extend image line
      Extension.Extend();
      //   Extension iterator
      const T *ExtIter = Extension.Buffer;
      //   Apply filter to every pixel of the line
      while (ExtIter < ExtStop)
         *(ResIter++) = Window.Apply(ExtIter++);
         //   Move to the next line
      ImIter += Image.Size.x;
   }
   //   Initialize image descriptor with filter result
   Image.Buffer = pResult ? pResult : pImage;
   //   Set vertical extension mode
   Extension.SetMode(CExtension::ModeVertical);
   //   Extension stop position
   ExtStop = Extension.Buffer + Extension.Size.y;
   //   Result column iterator
   T *ResColumnIter = pResult ? pResult : pImage;
   //   Image iterator
   ImIter = Image.Buffer;
   //   Image stop position
   ImStop = Image.Buffer + Image.Size.x;
   //   Filter column by column
   while (ImIter < ImStop)
   {
      //   Paste image column into extension
      Extension.Paste(ImIter++);
      //   Extend image column
      Extension.Extend();
      //   Extension iterator
      const T *ExtIter = Extension.Buffer;
      //   Result pixel iterator
      ResIter = ResColumnIter;
      //   Apply filter to every pixel of the column
      while (ExtIter < ExtStop)
      {
         *ResIter = Window.Apply(ExtIter++);
         ResIter += Image.Size.x;
      }
      //   Move to the next column
      ++ResColumnIter;
   }
   //   Succeeded
   return true;
}

//   EXTENSION MEMORY ALLOCATION
//     _Size   - signal/image size
//     _Margin - extension margins
template <class T> bool TGaussianBlur<T>::CExtension::_Allocate(
   const CSize &_Size, unsigned int _Margin)
{
   //   Allocate extension buffer
   Buffer = new T[(_Size.x > _Size.y ? _Size.x : _Size.y) + (_Margin << 1)];
   //   Check buffer allocation
   if (!Buffer)
      //   Buffer allocation failed
      return false;
   //   Initialize size descriptors
   Size = _Size;
   Margin = _Margin;
   //   Succeeded
   return true;
}

//   PASTING DATA LINE/COLUMN INTO EXTENSION FROM DATA ARRAY
//     _Start - start postion in image/signal to paste from
template <class T> void TGaussianBlur<T>::CExtension::Paste(const T *const _Start)
{
   if (Mode == ModeHorizontal)
   {
      //   Paste line
      memcpy(Buffer + Margin, _Start, Size.x * sizeof(T));
   }
   else
   {
      //   Stop position
      const T *const Stop = _Start + Size.Area();
      //   Array iterator
      const T *ArrIter = _Start;
      //   Extension iterator
      T *ExtIter = Buffer + Margin;
      //   Paste array column element by element
      while (ArrIter < Stop)
      {
         //   Copy line
         *(ExtIter++) = *ArrIter;
         //   Jump to the next line
         ArrIter += Size.x;
      }
   }
}

//   EXTENSION CREATION
template <class T> void TGaussianBlur<T>::CExtension::Extend()
{
   //   Line size
   const unsigned int Line = Mode == ModeHorizontal ? Size.x : Size.y;
   //   Stop position
   const T *const Stop = Buffer - 1;
   //   Left extension iterator
   T *ExtLeft = Buffer + Margin - 1;
   //   Left array iterator
   const T *ArrLeft = ExtLeft + 2;
   //   Right extension iterator
   T *ExtRight = ExtLeft + Line + 1;
   //   Left array iterator
   const T *ArrRight = ExtRight - 2;
   //   Create extension line element by element
   while (ExtLeft > Stop)
   {
      //   Left extension
      *(ExtLeft--) = *(ArrLeft++);
      //   Right extension
      *(ExtRight++) = *(ArrRight--);
   }
}

//   FILTER WINDOW CREATION
//     _Size - window size
template <class T> bool TGaussianBlur<T>::CWindow::Create(unsigned int _Size)
{
   //   Allocate window buffer
   Weights = new double[_Size];
   //   Check allocation
   if (!Weights)
      //   Window buffer allocation failed
      return false;
   //   Set window size
   Size = _Size;
   //   Window half
   const unsigned int Half = Size >> 1;
   //   Central weight
   Weights[Half] = 1.;
   //   The rest of weights
   for (unsigned int Weight = 1; Weight < Half + 1; ++Weight)
   {
      //   Support point
      const double x = 3.* double(Weight) / double(Half);
      //   Corresponding symmetric weights
      Weights[Half - Weight] = Weights[Half + Weight] = exp(-x * x / 2.);
   }
   //   Weight sum
   double k = 0.;
   for (unsigned int Weight = 0; Weight < Size; ++Weight)
      k += Weights[Weight];
   //   Weight scaling
   for (unsigned int Weight = 0; Weight < Size; ++Weight)
      Weights[Weight] /= k;
   //   Succeeded
   return true;
}

#endif