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

#include <Base/Image/ImageIO.hh>
#include <Base/Image/ImageConvert.hh>

#include <Base/Math/Vector2.hh>
#include <algorithm>

using namespace BIAS;
using namespace std;

template <class InputStorageType, class OutputStorageType> 
Bilateral<InputStorageType, OutputStorageType>::Bilateral()
  : FilterNToN<InputStorageType,OutputStorageType>()
{
  _lastBilateralSize=-1;
  _lastSecondSize=-1;
  _secondSize=0;
  _BilateralSize=0;
  _GaussSigma= _BilateralSigma=0.7; 
  _ignoreValue = (InputStorageType)0.0;
  SetSize(2,2);
}


template <class InputStorageType, class OutputStorageType> 
Bilateral<InputStorageType, OutputStorageType>::
Bilateral(const Bilateral<InputStorageType, OutputStorageType>& other) 
  : FilterNToN<InputStorageType, OutputStorageType>(other)
{

   _BilateralSize = other._BilateralSize;
  _secondSize = other._secondSize;
    _GaussSigma = other._GaussSigma;
  _ignoreValue = other._ignoreValue;
  _BilateralSigma= other._BilateralSigma;

  SetSize(_BilateralSize,_secondSize);
}


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

template <class InputStorageType, class OutputStorageType> 
void Bilateral<InputStorageType, OutputStorageType>::
SetSize(int newsize, int secondsize)
{
  if(newsize != _lastBilateralSize || secondsize!=_lastSecondSize){
    _BilateralSize = newsize;
    _lastBilateralSize =_BilateralSize;
    if(secondsize != -1){
      _secondSize = secondsize;
      _lastSecondSize = _secondSize;
    }
    else{
      _secondSize = _BilateralSize;
      _lastSecondSize = _secondSize;
    }    
  }
  _CalculateKernels(_GaussSigma);
}

template <class InputStorageType, class OutputStorageType>
void Bilateral<InputStorageType, OutputStorageType>::
_CalculateKernels(double gaussSigma)
{
  _GaussSigma = gaussSigma;
  //double fac= 1.0/(2.0*gaussSigma*gaussSigma);

  register double g=0.0;
   _gaussFilterMask.Init(_BilateralSize*2+1,_secondSize*2+1,1);
  float **fmida = _gaussFilterMask.GetImageDataArray();
  //create 2D gauss mask
  Vector2<double> zero(0,0);
  for(int xh=-_BilateralSize;xh<=_BilateralSize;xh++){      
    for(int yh=-_secondSize;yh<=_secondSize;yh++){
      Vector2<double> p(xh,yh); 
      double dist = zero.Dist(p)*zero.Dist(p);
      g  = exp(double(-0.5*(double)(dist/gaussSigma*gaussSigma)));           //(xh*xh)+(yh*yh))));
      fmida[yh+_secondSize][xh+_BilateralSize] = float(g);
    }    
  }
  //BIAS::ImageIO::Save("filtermask.mip",_gaussFilterMask);
}

