PyramidImage.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 <bias_config.h>
#include "PyramidImage.hh"
#include "Base/Image/ImageIO.hh"
#include "Base/Image/ImageConvert.hh"

#include "Base/Math/Matrix2x2.hh"
#include "Base/Math/Vector2.hh"

using namespace BIAS;
using namespace std;

/// dont create levels smaller than this in auto mode:
#define SMALLEST_SIZE 32.0

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


template <class StorageType>
PyramidImage<StorageType>::
PyramidImage()
  : _Images(), _rescale()
{
  //  resize(2); // NO: THIS CREATES NULL-POINTERS!
  _RescaleFactor=2.0;
  _rescale.SetFactor(_RescaleFactor);
  _PositionOffset=0.0;

  // setting sigma here is completely useless, since sigma value will be
  // overwritten in BIAS::Rescale::UpdateLowpass(...)!!
  // furthermore, using gauss for filtering leads to very blurry images
  //Gauss<StorageType, StorageType> myGauss;
  //myGauss.SetSigma(1.0);
  //SetLowPassFilter(myGauss);
}

template <class StorageType>
PyramidImage<StorageType>::
PyramidImage(const PyramidImage<StorageType>& pim)
  : _Images(), _rescale()
{
  _RescaleFactor=2.0;
  *this=pim;
}


template <class StorageType>
PyramidImage<StorageType>::
PyramidImage(const double factor,
             const vector<SharedPtr<Image<StorageType> > > &imgs)
  : _Images(imgs), _RescaleFactor(factor), _vFactor()
{
  const int size = imgs.size();
  _vFactor.resize(size);
  _vFactor[0] = factor;
  for (int i=1; i<size; i++){
    _vFactor[i] = _vFactor[i-1] * factor;
  }
  _PositionOffset=0.0;
}


template <class StorageType>
PyramidImage<StorageType>::
~PyramidImage()
{
  clear();
}


template <class StorageType>
void PyramidImage<StorageType>::
Init(const Image<StorageType>& image, const unsigned size)
{
  if (!IsEmpty()) {
#ifdef BIAS_DEBUG
    //    BIASERR("warning: re-initializing pyramid image");
#endif
    Clear();
  }

  int thesize = size;
  if (thesize==0) {
    // auto pyramid size
    double smallerimagesize = image.GetWidth();
    if (smallerimagesize>image.GetHeight())
      smallerimagesize = image.GetHeight();
    BIASASSERT(smallerimagesize>0.5);
    BIASASSERT(_RescaleFactor>1.1);
    thesize = int(log(smallerimagesize/SMALLEST_SIZE)/log(_RescaleFactor));
    if (thesize<1) thesize = 1;
    //cout<<"smaller image size "<<smallerimagesize<<" led to auto pyramid "
    //    <<" size of "<<thesize<<" for pyramid base "<<_RescaleFactor<<endl;
  }

  SharedPtr<Image<StorageType> > L0 =
    SharedPtr<Image<StorageType> >(new Image<StorageType>(image));
  push_back(L0);
  for (int i=1; i<thesize; i++)
    push_back(SharedPtr<Image<StorageType> >(new Image<StorageType>));
  _CreateLevels();
  if (image.GetROI()->GetROIType() == ROI_Rows)
    SetROI(*image.GetROI());
}


