itkVectorGradientMagnitudeImageFilter.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkVectorGradientMagnitudeImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2005/03/30 15:13:39 $
  Version:   $Revision: 1.11 $

  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 __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/vnl_vector_fixed.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 ConstNeighborhoodIteratorType::RadiusType RadiusType;  
  
  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);
  }

  static int CubicSolver(double *, double *);

protected:
  VectorGradientMagnitudeImageFilter();
  virtual ~VectorGradientMagnitudeImageFilter() {}

  void BeforeThreadedGenerateData ();

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

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

  typedef typename InputImageType::Superclass ImageBaseType;

  itkGetConstObjectMacro( RealValuedInputImage, ImageBaseType );

  itkGetConstReferenceMacro( NeighborhoodRadius, RadiusType );
  itkSetMacro( NeighborhoodRadius, RadiusType );
  

  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] * m_SqrtComponentWeights[j] 
          * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
        sum += dx * dx;
        }
      accum += sum;
      }
    return vcl_sqrt(accum);
  }

  TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType &it) const
  {
    // WARNING:  ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
    unsigned int i, j;
    double Lambda[3];
    double CharEqn[3];
    double ans;

    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] * m_SqrtComponentWeights[j]
             * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]); }
      }

    // 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 coefficients of the characteristic equation det(g - lambda I)=0
    //    CharEqn[3] = 1.0;

    CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );

    CharEqn[1] =(g[0][0]*g[1][1] + g[0][0]*g[2][2] + g[1][1]*g[2][2])
      - (g[0][1]*g[1][0] + g[0][2]*g[2][0] + g[1][2]*g[2][1]);

    CharEqn[0] = g[0][0] * ( g[1][2]*g[2][1] -  g[1][1]*g[2][2]  ) +
      g[1][0] * (  g[2][2]*g[0][1] - g[0][2]*g[2][1] ) +
      g[2][0] * ( g[1][1]*g[0][2] - g[0][1]*g[1][2] );

    //(g[0][0]*g[1][2]*g[2][1] + g[1][1]*g[0][2]*g[2][0] + g[2][2]*g[0][1]*g[1][0])
    //      - (g[0][0]*g[1][1]*g[2][2] + g[0][1]*g[2][0]*g[1][2] + g[0][2]*g[1][0]*g[2][1]);
    
    // Find the eigenvalues of g
    int numberOfDistinctRoots =  this->CubicSolver(CharEqn, Lambda);

    // Define gradient magnitude as the difference between two largest
    // eigenvalues.  Other definitions may be appropriate here as well.
    if (numberOfDistinctRoots == 3) // By far the most common case
      {
      if (Lambda[0] > Lambda[1])
        {
        if ( Lambda[1] > Lambda[2] )
          {  ans = Lambda[0] - Lambda[1]; } // Most common, guaranteed?
        else
          {
          if ( Lambda[0] > Lambda[2] )
            { ans = Lambda[0] - Lambda[2]; }
          else
            { ans = Lambda[2] - Lambda[0]; }
          }
        }
      else
        {
        if ( Lambda[0] > Lambda[2] )
          { ans = Lambda[1] - Lambda[0]; }
        else
          {
          if ( Lambda[1] > Lambda[2] )
            { ans = Lambda[1] - Lambda[2]; }
          else
            { ans = Lambda[2] - Lambda[1]; }
          }            
        }
      }
    else if (numberOfDistinctRoots == 2)
      {
      if ( Lambda[0] > Lambda[1] )
        { ans = Lambda[0] - Lambda[1]; }
      else
        { ans = Lambda[1] - Lambda[0]; }
      }
    else if (numberOfDistinctRoots == 1)
      {
      ans = 0.0;
      }
    else
      {
      itkExceptionMacro( << "Undefined condition. Cubic root solver returned "
                         << numberOfDistinctRoots << " distinct roots." );
      }

    return ans;  
  }
  
  // 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] * m_SqrtComponentWeights[j]
             * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j] ); }
      }
    
    // 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];
  TRealType m_SqrtComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];

private:
  bool m_UseImageSpacing;
  bool m_UsePrincipleComponents;
  int m_RequestedNumberOfThreads;


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

  RadiusType    m_NeighborhoodRadius;  
};
  
} // end namespace itk

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

#endif

---