clfGauss.cpp

/* 
 This file is part of the BIAS library (Basic ImageAlgorithmS).

 Copyright (C) 2003, 2004    (see file CONTACTS for details)
 Multimediale Systeme der Informationsverarbeitung
 Institut fuer Informatik
 Christian-Albrechts-Universitaet Kiel


 BIAS is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.

 BIAS is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.

 You should have received a copy of the GNU Lesser General Public License
 along with BIAS; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <iostream>
#include "clfGauss.hh"
#include <Base/Math/Vector.hh>
#include <Base/Common/BIASpragma.hh>
using namespace std;
using namespace BIAS;

namespace BIAS {

  template <class InputStorageType, class OutputStorageType>
  clfGauss<InputStorageType, OutputStorageType>::
  clfGauss(clfContext *ctx) : clfFilter<InputStorageType, OutputStorageType>(ctx) {
    
    GaussSigma_ = 0.7;
    GaussRatio_ = 0.01;

    LastSigma_ = 0.0;
    LastRatio_ = 0.0;
    
    CalculateKernels_(GaussSigma_, GaussRatio_);

    try {
      program_->AddKernel("ConvolveX");
      program_->AddKernel("ConvolveY");
    } catch (clfException &e) {
      cout << "no good:" << e.GetMessageString()<<endl;
    }

  }

  template <class InputStorageType, class OutputStorageType>
  clfGauss<InputStorageType, OutputStorageType>::
  ~clfGauss() {

  }
  
  template <class InputStorageType, class OutputStorageType>
  int clfGauss<InputStorageType, OutputStorageType>::
  Filter(clfImage2D *src, clfImage2D *dst) {
    return clfFilter<InputStorageType, OutputStorageType>::FilterSep_("ConvolveX", "ConvolveY", src, dst);
  }
  
  template <class InputStorageType, class OutputStorageType>
  int clfGauss<InputStorageType, OutputStorageType>::
  Build(unsigned int size) {
    SetWinSize(size, true);
    CalculateKernels_(GaussSigma_, GaussRatio_);
    return 0;
  }
  
  template <class InputStorageType, class OutputStorageType>
  void clfGauss<InputStorageType, OutputStorageType>::
  CalculateKernels_(double Sigma, double Ratio)
  {
    if (LastSigma_ == GaussSigma_ && LastRatio_ == GaussRatio_) return;
    LastSigma_ = GaussSigma_;
    LastRatio_ = GaussRatio_;
    
    const int hws = (int)(floor(sqrt(-2.0*Sigma*Sigma*log(Ratio))));
    const int size = (hws<<1)+1;
    
    BIAS::Vector<float> H(size), V(size);

    double sum=0.0;
    int i, ip, im;
    if (Sigma==0.0){
      H[0]=1.0; V[0]=1.0;
    } else { 
      double fac=1.0/(2.0*Sigma*Sigma);
      register double g;
      // calculate the filter coefficients as gauss function
      for (i=0; i<=hws; i++){
        ip = hws+i;
        im = hws-i;
        g  = expf((float)(-(double)(i*i)*fac));
        V[ip] = V[im] = H[im] = H[ip] = (float)g;
        sum += g*2.0;
      }
      sum -= V[hws]; // added twice

      // normalize filter, so that the sum is 1.0
      sum = 1.0/sum;
      for (i=0; i<=hws; i++){
        ip=hws+i;
        im=hws-i;
        H[ip] = H[im] = V[ip] = V[im] = V[im] * (float)sum;
      }
    }
//
//    sum = 0.0;
//    for (i=-hws; i<=hws; i++){
//      sum += H[i+hws];
//    }
//    sum = 0.0;
//    for (i=-hws; i<=hws; i++){
//      sum += V[i+hws];
//    }
//
    clfFilter<InputStorageType, OutputStorageType>::SetFilter_(size, H.GetData());
  }

  template class BIASOpenCLFramework_EXPORT clfGauss<float, float>;

} /* namespace DECS */