template <class StorageType>
void PyramidImage<StorageType>::
InitFromImageBase(const ImageBase& image, const unsigned size)
{
  if (!IsEmpty())  {
#ifdef BIAS_DEBUG
    BIASERR("warning: re-initializing pyramid image");
#endif
    Clear();
  }

  int thesize = size;
  if (thesize==0) {
    // auto pyramid size
    double smallerimagesize = image.GetWidth();
    if (smallerimagesize>image.GetHeight())
      smallerimagesize = image.GetHeight();
    BIASASSERT(smallerimagesize>0.5);
    BIASASSERT(_RescaleFactor>1.1);
    thesize = int(log(smallerimagesize/SMALLEST_SIZE)/log(_RescaleFactor));
    if (thesize<1) thesize = 1;
    cout<<"smaller image size is "<<smallerimagesize<<" lead to auto pyramid "
        <<" size of "<<thesize<<" for pyramid base "<<_RescaleFactor<<endl;
  } else {
#ifdef BIAS_DEBUG
      double smallerimagesize = image.GetWidth();
      if (smallerimagesize>image.GetHeight())
        smallerimagesize = image.GetHeight();
      BIASASSERT(smallerimagesize>0.5);
      BIASASSERT(_RescaleFactor>1.1);
      int maxsize = int(log(smallerimagesize/2.0)/log(_RescaleFactor));
      if (maxsize<1) maxsize = 1;
//      cout<<"smaller image size is "<<smallerimagesize<<" lead to max pyramid"
//          <<" size of "<<thesize<<" for pyramid base "<<_RescaleFactor<<endl;
      BIASASSERT(thesize<=maxsize);
#endif
  }


  Image<StorageType> *L0 =
    new Image<StorageType>(image.GetWidth(), image.GetHeight(),
                           image.GetChannelCount());
  ImageConvert::ConvertST(image, *L0, L0->GetStorageType());
  push_back(SharedPtr<Image<StorageType> >(L0));
  for (int i=1; i<thesize; i++)
    push_back(SharedPtr<Image<StorageType> >(new Image<StorageType>));
  _CreateLevels();
}


template <class StorageType>
void PyramidImage<StorageType>::
Init(const unsigned pyramid_size)
{
  if (!IsEmpty()) {
#ifdef BIAS_DEBUG
    BIASERR("warning: re-initializing pyramid image");
#endif
    Clear();
  }
  resize(pyramid_size);
}


template <class StorageType>
void PyramidImage<StorageType>::
Init(const unsigned width, const unsigned height,
     const unsigned channelcount, const unsigned pysize)
{
  if (!IsEmpty()) {
#ifdef BIAS_DEBUG
    BIASERR("warning: re-initializing pyramid image");
#endif
    Clear();
  }
  resize(pysize);
  unsigned int newwidth=0, newheight=0;
  double pyfac = 1.0;
  const double factor=_RescaleFactor;

  for (unsigned i=0; i<pysize; i++){
    newwidth = (unsigned int)rint((double)width / pyfac);
    newheight = (unsigned int)rint((double)height / pyfac);
    if (!_Images[i]->IsEmpty()) _Images[i]->Release();
    _Images[i]->Init(newwidth, newheight, channelcount);
    _vFactor[i] = pyfac;
    pyfac*=factor;
  }
}

template <class StorageType>
void PyramidImage<StorageType>::
Clear()
{
  _Images.clear();
  _vFactor.clear();
  _PositionOffset=0.0;
}

template <class StorageType>
void PyramidImage<StorageType>::
resize(const unsigned size)
{
  const unsigned current_size = _Images.size();
  if (size>current_size){
    for (unsigned i= current_size;i<size; i++){
      _Images.push_back(SharedPtr<Image<StorageType> >(new Image<StorageType>));
      _vFactor.push_back(0.0);
    }
  }
}

template <class StorageType>
PyramidImage<StorageType>& PyramidImage<StorageType>::
ShallowCopy(const PyramidImage<StorageType>& pim)
{
  _rescale=pim._rescale;
  _RescaleFactor=pim._RescaleFactor;
  _vFactor = pim.GetFactors();
  _PositionOffset=pim.GetPositionOffset();

  if (!IsEmpty()) Clear();
  resize(pim.size());
  unsigned size = pim.size();
  for (unsigned i=0; i<size; i++){
    (*this)[i] = ConstCast<Image<StorageType> >(pim[i]);
  }
  return *this;
}


template <class StorageType>
PyramidImageInterface<StorageType> *PyramidImage<StorageType>::
ShallowClone() const
{
  PyramidImage<StorageType> *res = new PyramidImage<StorageType>;
  res->ShallowCopy(*this);
  return res;
}


