filterbank.cc

用于计算矩阵的特征值,及矩阵的其他运算.可以用与稀疏矩阵

// Copyright (C) 2002 Charless C. Fowlkes <[e-mail]>
// Copyright (C) 2002 David R. Martin <[e-mail]>
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of the
// License, or (at your option) any later version.
//
// This program 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
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
// 02111-1307, USA, or see http://www.gnu.org/copyleft/gpl.html.


#include <math.h>
#include <assert.h>
#include "util.hh"
#include "string.hh"
#include "message.hh"
#include "image.hh"
#include "array.hh"
#include "filterbank.hh"

namespace Group
{

  ////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // some static functions useful for building filters
  //

  //
  //return a new image that is zero mean
  //
  static void makeZeroMean(Util::Image& im)
  {
    int width = im.size(0);
    int height = im.size(1);

    float sum = 0;
    for (int i = 0; i < width; i++)
    {
      for (int j = 0; j < height; j++)
      {
        sum += im(i,j);
      }
    }
                                            
    if (width*height > 0)
    {
      float mean = (sum / ((float)width*(float)height));
      im -= mean;
    }
  }

  //
  //return an image that is unit L1norm
  //
  static void makeL1Normalized(Util::Image& im)
  {
    float norm = 0;
    int width = im.size(0);
    int height = im.size(1);
    for (int i = 0; i < width; i++)
    {
      for (int j = 0; j < height; j++)
      {
          norm += fabs(im(i,j));
      }
    }
    if (norm > 0)
    {
      im /= norm;
    }
  }


  //
  // compute 1 dimensional DFT for odd length signal
  //
  static void DFT(const Util::Array1D<float> re, const Util::Array1D<float> im,
            Util::Array1D<float>& dft_re, Util::Array1D<float>& dft_im)
  {
    int N = re.size(0);
    assert((N % 2) == 1);
    int halfN = (int)(N/2);
    float Ninv = 1/(float)N;

    dft_re.resize(N);
    dft_im.resize(N);
    dft_re.init(0);
    dft_im.init(0);

    for(int k = -halfN; k <= halfN; k++)
    {
      float d = -k * 2 * M_PI * Ninv;
      for(int j = 0; j < N; j++)
      {
        float c = cos(d * (j+1));
        float s = sin(d * (j+1));
        dft_re(k + halfN) += (re(j)*c - im(j)*s);
        dft_im(k + halfN) += (re(j)*s + im(j)*c);
      }
    }
    dft_re *= Ninv;
    dft_im *= Ninv;
  }

  //
  // compute inverse of 1 dimensional DFT for odd length signal
  //
  static void IDFT(const Util::Array1D<float> re, const Util::Array1D<float> im,
            Util::Array1D<float>& dft_re, Util::Array1D<float>& dft_im)
  {
    int N = re.size(0);
    assert((N % 2) == 1);
    int halfN = (int)(N/2);
    float Ninv = 1/(float)N;

    dft_re.resize(N);
    dft_im.resize(N);
    dft_re.init(0);
    dft_im.init(0);

    for(int j = 0; j < N; j++)
    {
      float d = (j+1) * 2.0 * M_PI * Ninv;
      for(int k = -halfN; k <= halfN; k++)
      {
        float c = cos(d * k);
        float s = sin(d * k);
        dft_re(j) += (re(k+halfN)*c - im(k+halfN)*s);
        dft_im(j) += (re(k+halfN)*s + im(k+halfN)*c);
      }
    }
  }

  //
  // compute the hilbert transform of a discrete sequence of odd length. 
  // operates in place on yval
  //
  static void computeHilbert (Util::Array1D<float>& yval)
  {
      int N = yval.size(0);
      assert((N % 2) == 1);

      Util::Array1D<float> re = yval;
      Util::Array1D<float> im(N);
      im.init(0);

      Util::Array1D<float> dft_re(N);
      Util::Array1D<float> dft_im(N);
      
      DFT(re,im,dft_re,dft_im);

      //double positive frequencies
      //zero out negative frequencies
      //leave the dc component
      for (int i = 0; i < (N/2); i++)
      {
        dft_re(i) *= 2.0;    
        dft_im(i) *= 2.0;    
      }
      for (int i = (N/2)+1; i < N; i++)
      {
        dft_re(i) = 0.0;    
        dft_im(i) = 0.0;    
      }

      IDFT(dft_re,dft_im,im,yval);
  }

