[Insight-users] Q: Approaches for convolution with a general kernel

[Gunnar Farnebäck]

Stefan wrote:
> I have a relatively small-sized 2D kernel (no more than 40x40 pixels) 
> that I want to use for convolution with images of about 1000x1000 
> pixels. I can either calculate this by brute force or do the convolution 
> in Fourier transform space. I do not find a lot of references to this 
> problem in the documentation, but it seems like a common enough problem 
> that I must have missed something. My current plan is to work from the 
> Sobel edge-detection code and develop a direct-calculation routine from 
> that. Any suggestions for alternative places to start?

As Luis said 40x40 is quite a large kernel. But if it is Cartesian
separable and preferrably somewhat smaller than 40x40, two direct 1D
convolutions may well be faster than an fft solution. In that case you
may find the attached itkDirectionalOperator class useful. I can try
to extract a small usage example if needed.

/Gunnar Farnebäck

-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkDirectionalOperator.h,v $
  Language:  C++
  Date:      $Date: 2003/09/10 14:29:09 $
  Version:   $Revision: 1.22 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkDirectionalOperator_h
#define __itkDirectionalOperator_h

#include "itkNeighborhoodOperator.h"
namespace itk {

/**
 * \class DirectionalOperator
 * \brief A NeighborhoodOperator for one dimensional directional
 * kernels with explicitly specified coefficients.
 *
 * DirectionalOperator can be used to perform directional filtering
 * by taking its inner product with a Neighborhood
 * (NeighborhoodIterator) that is swept across an image region.
 *
 * DirectionalOperator takes two parameters:
 *
 * (1) The direction of the operator.
 *
 * (2) The coefficients of the filter.
 *
 * After setting the parameters it is necessary to call the
 * CreateDirectional() method, before the operator can be used.
 *
 * \sa NeighborhoodOperator
 * \sa NeighborhoodIterator
 * \sa Neighborhood
 *
 * \ingroup Operators
 */
template<class TPixel, unsigned int VDimension = 2,
         class TAllocator = NeighborhoodAllocator<TPixel> >
class ITK_EXPORT DirectionalOperator
  : public NeighborhoodOperator<TPixel, VDimension, TAllocator>
{
public:
  /** Standard class typedefs. */
  typedef DirectionalOperator Self;
  typedef NeighborhoodOperator<TPixel, VDimension, TAllocator>  Superclass;
  
  /** Constructor. */
  DirectionalOperator()
  {
    m_Coefficients = std::vector<double>(1, 0.0);
  }
  
  /** Copy constructor. */
  DirectionalOperator(const Self &other)
    : NeighborhoodOperator<TPixel, VDimension, TAllocator>(other)
  {
    m_Coefficients = other.m_Coefficients;
  }

  /** Assignment operator. */
  Self &operator=(const Self &other)
  {
    Superclass::operator=(other);
    m_Coefficients = other.m_Coefficients;
    return *this;
  }
  
  /** Sets the kernel coefficients. */
  void SetCoefficients(const std::vector<double> coefficients)
  {
    m_Coefficients = coefficients;
  }

  /** Returns the kernel coefficients. */
  std::vector<double> GetCoefficients() const
  {
    return m_Coefficients;
  }
  
  /** Prints some debugging information. */
  void PrintSelf(std::ostream &os, Indent i) const
  {
    os << i << "DirectionalOperator { this=" << this
       << ", m_Coefficients =";
    for (int k = 0; k < m_Coefficients.size(); k++)
      os << " " << m_Coefficients[k];
    os << "} " << std::endl;
    Superclass::PrintSelf(os, i.GetNextIndent());
  }
  
protected:
  typedef typename Superclass::CoefficientVector CoefficientVector;

  /** Calculates operator coefficients. */
  CoefficientVector GenerateCoefficients()
  {
    return m_Coefficients;
  }

  /** Arranges coefficients spatially in the memory buffer. */
  void Fill(const CoefficientVector& coeff)
  {
    this->FillCenteredDirectional(coeff);
  }

private:
  /** Filter coefficients. */
  CoefficientVector m_Coefficients;
  
  /** For compatibility with itkWarningMacro. */
  const char *GetNameOfClass()
  {
    return "itkDirectionalOperator";
  }
  
};

} // namespace itk


#ifndef ITK_MANUAL_INSTANTIATION
#endif

#endif