template <class StorageType>
PyramidImage<StorageType>& PyramidImage<StorageType>::
operator=(const PyramidImage<StorageType>& pim)
{
  _rescale=pim._rescale;
  _RescaleFactor=pim._RescaleFactor;
  _vFactor = pim.GetFactors();

  if (!IsEmpty()) Clear();
  resize(pim.size());
  typename std::vector<SharedPtr<Image<StorageType> > >::iterator it;
  typename std::vector<SharedPtr<Image<StorageType> > >::const_iterator cit;

  for (it=_Images.begin(), cit=pim._Images.begin();
       it!=_Images.end() && cit!=pim._Images.end(); it++, cit++){
    **it = **cit;
  }
  _PositionOffset=pim.GetPositionOffset();

  return *this;
}

/////////////////////////////////////////////////////////////////////////////

template <class StorageType>
int PyramidImage<StorageType>::
Downsample()
{
  int res=0;
  typename vector<SharedPtr<Image<StorageType> > >::iterator it, lastit;
  for (lastit=_Images.begin(), it=lastit+1;
       it!=_Images.end() && lastit!=_Images.end(); it++, lastit++){
    res += _rescale.Downsample(*(*lastit), *(*it));
    (*it)->SetOutsideROIZero();
//    cout << "setting pixels outside roi to zero--------------- " << (*it)->GetROI() << endl;

  }

  _PositionOffset= _rescale.GetPositionOffset();
  return res;
}

template <class StorageType>
void PyramidImage<StorageType>::
SetUID(BIAS::UUID uid)
{
  typename vector<SharedPtr<Image<StorageType> > >::iterator it;
  for (it=_Images.begin(); it!=_Images.end(); it++)
    (*it)->SetUID(uid);
}

template <class StorageType>
int PyramidImage<StorageType>::
SetROI(const ROI &roi)
{
  switch (roi.GetROIType()){
  case ROI_Rows:
    {
      vector<unsigned> OrgStart, OrgEnd;
      roi.GetRows(OrgStart,OrgEnd);
      for (int i=0; i<(int)_Images.size(); i++){
        _Images[i]->GetROI()->SetROIType(ROI_Rows);
        vector<unsigned> start,end;
        const int height = _Images[i]->GetHeight();
        start.resize(height); end.resize(height);
        for (int y=0;y<height;y++){
          const int pos = (int)floor((double)y * pow(_RescaleFactor,(double)i));
          start[y] = (int)floor((double)OrgStart[pos] / 
                                pow(_RescaleFactor,(double)i));
          end[y] = (int)floor(((double)OrgEnd[pos] /
                               pow(_RescaleFactor,(double)i)));
        }
        _Images[i]->GetROI()->SetRows(start,end);
      }
    }
    break;
  default:
    BIASERR("PyramidImage<StorageType>::SetROI() unfinished");
    return -1;
  }
  return 0;
}


template <class StorageType>
int PyramidImage<StorageType>::
SetROI(unsigned minx, unsigned miny, unsigned maxx, unsigned maxy)
{
  typename vector<SharedPtr<Image<StorageType> > >::iterator it;
  for (it=_Images.begin(); it!=_Images.end(); it++){
    if ((*it)->GetROI()->SetCorners(minx, miny, maxx, maxy)!=0){
      return -1;
    }
    minx=(unsigned)floor(minx / _RescaleFactor);
    miny=(unsigned)floor(miny / _RescaleFactor);
    maxx=(unsigned)floor(maxx / _RescaleFactor);
    maxy=(unsigned)floor(maxy / _RescaleFactor);
  }
  return 0;
}

template <class StorageType>
void PyramidImage<StorageType>::
SetZero()
{
  typename vector<SharedPtr<Image<StorageType> > >::iterator it;
  for (it=_Images.begin(); it!=_Images.end(); it++){
    (*it)->SetZero();
  }
}

