GradientGaussAsymmetric.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 "GradientGaussAsymmetric.hh"

using namespace BIAS;
using namespace std;

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

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

template <class InputStorageType, class OutputStorageType>
GradientGaussAsymmetric<InputStorageType, OutputStorageType>::
GradientGaussAsymmetric(const GradientGaussAsymmetric<InputStorageType, OutputStorageType>& other)
{}

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

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

  if ((fabs(this->_GradGaussSigma-this->_LastSigma)>=DBL_EPSILON) || 
      (fabs(this->_GradGaussRatio-this->_LastRatio)>=DBL_EPSILON)){
    _CalculateKernels(this->_GradGaussSigma, this->_GradGaussRatio);
  }
  if (!src.SamePixelAndChannelCount(gx)){
    if (!gx.IsEmpty()) gx.Release();
    gx.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }
  if (!src.SamePixelAndChannelCount(gy)){
    if (!gy.IsEmpty()) gy.Release();
    gy.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }

 if (this->_Conv.Filter(src, gx)!=0) {
    BIASERR("Error in gradient x computation...");
    return -1;
  }
  FilterMask f;
  this->_Conv.GetKernel(f);
  f.Transpose();
  this->_Conv.SetKernel(f);
  int res = this->_Conv.Filter(src, gy);
  f.Transpose();
  this->_Conv.SetKernel(f);  
  return res;
}

template <class InputStorageType, class OutputStorageType>
int GradientGaussAsymmetric<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 GradientGaussAsymmetric<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, abs)\n");
  if ((fabs(this->_GradGaussSigma-this->_LastSigma)>=DBL_EPSILON) || 
      (fabs(this->_GradGaussRatio-this->_LastRatio)>=DBL_EPSILON)){
    _CalculateKernels(this->_GradGaussSigma, this->_GradGaussRatio);
  }
  if (!src.SamePixelAndChannelCount(absg)){
    if (!absg.IsEmpty()) absg.Release();
    absg.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }
  
  int res, absres=0;
  res = Filter(src, gx, gy);
  unsigned htlx, htly, hbrx, hbry, vtlx, vtly, vbrx, vbry;
  gx.GetROI()->GetCorners(htlx, htly, hbrx, hbry);
  gy.GetROI()->GetCorners(vtlx, vtly, vbrx, vbry);
  const int minx=(htlx<vtlx)?vtlx:htlx;
  const int miny=(htly<vtly)?vtly:htly;
  const int maxx=(hbrx>vbrx)?vbrx:hbrx;
  const int maxy=(hbry>vbry)?vbry:hbry;
  
  register int x, y;
  register OutputStorageType **gxi=gx.GetImageDataArray();
  register OutputStorageType **gyi=gy.GetImageDataArray();
  register OutputStorageType **absi=absg.GetImageDataArray();

  int mgx, mgy;
  for (y=miny; y<maxy; y++){
    for (x=minx; x<maxx; x++){
      mgx=(int)gxi[y][x];
      mgy=(int)gyi[y][x];
      absi[y][x]=(OutputStorageType) sqrt((double)(mgx)*(double)(mgx)+
                                          (double)(mgy)*(double)(mgy));
    }
  }
  
  absg.GetROI()->SetCorners(minx, miny, maxx, maxy);
  
  return (res==0 && absres==0)?0:-1;
  
}


template <class InputStorageType, class OutputStorageType>
void GradientGaussAsymmetric<InputStorageType, OutputStorageType>::
_CalculateKernels(double Sigma, double Ratio)
{
  /// \todo check if this is a proper filter mask in the
  /// sense of convolution -> reflection at center. We
  /// might have introduced an inconsistency here when switching
  /// to generic convolution.
  const double inv_max_gauss_deriv=1.0/(Sigma*expf(-0.5f));
  BIASCDOUT(D_GRADGAUSS_KERNEL, "max_gauss_deriv = "
            <<1.0/inv_max_gauss_deriv<<endl);
  const double fac=1.0/(2.0*Sigma*Sigma);
  const int vhws=(int)(floor(sqrt(-2.0*Sigma*Sigma*log(Ratio))));
  int k;
  for (k=vhws; k*expf((float)(-k*k*fac))*(float)inv_max_gauss_deriv>=Ratio; k++);
  const int hhws=k-1;
  const int vsize=(vhws<<1)+1;
  const int hsize=(hhws<<1)+1;

  Vector<FM_FLOAT> H(hsize), V(vsize);
  
  BIASCDOUT(D_GRADGAUSS_KERNEL, "Sigma = "<<Sigma<<"\tRatio = "<<Ratio<<endl);
  BIASCDOUT(D_GRADGAUSS_KERNEL, "vsize = "<<vsize<<"\tvhws = "<<vhws<<endl);
  BIASCDOUT(D_GRADGAUSS_KERNEL, "hsize = "<<hsize<<"\thhws = "<<hhws<<endl);

  int i, ip, im;
  register double g;
  double sum=0.0, dsum=0.0;
  // calculate the filter coefficients as gauss function
  for (i=0; i<=vhws; i++){
    ip=vhws+i;
    im=vhws-i;
    g=expf((float)(-(double)(i*i)*fac));
    V[ip] = V[im] = (FM_FLOAT)g;
    //BIASCDOUT(D_GRADGAUSS_KERNEL, i<<"v:"<<_VVec[im]<<"  "<<endl);
    sum+=g*2.0;
    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:"<<_HVec[im]<<"  "<<endl);
  }
  for (i=vhws+1; i<=hhws; i++){
    ip=hhws+i;
    im=hhws-i;
    g=double(i)*expf((float)(-(double)(i*i)*fac));
    H[im]=(FM_FLOAT)g;
    H[ip]=(FM_FLOAT)-g;
    dsum+=fabs(2.0*double(i)*g);
    //BIASCDOUT(D_GRADGAUSS_KERNEL, i<<"h:"<<_HVec[im]<<"  "<<endl);
  }

  sum-=1.0;
  //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<=vhws; i++){
    ip=vhws+i;
    im=vhws-i;
    H[ip]*=(float)dsum;
    H[im]*=(float)dsum;
    V[ip] = V[im] = V[im]*(float)sum;
  }
  for (i=vhws+1; i<=hhws; i++){
    ip=hhws+i;
    im=hhws-i;
    H[ip] *= (float)dsum;
    H[im] *= (float)dsum;
  }
  this->_LastSigma=Sigma;
  this->_LastRatio=Ratio;
  FilterMask F(H,V);
  // compute the int approximation
  int thebits = F.ComputeIntPrecisionBits(sizeof(InputStorageType)*8,
                                            sizeof(CONV_INT)*8-1);
  F.CreateIntFilter(thebits,thebits,thebits);
  this->_Conv.SetKernel(F);

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

  //#define COMPARE
#ifdef COMPARE
  fprintf(stdout, "gauss: \n");
  for (i=-hhws; i<=hhws; i++){
    fprintf(stdout, "%f ",H[i+hhws]); 
  }
  fprintf(stdout, "\n");
  fprintf(stdout, "gaussderiv: \n");
  for (i=-vhws; i<=vhws; i++){
    fprintf(stdout, "%f ",V[i+vhws]); 
  }
  fprintf(stdout, "\n");
#endif

}

//////////////////////////////////////////////////////////////////////////
//                 instantiation
//////////////////////////////////////////////////////////////////////////
namespace BIAS{
#define FILTER_INSTANTIATION_CLASS GradientGaussAsymmetric
#define FILTER_INSTANTIATION_NO_UNSIGNED_OUTPUT
#include "Filterinst.hh"
}