////////////////////////////////////////////////////////////////
template <class InputStorageType, class OutputStorageType> 
int Bilateral<InputStorageType, OutputStorageType>::
Filter(const Image<InputStorageType>& src, Image<OutputStorageType>& dst) 
{
  if (src.GetChannelCount()!=1) {
    return FilterColorImg(src,dst);   
  }
  
  Image<OutputStorageType> tmp;
  OutputStorageType **dstida = dst.GetImageDataArray();
  const InputStorageType **srcida = src.GetImageDataArray();
  bool identimg = false;
 
  if ((void*)src.GetImageData() != (void*)dst.GetImageData()) {   
    if (!src.SamePixelAndChannelCount(dst)) {
      if (!dst.IsEmpty())
        dst.Release();
      dst.Init(src.GetWidth(), src.GetHeight(), 1);
      dstida = dst.GetImageDataArray();
    }
  } else {
    identimg = true;
    tmp.Init(src.GetWidth(), src.GetHeight(), 1);
    dstida = tmp.GetImageDataArray();
  }

  register double g=0.0,gb=0.0;
  // double facB= 1.0/(2.0*_BilateralSigma*_BilateralSigma);
  float ** fmida =  _gaussFilterMask.GetImageDataArray();
  int y=0,x=0,xh=0,yh=0;
  int width = (int)src.GetWidth();
  int height = (int)src.GetHeight();
  double sumGB=0.0,sumGBMulI=0.0;
  double scrSq=0.0;
  double bilSigSq = _BilateralSigma*_BilateralSigma;

  // BIAS::Image<float> bilFilter(width*(_BilateralSize*2+1),height*(_secondSize*2+1),1) ;
  // float ** bilida = bilFilter.GetImageDataArray();

#ifdef BIAS_HAVE_OPENMP
#pragma omp parallel for private(y,x,xh,yh,g,gb,scrSq,sumGB,sumGBMulI)
#endif
  for ( y=0; y < height; y++) {
    for ( x=0; x < width; x++) {   
      if(srcida[y][x] != 0.0){
        sumGB =sumGBMulI = 0.0;
        unsigned count=0;
        for(yh = -_secondSize; yh <= _secondSize; yh++){   
          for(xh = -_BilateralSize; xh <= _BilateralSize; xh++){
            if((y+yh >= 0) && (y+yh < height) && 
               (x+xh >= 0) && (x+xh < width)  &&
               srcida[y][x] != _ignoreValue && srcida[y+yh][x+xh] != _ignoreValue)
              {
                //standard gauss
                g  = (double)fmida[_secondSize+yh][_BilateralSize+xh];
      
                //bilateral term;         
                scrSq = (double)(srcida[y][x]-srcida[y+yh][x+xh]);               
                scrSq *=scrSq;
                gb = exp(double(-0.5*(scrSq/bilSigSq)) );

                //bilida[y*(_secondSize*2+1)+yh][x*(_BilateralSize*2+1)+xh] = gb;
            
                sumGB+= g*gb;              
                sumGBMulI += g*gb*(double)srcida[y+yh][x+xh]; 
                count++;
              }
          }
        } 
     
        if(sumGB != 0.0 && count != 0 && !BIAS_ISNAN(sumGB) && !BIAS_ISINF(sumGB)){
          //this is a data dependent normalization, 
          //better divide by count.
          dstida[y][x] = (OutputStorageType)(sumGBMulI / sumGB);                
          //dstida[y][x] = (OutputStorageType)(sumGBMulI / count);      
        }
        else{       
          dstida[y][x] = (OutputStorageType) srcida[y][x];
        }
      }
    }
  }
  // BIAS::ImageIO::Save("bilateralmask.mip",bilFilter);
  if (identimg)  dst.StealImage(tmp);
  return 0;
}