template <class StorageType>
void PyramidImage<StorageType>::
GetSingleImage(Image<StorageType>& im) const
{
  unsigned size=(*this).size();
  if (size==1){
    im=*(*this)[0];
  } else if (size>1) {
    //SetDebugLevel(GetDebugLevel() | D_PYRAMID_SINGLE_IM);
    register unsigned int i=0, j=0, h=0;
    unsigned int width=(*this)[0]->GetWidth()+
      (*this)[1]->GetWidth()+PYRAMID_BORDER;
    unsigned int height=(*this)[0]->GetHeight();
    unsigned int height2=0;

    //    StorageType ***pida = new StorageType**[size];
    vector<const StorageType **> pida;
    unsigned *widthstep = new unsigned[size];
    unsigned *pwidth = new unsigned[size];
    unsigned *pheight = new unsigned[size];
    unsigned cc=(*this)[0]->GetChannelCount();

    for (i=0; i<size; i++){
      pida.push_back((*this)[i]->GetImageDataArray());
      //      pida[i]=(*this)[i]->GetImageDataArray();
      widthstep[i]=(*this)[i]->GetWidthStep();
      pwidth[i]=(*this)[i]->GetWidth();
      pheight[i]=(*this)[i]->GetHeight();
      BIASCDOUT(D_PYRAMID_SINGLE_IM, "pida[i]: "<< pida[i] << "\t"
                <<(*this)[i]->GetImageDataArray()<<"\t"<<widthstep[i]
                <<"\t"<<pheight[i]<<endl);
    }
    // calculate the height of images 1- including border
    for (i=1; i<size; i++)
      height2+=(pheight[i]+PYRAMID_BORDER);
    if (height2>height) height=height2;
    // init image
    if (im.GetWidth()!=width || im.GetHeight()!=height ||
        im.GetChannelCount()!=(*this)[0]->GetChannelCount()) {
      im.Release();
      im.Init(width, height, (*this)[0]->GetChannelCount());
    }
    im.FillImageWithConstValue((StorageType)PYRAMID_BACKGROUND);
    StorageType **ida = im.GetImageDataArray();

    BIASCDOUT(D_PYRAMID_SINGLE_IM,"width: "<<width<<"\theight: "
              <<height<<endl);

    h=0;
    for (j=1; j<size; j++){
      for (i=0; i<pheight[j]; i++){
        // copy biggest image
        if (i+h<pheight[0])
          memcpy((void *)(ida[i+h]), (void *)(pida[0][i+h]), widthstep[0]);
        // copy j-biggest image
        memcpy((void *)(ida[i+h]+pwidth[0]+PYRAMID_BORDER*cc),
               (void *)(pida[j][i]),widthstep[j]);
      }
      for (i=0; i< PYRAMID_BORDER; i++)
        if (i+pheight[j]+h<pheight[0])
          memcpy((void *)(ida[i+pheight[j]+h]),
                 (void *)(pida[0][i+pheight[j]+h]),
                 widthstep[0]);
      h+=pheight[j]+PYRAMID_BORDER;
    }
    for (i=h; i<(*this)[0]->GetHeight(); i++){
      memcpy((void *)ida[i], pida[0][i], widthstep[0]);
    }
    //delete[] pida;
    delete[] widthstep;
    delete[] pwidth;
    delete[] pheight;
  }
}


template <class StorageType>
int PyramidImage<StorageType>::
CreateAdditionalLayer(unsigned int numextralayers,
                      unsigned int minImageWidth)
{
  //BIASWARN("creating "<<numextralayers<<" new levels ");
  BIASASSERT(numextralayers>0);
  BIASASSERT(!IsEmpty());
  BIASASSERT(_Images.back() != NULL);
  int res = 0;
  double pyfac = _vFactor.back();
  _Images.reserve(_Images.size()+numextralayers);
  for (unsigned int i=0; i< numextralayers; i++) {
    // check if image would be too small:
    if (double(_Images.back()->GetWidth()/_RescaleFactor<double(minImageWidth)))
      return 1;

    // create new image and downsample from top of pyramid
    Image<StorageType> *pImage = new Image<StorageType>;
    if ((res=_rescale.Filter(*(_Images.back()), *pImage))!=0){
      BIASERR("filter did not succeed "<<res);
      BIASABORT;
    }

    // update vectors and statistics
    push_back(SharedPtr<Image<StorageType> >(pImage));
    pyfac *= _RescaleFactor;
    _vFactor.push_back(pyfac);
  }
  return res;
}



