clfGradientGauss.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 <Base/Common/BIASpragma.hh>
#include <OpenCLFramework/Filter/clfGradientGauss.hh>
#include <Base/Math/Vector.hh>
#include <Filter/FilterMask.hh>

using BIAS::FilterMask;

namespace BIAS {

  template <class InputStorageType, class OutputStorageType>
  clfGradientGauss<InputStorageType, OutputStorageType>::
  clfGradientGauss(clfContext *ctx) : clfGauss<InputStorageType, OutputStorageType>(ctx) {
      GaussSigma_=1.0;
      GaussRatio_=0.01;
      filterMaskVert_ = ctx->CreateBuffer();
      try {
        program_->AddKernel("ConvolveGradGaussXY");
      } catch (clfException &e) {
//        cout << "no good" << endl;
      }
  }

  template <class InputStorageType, class OutputStorageType>
  clfGradientGauss<InputStorageType, OutputStorageType>::
  ~clfGradientGauss() {
    delete filterMaskVert_;
  }
  
  template <class InputStorageType, class OutputStorageType>
  int clfGradientGauss<InputStorageType, OutputStorageType>::
  Filter(clfImage2D *src, clfImage2D *dst) {
    CalculateKernels_(GaussSigma_, GaussRatio_);
    tmpImage_->AllocateFromTemplate( *src, false, false);
    program_->KernelSetArgument( "ConvolveGradGaussXY", 0, *src );
    program_->KernelSetArgument( "ConvolveGradGaussXY", 1, *filterMask_ );
    program_->KernelSetArgument( "ConvolveGradGaussXY", 2, *filterMaskVert_ );
    program_->KernelSetArgument( "ConvolveGradGaussXY", 3, (int)winSize_ );
    program_->KernelSetArgument( "ConvolveGradGaussXY", 4, *tmpImage_ );
    unsigned int w,h;
    dst->GetImageDim(w,h);

    unsigned int localx=1;
    unsigned int localy=1;
    if (w%8==0) localx = 8;
    if (w%16==0) localx = 16;
    if (h%8==0) localy = 8;
    if (h%16==0) localy = 16;
    context_->RunOn2DRange( *program_, "ConvolveGradGaussXY", w, h, localx, localy);

    program_->KernelSetArgument( "ConvolveGradGaussXY", 0, *tmpImage_ );
    program_->KernelSetArgument( "ConvolveGradGaussXY", 2, *filterMask_ );
    program_->KernelSetArgument( "ConvolveGradGaussXY", 1, *filterMaskVert_ );
    program_->KernelSetArgument( "ConvolveGradGaussXY", 4, *dst );

    context_->RunOn2DRange( *program_, "ConvolveGradGaussXY", w, h, localx, localy);

    return 0;
  }
  
  template <class InputStorageType, class OutputStorageType>
  int clfGradientGauss<InputStorageType, OutputStorageType>::
  SetFilterVert_(unsigned int size, const float *mask) {
    vertSize_ = size;
//    std::cout << vertSize_ << std::endl;
    try {
      filterMaskVert_->Allocate( vertSize_ * sizeof(float), true, false );
      filterMaskVert_->WriteToBuffer( (float*)mask );
    } catch (clfException &e) {
      BIASERR(e.GetDetailedString());
      return 1;
    }
    return 0;
  }
  
  template <class InputStorageType, class OutputStorageType>
  int clfGradientGauss<InputStorageType, OutputStorageType>::
  SetFilterHoriz_(unsigned int size, const float *mask) {
    return clfFilter<InputStorageType, OutputStorageType>::SetFilter1D_(size, mask);
  }
  
  template <class InputStorageType, class OutputStorageType>
  void clfGradientGauss<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);

    int i, ip, im;
    double fac=1.0/(2.0*Sigma*Sigma);
    register double g;
    double sum=0.0, dsum=0.0;
    // calculate the filter coefficients as gauss function
    // this is a filtermask in the sense of convolution
    // it is reflected at its symmetry center before being "dropped"
    // onto the image !
    for (i=0; i<=hws; i++){
      ip=hws-i;
      im=hws+i;
      g=expf((float)(-(double)(i*i)*fac));
      V[ip] = V[im] = (float) g;
      sum += g*2.0;
      // why, it is now not a symmetric filter??
      g = double(i)*g;
      H[im] = (float)  g;
      H[ip] = (float) -g;
      dsum+=fabs(2.0*double(i)*g);
    }
    // added double, sub one here for consistency
    sum -= H[hws]; 

    // normalize filter, so that the sum is 1.0
    sum  = 1.0 / sum;
    dsum = 1.0 / dsum;
    for (i=0; i<=hws; i++){
      ip = hws+i;
      im = hws-i;
      H[ip] *= (float)dsum;
      H[im] *= (float)dsum;
      V[ip] = V[im] = V[im]*(float)sum;
    }
    
//    std::cout << "horiz: " << H << std::endl << "vert: "<< V << std::endl;
//
    SetFilterHoriz_(size, H.GetData());
    SetFilterVert_(size, V.GetData());
  }
  
  template class BIASOpenCLFramework_EXPORT clfGradientGauss<float, float>;

} /* namespace DECS */