////////////////////////////////////////////////////////////////
template <class InputStorageType, class OutputStorageType> 
int Bilateral<InputStorageType, OutputStorageType>::
Filter(const Image<InputStorageType>& src, 
       const Image<InputStorageType>& support, 
       Image<OutputStorageType>& dst) 
{
  if (!src.SamePixelAndChannelCount(support)) {
    BIASERR("Source and support image do not have same pixel and channel count!");
    return -1;
  }

  if (src.GetChannelCount()!=1) {
    return FilterColorImg(src,support,dst);   
  }
  
  Image<OutputStorageType> tmp;
  OutputStorageType **dstida = dst.GetImageDataArray();
  const InputStorageType **srcida = src.GetImageDataArray();
  const InputStorageType **suppida= support.GetImageDataArray();
  bool identimg = false;
 
  if ((void*)src.GetImageData() != (void*)dst.GetImageData()) {   
    if (!src.SamePixelAndChannelCount(dst)) {
      if (!dst.IsEmpty())
        dst.Release();
      dst.Init(src.GetWidth(), src.GetHeight(), 1);
      dstida = dst.GetImageDataArray();
    }
  } else {
    identimg = true;
    tmp.Init(src.GetWidth(), src.GetHeight(), 1);
    dstida = tmp.GetImageDataArray();
  }

  register double g=0.0,gb=0.0; 
  float ** fmida =  _gaussFilterMask.GetImageDataArray();
  int y=0,x=0,xh=0,yh=0;
 
  int width = (int)src.GetWidth();
  int height = (int)src.GetHeight();
  double sumGB=0.0,sumGBMulI=0.0;
  double scrSq=0.0;
  double bilSigSq = _BilateralSigma*_BilateralSigma;
    
#ifdef BIAS_HAVE_OPENMP
#pragma omp parallel for private(y,x,xh,yh,g,gb,scrSq,sumGB,sumGBMulI)
#endif
  for ( y=0; y < height; y++) {
    for ( x=0; x < width; x++) {   
      if(srcida[y][x] != 0.0){
        sumGB =sumGBMulI = 0.0;
      
        for(xh=-_BilateralSize;xh<=_BilateralSize;xh++){
          for(yh=-_secondSize;yh<=_secondSize;yh++){
            if(y+yh >= 0 && y+yh < height &&
               x+xh >= 0 && x+xh < width &&
               srcida[y][x] != _ignoreValue && srcida[y+yh][x+xh] != _ignoreValue)
              {
                //standard gauss             
                g  = (double)fmida[_secondSize+yh][_BilateralSize+xh];

                //bilateral term;         
                scrSq = (double)(suppida[y][x]-suppida[y+yh][x+xh]);               
                scrSq *=scrSq;
                gb = exp(double(-0.5*(scrSq/bilSigSq)) );
              
                //add up
                sumGB+= g*gb;              
                sumGBMulI += g*gb*srcida[y+yh][x+xh];
              }
          }     
        }   
      
        if(sumGB != 0.0 && !BIAS_ISNAN(sumGB) && !BIAS_ISINF(sumGB)){
          dstida[y][x] = (OutputStorageType)(sumGBMulI / sumGB);      
        }
        else{ 
          dstida[y][x] = (OutputStorageType) srcida[y][x];
        }
      }
    }
  }

  if (identimg)  dst.StealImage(tmp);
  return 0;
}

template <class InputStorageType, class OutputStorageType> 
int Bilateral<InputStorageType, OutputStorageType>::
FilterInt(const Image<InputStorageType>& src, Image<OutputStorageType>& dst) 
{
  BIASERR("FilterInt makes no sense here, using general Filter interface.");
  return Filter(src, dst);
}



template <class InputStorageType, class OutputStorageType> 
int Bilateral<InputStorageType, OutputStorageType>::
FilterFloat(const Image<InputStorageType>& src, Image<OutputStorageType>& dst) 
{
  BIASERR("FilterFloat makes no sense here, using general Filter interface.");
  return Filter(src, dst);
}


template <class InputStorageType, class OutputStorageType> 
void Bilateral<InputStorageType, OutputStorageType>::
GetBordersValid_(int &border_x, int &border_y) const 
{
  if (_secondSize > _BilateralSize) {
    border_x = border_y = _secondSize;
  } else {
    border_x = border_y = _BilateralSize;
  }
}