////////////////////////////////////////////////////////////////////////////
//                private functions
////////////////////////////////////////////////////////////////////////////

template <class StorageType>
void PyramidImage<StorageType>::
_CreateLevels()
{
  const double factor=_RescaleFactor;

  int res = 0;
  double pyfac=1.0;
  _vFactor.push_back(pyfac);
  typename vector<SharedPtr<Image<StorageType> > >::iterator it, lastit;
  for (lastit=_Images.begin(), it=lastit+1;
         it!=_Images.end() && lastit!=_Images.end(); it++, lastit++){
    if ((*it)==NULL) (*it) = SharedPtr<Image<StorageType> >(new Image<StorageType>());
    else if (!(*it)->IsEmpty()) (*it)->Release();
    BIASASSERT(*lastit);


    if ((res=_rescale.Filter(*(*lastit), *(*it)))!=0){
      BIASERR("filter did not succeed "<<res);
      BIASABORT;
    }
    (*it)->SetOutsideROIZero();
//    cout << "setting pixels outside roi to zero--------------- " << (*it)->GetROI() << endl;
    pyfac*=factor;
    _vFactor.push_back(pyfac);
  }
  _PositionOffset = _rescale.GetPositionOffset();
}

template <class StorageType>
int PyramidImage<StorageType>::
WriteImages(const std::string& prefix) const
{
  int res=0;
  const int psize=_Images.size();
  ostringstream os;
  for (int i=0; i<psize; i++){
    os.str("");
    os << prefix << "-pl"<<i<<".mip";
    //if ((res=ImageIO::Save(os.str(), *((*this)[i])))!=0){
    if ((res=ImageIO::Save(os.str(), *((*this)[i])))!=0){
      BIASERR("error writing "<<os.str());
      break;
    }
  }
  return res;
}

template <class StorageType>
double PyramidImage<StorageType>::
GetImageValue(const double& x, const double& y,
              unsigned int scale1, int channel) const
{
  if (scale1>_Images.size()-1) scale1 = _Images.size()-1;

  // compute coordinates in smaller image
  double x1 = x, y1 = y;
  const double offset = GetPositionOffset();
  for (unsigned int s=0; s<scale1; s++) {
    // -0.5 necessary because of offset in 2^n downsampling
    x1 = (x1 - offset)/ _RescaleFactor;
    y1 = (y1 - offset)/ _RescaleFactor;
  }

/*
  unsigned int UpperLeftX, UpperLeftY, LowerRightX, LowerRightY;
  (*this)[scale1]->GetROICorners(UpperLeftX, UpperLeftY,
                                 LowerRightX, LowerRightY);
  if (x1 < (double)LowerRightX-1.0 &&
      y1 < (double)LowerRightY-1.0 &&
      x1 > double(UpperLeftX) && y1 > double(UpperLeftY)) {
*/

  // compute value only in image (no pyramid)
  if ((*this)[scale1]->GetROI()->CheckBilinearInterpolation(x1, y1)) {
    return (*this)[scale1]->
      BilinearInterpolation(x1, y1, (unsigned short int)channel);
  }
  return StorageType(0);
}


