GradientGauss.cpp

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

Copyright (C) 2003-2009    (see file CONTACT 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 "GradientGauss.hh"
#include "Convolution.hh"

using namespace BIAS;
using namespace std;

//////////////////////////////////////////////////////////////////////////
//                 implementation
//////////////////////////////////////////////////////////////////////////

template <class StorageType, class OutputStorageType>
GradientGauss<StorageType, OutputStorageType>::
GradientGauss() 
{
  _LastSigma=-DBL_MAX;
  _GradGaussSigma=1.0;
  _GradGaussRatio=0.01;
}

template <class StorageType, class OutputStorageType>
GradientGauss<StorageType, OutputStorageType>::
GradientGauss(const GradientGauss<StorageType, OutputStorageType>& other)
  : _GradGaussSigma(other._GradGaussSigma),
    _GradGaussRatio(other._GradGaussRatio),
    _LastSigma(other._LastSigma),
    _LastRatio(other._LastRatio),
    _Conv(other._Conv)
{
}

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

template <class InputStorageType, class OutputStorageType>
int GradientGauss<InputStorageType,OutputStorageType>::
Filter(const Image<InputStorageType>& src, Image<OutputStorageType>& grad) 
{
  BIASCDOUT(D_FILTERBASE_CALLSTACK, "GradientGauss::Filter(src, dst)\n");

  Image<OutputStorageType> gy;
  int res = Filter(src, grad, gy);
  grad.AppendChannel(gy);
  return res;
}

template <class InputStorageType, class OutputStorageType>
int GradientGauss<InputStorageType,OutputStorageType>::
Filter(const Image<InputStorageType>& src, Image<OutputStorageType>& gx, 
       Image<OutputStorageType>& gy)
{
  BIASCDOUT(D_FILTERBASE_CALLSTACK, "GradientGauss::Filter(src, gx, gy)\n");
  //  SetDebugLevel(D_GRADGAUSS_KERNEL);
  //cout <<"GradientGauss values are "<< *_GradGaussSigma<< " "<<
  //  *_GradGaussRatio<<endl;
  if ((fabs(_GradGaussSigma-_LastSigma)>=DBL_EPSILON) || 
      (fabs(_GradGaussRatio-_LastRatio)>=DBL_EPSILON)){
    _CalculateKernels(_GradGaussSigma, _GradGaussRatio);
  }
  if (!gx.SamePixelAndChannelCount(src)){
    gx.Release();
    gx.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }
  if (!gy.SamePixelAndChannelCount(src)){
    gy.Release();
    gy.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }
  
  if (_Conv.Filter(src, gx)!=0) {
    BIASERR("Error in gradient x computation...");
    return -1;
  }
  FilterMask f;
  _Conv.GetKernel(f);
  f.Transpose();
  _Conv.SetKernel(f);
  int res = _Conv.Filter(src, gy);
  f.Transpose();
  _Conv.SetKernel(f);
  return res;
}

template <class InputStorageType, class OutputStorageType>
int GradientGauss<InputStorageType,OutputStorageType>::
Filter(const Image<InputStorageType>& src, Image<OutputStorageType>& gx, 
       Image<OutputStorageType>& gy, Image<OutputStorageType>& absg) 
{
  BIASCDOUT(D_FILTERBASE_CALLSTACK, "GradientGauss::Filter(src, gx, gy, absg)\n");

  if (!src.SamePixelAndChannelCount(absg)){
    if (!absg.IsEmpty()) absg.Release();
    absg.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }
  int res = Filter(src,gx,gy);
  int absres = this->VecLen(gx, gy, absg);
  return (res==0 && absres ==0)?0:-1;
}

template <class InputStorageType, class OutputStorageType>
void GradientGauss<InputStorageType, OutputStorageType>::
_CalculateKernels(double Sigma, double Ratio)
{
  const int hws=(int)(floor(sqrt(-2.0*Sigma*Sigma*log(Ratio))));
  const int size=(hws<<1)+1;
  Vector<FM_FLOAT> H(size), V(size);

  BIASCDOUT(D_GRADGAUSS_KERNEL, "Sigma = "<<Sigma<<"\tRatio = "<<Ratio<<endl);
  BIASCDOUT(D_GRADGAUSS_KERNEL, "size = "<<size<<"\thws = "<<hws<<endl);

  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] = (FM_FLOAT) g;
    BIASCDOUT(D_GRADGAUSS_KERNEL, i<<" v:"<<V[im]<<"  "<<endl);
    sum += g*2.0;
    // why, it is now not a symmetric filter??
    g = double(i)*g;
    H[im] = (FM_FLOAT)  g;
    H[ip] = (FM_FLOAT) -g;
    dsum+=fabs(2.0*double(i)*g);
    BIASCDOUT(D_GRADGAUSS_KERNEL, i<<" h:"<< H[im]<<"  "<<endl);
  }
  // added double, sub one here for consistency
  sum -= H[hws]; 
  BIASCDOUT(D_GRADGAUSS_KERNEL, "sum = "<<sum<<"\tdsum = "<<dsum<<endl);

  // normalize filter, so that the sum is 1.0
  sum  = 1.0 / sum;
  dsum = 1.0 / dsum;
  BIASCDOUT(D_GRADGAUSS_KERNEL, "sum = "<<sum<<"\tdsum = "<<dsum<<endl);
  for (i=0; i<=hws; i++){
    ip = hws+i;
    im = hws-i;
    H[ip] *= (FM_FLOAT)dsum;
    H[im] *= (FM_FLOAT)dsum;
    V[ip] = V[im] = V[im]*(FM_FLOAT)sum;
    BIASCDOUT(D_GRADGAUSS_KERNEL, "h:"<< H[im]<<"  "<<endl);
    //BIASCDOUT(D_GRADGAUSS_KERNEL, "v:"<<_VVec[im]<<"  "<<endl);
  }
  FilterMask F(H,V);
  int thebits = F.ComputeIntPrecisionBits(sizeof(InputStorageType)*8,
                                          sizeof(CONV_INT)*8-1);
  F.CreateIntFilter(thebits,thebits,thebits);

  _Conv.SetKernel(F);
  _LastSigma = Sigma;
  _LastRatio = Ratio;

  
#ifdef BIAS_DEBUG
  sum=dsum=0.0;
  for (i=-hws; i<=hws; i++){
    sum += H[i+hws];
    BIASCDOUT(D_GRADGAUSS_KERNEL, H[i+hws]<<"  ");
  }
  BIASCDOUT(D_GRADGAUSS_KERNEL, "  sum = "<<sum<<endl);
  for (i=-hws; i<=hws; i++){
    dsum += V[i+hws];
    BIASCDOUT(D_GRADGAUSS_KERNEL, V[i+hws]<<"  ");
  }
  BIASCDOUT(D_GRADGAUSS_KERNEL, "  sum = "<<dsum<<endl);
#endif


}

template <class InputStorageType, class OutputStorageType>
void GradientGauss<InputStorageType, OutputStorageType>::
GetBordersValid_(int &border_x, int &border_y) const
{
  _Conv.GetBorders(border_x, border_y);
}

//////////////////////////////////////////////////////////////////////////
//                 instantiation
//////////////////////////////////////////////////////////////////////////

namespace BIAS{
#define FILTER_INSTANTIATION_CLASS GradientGauss
#define FILTER_INSTANTIATION_NO_UNSIGNED_OUTPUT
#include "Filterinst.hh"
}