itkVectorGradientMagnitudeImageFilter.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkVectorGradientMagnitudeImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2003/01/30 21:02:04 $
  Version:   $Revision: 1.5 $

  Copyright (c) 2002 Insight 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 __itkVectorGradientMagnitudeImageFilter_h
#define __itkVectorGradientMagnitudeImageFilter_h

#include "itkConstNeighborhoodIterator.h"
#include "itkNeighborhoodIterator.h"
#include "itkImageToImageFilter.h"
#include "itkImage.h"
#include "itkVector.h"
#include "vnl/vnl_matrix.h"
#include "vnl/algo/vnl_symmetric_eigensystem.h"
#include "vnl/vnl_math.h"

namespace itk
{
template < typename TInputImage,
           typename TRealType = float,
           typename TOutputImage = Image< TRealType,
                            ::itk::GetImageDimension<TInputImage>::ImageDimension >
>
class ITK_EXPORT VectorGradientMagnitudeImageFilter :
    public ImageToImageFilter< TInputImage, TOutputImage >
{
public:
  typedef VectorGradientMagnitudeImageFilter Self;
  typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
  typedef SmartPointer<Self> Pointer;
  typedef SmartPointer<const Self>  ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(VectorGradientMagnitudeImageFilter, ImageToImageFilter);
  
  typedef typename TOutputImage::PixelType OutputPixelType;
  typedef typename TInputImage::PixelType  InputPixelType;

  typedef TInputImage  InputImageType;
  typedef TOutputImage OutputImageType;
  typedef typename InputImageType::Pointer  InputImagePointer;
  typedef typename OutputImageType::Pointer OutputImagePointer;
  
  itkStaticConstMacro(ImageDimension, unsigned int,
                      TOutputImage::ImageDimension);
  
  itkStaticConstMacro(VectorDimension, unsigned int,
                      InputPixelType::Dimension);

  typedef TRealType RealType;
  typedef Vector<TRealType, ::itk::GetVectorDimension<InputPixelType>::VectorDimension> RealVectorType;
  typedef Image<RealVectorType, ::itk::GetImageDimension<TInputImage>::ImageDimension>  RealVectorImageType;
  

  typedef ConstNeighborhoodIterator<RealVectorImageType> ConstNeighborhoodIteratorType;
  
  typedef typename Superclass::OutputImageRegionType OutputImageRegionType;

  virtual void GenerateInputRequestedRegion() throw(InvalidRequestedRegionError);

  void SetUseImageSpacingOn()
  { this->SetUseImageSpacing(true); }

  void SetUseImageSpacingOff()
  { this->SetUseImageSpacing(false); }

  void SetUseImageSpacing(bool);                         
  itkGetMacro(UseImageSpacing, bool);

  void SetDerivativeWeights(TRealType data[]);
  itkGetVectorMacro(DerivativeWeights, const TRealType, itk::GetImageDimension<TInputImage>::ImageDimension);

  itkSetVectorMacro(ComponentWeights, TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
  itkGetVectorMacro(ComponentWeights, const TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
  
  itkSetMacro(UsePrincipleComponents, bool);
  itkGetMacro(UsePrincipleComponents, bool);
  void SetUsePrincipleComponentsOn()
  {
    this->SetUsePrincipleComponents(true);
  }
  void SetUsePrincipleComponentsOff()
  {
    this->SetUsePrincipleComponents(false);
  }
  
protected:
  VectorGradientMagnitudeImageFilter();
  virtual ~VectorGradientMagnitudeImageFilter() {}

  void BeforeThreadedGenerateData ();

  void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
                            int threadId );

  void PrintSelf(std::ostream& os, Indent indent) const;

  TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
  {
    unsigned i, j;
    TRealType dx, sum, accum;
    accum = NumericTraits<TRealType>::Zero;
    for (i = 0; i < ImageDimension; ++i)
      {
        sum = NumericTraits<TRealType>::Zero;
        for (j = 0; j < VectorDimension; ++j)
          {
            dx = m_DerivativeWeights[i]
              * (0.5 * it.GetNext(i)[j] - 0.5 * it.GetPrevious(i)[j]);
            sum += m_ComponentWeights[j] * (dx * dx);
          }
        accum += sum;
      }
    return ::vnl_math_sqrt(accum);
  }

  // Function is defined here because the templating confuses gcc 2.96 when defined
  // in .txx file.  jc 1/29/03
  TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
  {
    unsigned int i, j;
    vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
    vnl_vector_fixed<TRealType, VectorDimension>
      d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
    
    // Calculate the directional derivatives for each vector component using
    // central differences.
    for (i = 0; i < ImageDimension; i++)
      {
        for (j = 0; j < VectorDimension; j++)
          {  d_phi_du[i][j] = m_DerivativeWeights[i]
               * (it.GetNext(i)[j] * 0.5 - it.GetPrevious(i)[j] * 0.5); }
      }
    
    // Calculate the symmetric metric tensor g
    for (i = 0; i < ImageDimension; i++)
      {
        for (j = i; j < ImageDimension; j++)
          {
            g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
          }
      }
    
    // Find the eigenvalues of g
    vnl_symmetric_eigensystem<TRealType> E(g);
    
    // Return the difference in length between the first two principle axes.
    // Note that other edge strength metrics may be appropriate here instead..
    return ( E.get_eigenvalue(ImageDimension - 1)
             - E.get_eigenvalue(ImageDimension - 2) ); 
  }
  
  TRealType m_DerivativeWeights[itk::GetImageDimension<TInputImage>::ImageDimension];

  TRealType m_ComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
  
private:
  bool m_UseImageSpacing;
  bool m_UsePrincipleComponents;
  int m_RequestedNumberOfThreads;

  typedef typename InputImageType::Superclass ImageBaseType;

  typename ImageBaseType::ConstPointer m_RealValuedInputImage;
  
  VectorGradientMagnitudeImageFilter(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented

  typename ConstNeighborhoodIteratorType::RadiusType m_NeighborhoodRadius;  
};
  
} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkVectorGradientMagnitudeImageFilter.txx"
#endif

#endif

---