template <class StorageType>
int PyramidImage<StorageType>::
GetTrilinearImageValue(const double& xsource,
                       const double& ysource,
                       const double& scale,
                       double& T,
                       int channel) const {
  BIASASSERT(scale>=0.0);
  // compute a weight for the smaller and a weight for the larger pyramid level
  // then compute the colors at these levels and average them according
  // to their weights

  double value1 =  0, value2 = 0;
  double scale_error = scale - int(scale);

  int scale1 = int(scale);
  if (scale1>=(int)_Images.size()-1) {
    // we might get aliasing, pyramid is too small !!!
    scale1 = _Images.size()-1;
    scale_error = 0.0;
  } else if (scale1<0) {
    // use max resolution image, dont need to upsample in pyramid
    scale1 = 0;
    scale_error = 0.0;
  }

  // compute coordinates in smaller image
  double x1 = xsource, y1 = ysource;
  const double offset = GetPositionOffset();
  for (int s=0; s<scale1; s++) {
    // -0.5 necessary because of offset in 2^n downsampling
    x1 = (x1 - offset) / _RescaleFactor;
    y1 = (y1 - offset) / _RescaleFactor;
  }
  // if scale 1 fits perfectly
  // we dont need scale2, set it to the same value to avoid access violations
  const int scale2 = (scale_error==0.0)? scale1 : scale1 + 1;
  const double x2  = (scale_error==0.0)? x1     : (x1 - offset)/_RescaleFactor;
  const double y2  = (scale_error==0.0)? y1     : (y1 - offset)/_RescaleFactor;

  double weight1 = 1.0 - scale_error, weight2 = scale_error;

/*
  unsigned int UpperLeftX, UpperLeftY, LowerRightX, LowerRightY;
  (*this)[scale1]->GetROICorners(UpperLeftX, UpperLeftY,
                                 LowerRightX, LowerRightY);

  if (x1 < (double)LowerRightX-1.0 &&
      y1 < (double)LowerRightY-1.0 &&
      x1 > double(UpperLeftX) && y1 > double(UpperLeftY) ) {
*/

  // compute value only in image (no pyramid)
  if ((*this)[scale1]->GetROI()->CheckBilinearInterpolation(x1, y1)) {
    value1 = (*this)[scale1]->
      BilinearInterpolation(x1, y1, (unsigned short int)channel);
  } else {
    weight1 = 0;
  }

  // compute second scale only if necessary
  if (weight2>0.0)
  {
/*
    (*this)[scale2]->GetROICorners(UpperLeftX, UpperLeftY,
                                   LowerRightX, LowerRightY);

    if (x2 < (double)LowerRightX-1.0 &&
        y2 < (double)LowerRightY-1.0 &&
        x2 > double(UpperLeftX) && y2 > double(UpperLeftY) ) {
*/
    // compute value only in image (no pyramid)
    if ((*this)[scale2]->GetROI()->CheckBilinearInterpolation(x2, y2)) {
      value2 = (*this)[scale2]->
        BilinearInterpolation(x2, y2, (unsigned short int)channel);
    } else {
      weight2 = 0;
    }
  }
  // check if we computed any value at all
  const double weightsum = weight1 + weight2;

  // this is the tri in trilinear:
  T = (weight1 * value1 + weight2 * value2) / weightsum;
  if (weightsum>0.5) return 0;
  return -1;
}


// (1.0/log(2.0));
#define ONE_OVER_LOG2 1.442695

#define SQRT_2PI 2.5066

