HessianGauss.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 "HessianGauss.hh"
#include "Convolution.hh"

using namespace BIAS;
using namespace std;

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

template <class StorageType, class OutputStorageType>
HessianGauss<StorageType, OutputStorageType>::
HessianGauss() : _ConvXX(), _ConvYY(), _ConvXY() 
{
  _LastSigma=-DBL_MAX;
  _HessGaussSigma=1.0;
  _HessGaussRatio=0.01;

}

template <class StorageType, class OutputStorageType>
HessianGauss<StorageType, OutputStorageType>::
HessianGauss(const HessianGauss<StorageType, OutputStorageType>& other)
  : _HessGaussSigma(other._HessGaussSigma),
    _HessGaussRatio(other._HessGaussRatio),
    _LastSigma(other._LastSigma),
    _LastRatio(other._LastRatio),
    _ConvXX(other._ConvXX),
    _ConvYY(other._ConvYY),
    _ConvXY(other._ConvXY)
{
}

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


template <class InputStorageType, class OutputStorageType>
int HessianGauss<InputStorageType,OutputStorageType>::
Filter(const Image<InputStorageType>& src, Image<OutputStorageType>& hxx, 
       Image<OutputStorageType>& hyy, Image<OutputStorageType>& hxy) 
{
 //  _ConvXX.SetDebugLevel(D_FILTERBASE_CALLSTACK|D_CONV_TYPES|D_WRITE_IMAGES);
//   _ConvYY.SetDebugLevel(D_FILTERBASE_CALLSTACK|D_CONV_TYPES|D_WRITE_IMAGES);
//   _ConvXY.SetDebugLevel(D_FILTERBASE_CALLSTACK|D_CONV_TYPES|D_WRITE_IMAGES);
  BIASCDOUT(D_FILTERBASE_CALLSTACK, "HessianGauss::Filter(src, hxx, gyy, hxy)\n");

  if ((fabs(_HessGaussSigma-_LastSigma)>=DBL_EPSILON) || 
      (fabs(_HessGaussRatio-_LastRatio)>=DBL_EPSILON)){
    _CalculateKernels(_HessGaussSigma, _HessGaussRatio);
  }
  if (!src.SamePixelAndChannelCount(hxy)){
    if (!hxy.IsEmpty()) hxy.Release();
    hxy.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }
   
  if (!hxx.SamePixelAndChannelCount(src)){
    hxx.Release();
    hxx.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }
  if (!hyy.SamePixelAndChannelCount(src)){
    hyy.Release();
    hyy.Init(src.GetWidth(), src.GetHeight(), src.GetChannelCount());
  }
  
  int res =0;
  if ((res =_ConvXX.Filter(src, hxx))!=0) {
    BIASERR("Error in hessian xx computation...");
    return -1;
  }
  if ((res = _ConvYY.Filter(src, hyy))!=0) {
    BIASERR("Error in hessian yy computation...");
    return -1;
  }
  //now convolute with xy kernel for 3rd image
  if ((res = _ConvXY.Filter(src, hxy))!=0) {
    BIASERR("Error in hessian xy computation...");
    return -1;
  }

  return (res==0)?0:-1;
}

template <class InputStorageType, class OutputStorageType>
void HessianGauss<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;
  Matrix<FM_FLOAT> Mxx(size, size);
  Matrix<FM_FLOAT> Myy(size, size);
  Matrix<FM_FLOAT> Mxy(size, size);

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

  // calculate the filter coefficients as second derivative of gauss function
  // this is a filtermask in the sense of convolution 
  double fac=-1.0/(2.0*Sigma*Sigma);
  //normalization factor
  register double N = 1.0/(Sigma*sqrt(2*M_PI));
  register double gxx=0,gxy=0,gyy=0;
  register double g1=0,gxx2=0,gyy2=0;
  register double eTerm=0;
  //int xp,xm,yp,ym;
  // first compute hessian
  for (int x=-hws; x<=hws; x++){
    for (int y=-hws; y<=hws; y++ ){    
      //xp= hws+x;
      //xm= hws-x;
      //yp= hws+y;
      //ym= hws-y;
      eTerm = expf((float)((double)((x*x)+(y*y))*fac));
      g1=-N/(Sigma*Sigma)*eTerm;
      gxx2=(N*x*x)/(Sigma*Sigma*Sigma*Sigma)*eTerm;
      gyy2=(N*y*y)/(Sigma*Sigma*Sigma*Sigma)*eTerm;
      gxx=g1+gxx2;
      gyy=g1+gyy2;
      gxy=(N*x*y)/(Sigma*Sigma*Sigma*Sigma)*eTerm;
      Mxx[x+hws][y+hws]= (FM_FLOAT)gxx; 
      Myy[x+hws][y+hws]= (FM_FLOAT)gyy; 
      Mxy[x+hws][y+hws]= (FM_FLOAT)gxy; 
  //     cout<<"Gxy: "<<gxy<<endl;
//       cout<<"Pixels: ("<<xp<<","<<yp<<")"<<endl;
    }
  }
  //set the filter masks
  FilterMask F(Mxx);
  BIASCDOUT(D_HESSGAUSS_KERNEL,"Filter xx:"<<Mxx);
  int thebits = F.ComputeIntPrecisionBits(sizeof(InputStorageType)*8,
                                          sizeof(CONV_INT)*8-1);
  F.CreateIntFilter(thebits,thebits,thebits);
  //  F.CreateFloatFilter();
  _ConvXX.SetKernel(F);
  FilterMask F2(Myy);
  BIASCDOUT(D_HESSGAUSS_KERNEL,"Filter yy:"<<Myy);

  thebits = F2.ComputeIntPrecisionBits(sizeof(InputStorageType)*8,
                                      sizeof(CONV_INT)*8-1);
  F2.CreateIntFilter(thebits,thebits,thebits);
  //  F2.CreateFloatFilter();
  _ConvYY.SetKernel(F2);

  FilterMask F3(Mxy);
  BIASCDOUT(D_HESSGAUSS_KERNEL,"Filter xy:"<<Mxy); 
  thebits = F3.ComputeIntPrecisionBits(sizeof(InputStorageType)*8,
                                      sizeof(CONV_INT)*8-1);
  F3.CreateIntFilter(thebits,thebits,thebits);
  //  F3.CreateFloatFilter();
  _ConvXY.SetKernel(F3);

  _LastSigma = Sigma;
  _LastRatio = Ratio;  
}

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

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

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