  ////////////////////////////////////////////////////////////////////////////////////////////////
  //
  // create an empty filterbank
  //
  FilterBank::FilterBank ()
  {
      Util::Message::debug("creating empty filterbank");
      m_nfilters = 0;
      m_filters.resize(m_nfilters);
      m_sigmasX.resize(m_nfilters);
      m_sigmasY.resize(m_nfilters);
      m_orients.resize(m_nfilters);
  }


  //
  // create a filterbank
  //
  FilterBank::FilterBank(const int numOrientations, const int numScales,
                          const float startScale, const float scaleFactor)
  {
      m_nscales = numScales;
      m_norient = numOrientations;
      m_nfilters = (2*numOrientations+1) * numScales;
      m_filters.resize(m_nfilters);
      m_sigmasX.resize(m_nfilters);
      m_sigmasY.resize(m_nfilters);
      m_orients.resize(m_nfilters);
      int filterIndex = 0;

      float sigmaRatio = 3.0;
      float sigma = startScale;
      for (int scale = 0; scale < numScales; scale++)
      {
        float sigmaX = sigma;
        float sigmaY = sigma / sigmaRatio;
        Util::Message::debug(Util::String("sx = %f, sy = %f",sigmaX,sigmaY),1);
        Util::Message::startBlock(m_norient,"filterbank construction");

        //create oriented derivative of gaussian filter pairs
        for (int i = 0; i < m_norient; i++)
        {
            Util::Message::stepBlock();
            float orientation = ((i * M_PI) / (float) m_norient);

            createGaussianKernel(sigmaX, sigmaY, 3, orientation, 2, true, m_filters(filterIndex));
            m_sigmasX(filterIndex) = sigmaX;
            m_sigmasY(filterIndex) = sigmaY;
            m_orients(filterIndex) = orientation;
            filterIndex++;

            createGaussianKernel(sigmaX, sigmaY, 3, orientation, 2, false, m_filters(filterIndex));
            m_sigmasX(filterIndex) = sigmaX;
            m_sigmasY(filterIndex) = sigmaY;
            m_orients(filterIndex) = orientation;
            filterIndex++;
        }
        Util::Message::endBlock();

        //create the center-surround filter
        //make sure it's zero-mean and L1 normalized
        float support = 4;
        float sigmaC = sigma / sigmaRatio;
        float sigmaS = sigmaC * scaleFactor;
        float supportC = support * sigmaS / sigmaC;
        float supportS = support;

        Util::Image center;
        createGaussianKernel(sigmaC, sigmaC, supportC, 0, 0, false,center);
        Util::Image surround;
        createGaussianKernel (sigmaS, sigmaS, supportS, 0, 0, false,surround);
        int hC = center.size(1);
        int wC = center.size(0);
        int hS = surround.size(1);
        int wS = surround.size(0);
        assert(hC == hS);
        assert(wC == wS);

        m_filters(filterIndex) = surround-center;
        makeZeroMean(m_filters(filterIndex));
        makeL1Normalized(m_filters(filterIndex));

        m_sigmasX(filterIndex) = sigmaS;
        m_sigmasY(filterIndex) = sigmaS;
        m_orients(filterIndex) = 0;
        filterIndex++;

        sigma = sigma * scaleFactor;
      }
      assert (filterIndex == m_nfilters);
  }

  FilterBank::~FilterBank ()
  {
    Util::Message::debug("destroying filterbank",1);
  }


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

  const Util::Image& FilterBank::getFilter (int i) const
  {
    assert (i >= 0 && i < m_nfilters);
    return m_filters(i);
  }

  int FilterBank::getNumFilters () const
  {
    return m_nfilters;
  }

  float
  FilterBank::getSigmaX (int index) const
  {
    assert (index >= 0 && index < m_nfilters);
    return m_sigmasX(index);
  }