template <class StorageType>
int PyramidImage<StorageType>::
GetAnisotropicImageValue(const double& xsource,
                         const double& ysource,
                         const Matrix2x2<double>& thecov,
                         double& T,
                         unsigned int channel) const {
  Matrix2x2<double> Smoother(thecov);
  // compute eigenvalues of jacobian, which determine the largest step sizes
  Vector2<double> ev1, ev2;

  // compute eigenvalues = squared step sizes in source
  double d1 = 0.0, d2 = 0.0;
  Smoother.EigenvalueDecomposition(d1, ev1, d2, ev2);


  //cout<<"Eigenvalue decomposition: "<<d1<<" "<<ev1<<" "<<d2<<" "<<ev2<<endl;

  // remove square, require min eigenvalue 1, because we have no finer image
  // information
  // (should rarely happen for meaningful mappings)
  if (d1<0.25) d1 = 0.5; else d1 = sqrt(d1);
  if (d2<0.25) d2 = 0.5; else d2 = sqrt(d2);

  // make sure d1 is not smaller than d2
  if (d1<d2) {
    std::swap(d1, d2);
    ev1[0] = ev2[0];
    ev1[1] = ev2[1];
  }

  // compute pyramid levels
  const double smallscale = log(d2)*ONE_OVER_LOG2+1.0;
  const double largescale = log(d1)*ONE_OVER_LOG2+1.0;

  if (largescale<=0.1) {
    // bilinear is sufficient since we are only enlarging
    if (!(*this)[0]->IsInROI(xsource,ysource)) return -1;
    T = GetImageValue(xsource, ysource, 0, channel);
    return 0;
  }

  // compute real anisotropy which needs to be filtered away:
  double smoothsigma = sqrt(0.25*(d1/d2*d1/d2)-0.25);

  //cout<<"required anisotropic smoothing is "
  //    <<smoothsigma<<endl;

  // at small pyramid level, when do we get "new" information
  const double stepsize = pow(2.0, smallscale);

  // use trilinear if nearly isotropic
  if (GetTrilinearImageValue(xsource, ysource,
                             smallscale, T, channel)!=0) return -1;

  // how large a mask do we need at this scale
  const int MaskHalfSize = (int)(floor(2.5*smoothsigma/stepsize));
  // isotropic ?
  if (MaskHalfSize<1) return 0;

  // in which direction do we smooth
  const double stepdirx = ev1[0]*stepsize;
  const double stepdiry = ev1[1]*stepsize;
  // precompute constant
  const double invsigmasigma = -0.5 / (smoothsigma*smoothsigma);

  // really used weights and values
  // already have center value with relative weight 1
  double weightsum(1.0), valuesum(T);

  for (int i=1; i<=MaskHalfSize; i++) {
    // current relative weight
    const double curWeight = exp(i*i*invsigmasigma);

    double value = 0;
    if (GetTrilinearImageValue(xsource-double(i)*stepdirx,
                               ysource-double(i)*stepdiry,
                               smallscale, value, channel)==0) {
      weightsum += curWeight;
      valuesum += curWeight * value;
    } else break;
    if (GetTrilinearImageValue(xsource+double(i)*stepdirx,
                               ysource+double(i)*stepdiry,
                               smallscale, value, channel)==0) {
      weightsum += curWeight;
      valuesum += curWeight * value;
    } else break;
  }
  // check if we covered at least 50% of the gaussian
  if (weightsum / (smoothsigma*SQRT_2PI) < 0.5) return -1;
  T = valuesum/weightsum;

  return 0;
}

template <class StorageType>
void PyramidImage<StorageType>::
Dump(std::ostream& os) const
{
  os << "pyramid image of size "<<_Images.size()<<" with factor "
     << _RescaleFactor<<":\n";
  for (unsigned i=0; i<_Images.size(); i++){
    os << _Images[i]<<" : ";
    if (_Images[i]){
      os << _Images[i]->GetWidth()<<" x "<<_Images[i]->GetHeight();
    }
    os << endl;
  }
}


template <class StorageType>
bool PyramidImage<StorageType>::
IsInROI(double x, double y, int layer) const {

  for (int l=layer;l>0; l--) {
    x = (x - _PositionOffset)/ _RescaleFactor;
    y = (y - _PositionOffset)/ _RescaleFactor;
  }
  return (*this)[layer]->IsInROI(x,y);
}


//////////////////////////////////////////////////////////////////////////
//                 instantiation
//////////////////////////////////////////////////////////////////////////
namespace BIAS{
template class PyramidImage<unsigned char>;
template class PyramidImage<float>;

#ifdef BUILD_IMAGE_INT
template class PyramidImage<int>;
#endif
#ifdef BUILD_IMAGE_CHAR
template class PyramidImage<char>;
#endif
#ifdef BUILD_IMAGE_SHORT
template class PyramidImage<short>;
#endif
#ifdef BUILD_IMAGE_USHORT
template class PyramidImage<unsigned short>;
#endif
#ifdef BUILD_IMAGE_UINT
template class PyramidImage<unsigned int>;
#endif
#ifdef BUILD_IMAGE_DOUBLE
template class PyramidImage<double>;
#endif
}