template <class InputStorageType, class OutputStorageType> 
int Bilateral<InputStorageType, OutputStorageType>::
FilterColorImg(const Image<InputStorageType>& src,
               Image<OutputStorageType>& dst) 
{
  Image<InputStorageType> srcCopy, temp;
 
  int width = src.GetWidth();
  int height = src.GetHeight();
  dst =  Image<OutputStorageType>(width, height, 3, true);
   
  Image<InputStorageType> red(width, height, 1);
  Image<OutputStorageType> redG(width, height, 1);
  Image<InputStorageType> green(width, height, 1);
  Image<OutputStorageType> greenG(width, height, 1);
  Image<InputStorageType> blue(width, height, 1);
  Image<OutputStorageType> blueG(width, height, 1);

  /// use planar image in order to apply median filter seperately to each channel for now
  if (!((src.GetColorModel() == ImageBase::CM_RGB) && src.IsPlanar())) {
    temp = src;
    ImageConvert::Convert(temp, srcCopy, ImageBase::CM_RGB, true);
    red.GetChannel(srcCopy, 0);
    green.GetChannel(srcCopy, 1);
    blue.GetChannel(srcCopy, 2);
  }else{
    red.GetChannel(src, 0);
    green.GetChannel(src, 1);
    blue.GetChannel(src, 2);
  }

  Filter(red, redG);
  Filter(green, greenG);
  Filter(blue, blueG);

  BIAS::ImageBase::SetChannel(dst,0,redG.GetImageData());
  BIAS::ImageBase::SetChannel(dst,1,greenG.GetImageData());
  BIAS::ImageBase::SetChannel(dst,2,blueG.GetImageData());


  return 0 ;
}

template <class InputStorageType, class OutputStorageType> 
int Bilateral<InputStorageType, OutputStorageType>::
FilterColorImg(const Image<InputStorageType>& src,
               const Image<InputStorageType>& support,
               Image<OutputStorageType>& dst) 
{
  Image<InputStorageType> srcCopy,supportCopy;
  int width = src.GetWidth();
  int height = src.GetHeight();
  dst =  Image<OutputStorageType>(width, height, 3, true);

  Image<InputStorageType> red(width, height, 1);
  Image<InputStorageType> redSupp(width, height, 1);
  Image<OutputStorageType> redG(width, height, 1);
  Image<InputStorageType> green(width, height, 1);
  Image<InputStorageType> greenSupp(width, height, 1);
  Image<OutputStorageType> greenG(width, height, 1);
  Image<InputStorageType> blue(width, height, 1);
  Image<InputStorageType> blueSupp(width, height, 1);
  Image<OutputStorageType> blueG(width, height, 1);
  
  /// use planar image in order to apply median filter seperately to each channel for now
  if (!((src.GetColorModel() == ImageBase::CM_RGB) && src.IsPlanar())) {
    Image<InputStorageType> temp = src;
    ImageConvert::Convert(temp, srcCopy, ImageBase::CM_RGB, true);
    red.GetChannel(srcCopy, 0);
    green.GetChannel(srcCopy, 1);
    blue.GetChannel(srcCopy, 2);
  }else{
    red.GetChannel(src, 0);
    green.GetChannel(src, 1);
    blue.GetChannel(src, 2);
  }
  
  if (!((support.GetColorModel() == ImageBase::CM_RGB) && support.IsPlanar())) {
    Image<InputStorageType> temp = support;
    ImageConvert::Convert(temp, supportCopy, ImageBase::CM_RGB, true);
    redSupp.GetChannel(supportCopy, 0);
    greenSupp.GetChannel(supportCopy, 1);
    blueSupp.GetChannel(supportCopy, 2);
  }else{
    redSupp.GetChannel(support, 0);
    greenSupp.GetChannel(support, 1);
    blueSupp.GetChannel(support, 2);
  }

 
  Filter(red, redSupp,redG);
  Filter(green, greenSupp, greenG);
  Filter(blue, blueSupp, blueG);


  BIAS::ImageBase::SetChannel(dst,0,redG.GetImageData());
  BIAS::ImageBase::SetChannel(dst,1,greenG.GetImageData());
  BIAS::ImageBase::SetChannel(dst,2,blueG.GetImageData());  

  return 0 ;
}
//////////////////////////////////////////////////////////////////////////
//                 instantiation
//////////////////////////////////////////////////////////////////////////

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