Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkVectorGradientMagnitudeImageFilter.h

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkVectorGradientMagnitudeImageFilter.h,v $
00005   Language:  C++
00006   Date:      $Date: 2003/09/10 14:28:59 $
00007   Version:   $Revision: 1.9 $
00008 
00009   Copyright (c) Insight Software Consortium. All rights reserved.
00010   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
00011 
00012      This software is distributed WITHOUT ANY WARRANTY; without even 
00013      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
00014      PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 #ifndef __itkVectorGradientMagnitudeImageFilter_h
00018 #define __itkVectorGradientMagnitudeImageFilter_h
00019 
00020 #include "itkConstNeighborhoodIterator.h"
00021 #include "itkNeighborhoodIterator.h"
00022 #include "itkImageToImageFilter.h"
00023 #include "itkImage.h"
00024 #include "itkVector.h"
00025 #include "vnl/vnl_matrix.h"
00026 #include "vnl/algo/vnl_symmetric_eigensystem.h"
00027 #include "vnl/vnl_math.h"
00028 
00029 namespace itk
00030 {
00132 template < typename TInputImage,
00133            typename TRealType = float,
00134            typename TOutputImage = Image< TRealType,
00135                                           ::itk::GetImageDimension<TInputImage>::ImageDimension >
00136 >
00137 class ITK_EXPORT VectorGradientMagnitudeImageFilter :
00138     public ImageToImageFilter< TInputImage, TOutputImage >
00139 {
00140 public:
00142   typedef VectorGradientMagnitudeImageFilter Self;
00143   typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
00144   typedef SmartPointer<Self> Pointer;
00145   typedef SmartPointer<const Self>  ConstPointer;
00146 
00148   itkNewMacro(Self);
00149 
00151   itkTypeMacro(VectorGradientMagnitudeImageFilter, ImageToImageFilter);
00152   
00155   typedef typename TOutputImage::PixelType OutputPixelType;
00156   typedef typename TInputImage::PixelType  InputPixelType;
00157 
00159   typedef TInputImage  InputImageType;
00160   typedef TOutputImage OutputImageType;
00161   typedef typename InputImageType::Pointer  InputImagePointer;
00162   typedef typename OutputImageType::Pointer OutputImagePointer;
00163   
00165   itkStaticConstMacro(ImageDimension, unsigned int,
00166                       TOutputImage::ImageDimension);
00167   
00169   itkStaticConstMacro(VectorDimension, unsigned int,
00170                       InputPixelType::Dimension);
00171 
00173   typedef TRealType RealType;
00174   typedef Vector<TRealType, ::itk::GetVectorDimension<InputPixelType>::VectorDimension> RealVectorType;
00175   typedef Image<RealVectorType, ::itk::GetImageDimension<TInputImage>::ImageDimension>  RealVectorImageType;
00176   
00177 
00180   typedef ConstNeighborhoodIterator<RealVectorImageType> ConstNeighborhoodIteratorType;
00181   
00183   typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
00184 
00193   virtual void GenerateInputRequestedRegion() throw(InvalidRequestedRegionError);
00194 
00198   void SetUseImageSpacingOn()
00199   { this->SetUseImageSpacing(true); }
00200 
00204   void SetUseImageSpacingOff()
00205   { this->SetUseImageSpacing(false); }
00206 
00209   void SetUseImageSpacing(bool);                         
00210   itkGetMacro(UseImageSpacing, bool);
00211 
00214   void SetDerivativeWeights(TRealType data[]);
00215   itkGetVectorMacro(DerivativeWeights, const TRealType, itk::GetImageDimension<TInputImage>::ImageDimension);
00216 
00219   itkSetVectorMacro(ComponentWeights, TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
00220   itkGetVectorMacro(ComponentWeights, const TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
00221   
00227   itkSetMacro(UsePrincipleComponents, bool);
00228   itkGetMacro(UsePrincipleComponents, bool);
00229   void SetUsePrincipleComponentsOn()
00230   {
00231     this->SetUsePrincipleComponents(true);
00232   }
00233   void SetUsePrincipleComponentsOff()
00234   {
00235     this->SetUsePrincipleComponents(false);
00236   }
00237 
00240   static int CubicSolver(double *, double *);
00241 
00242 protected:
00243   VectorGradientMagnitudeImageFilter();
00244   virtual ~VectorGradientMagnitudeImageFilter() {}
00245 
00249   void BeforeThreadedGenerateData ();
00250 
00262   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
00263                             int threadId );
00264 
00265   void PrintSelf(std::ostream& os, Indent indent) const;
00266 
00267   TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
00268   {
00269     unsigned i, j;
00270     TRealType dx, sum, accum;
00271     accum = NumericTraits<TRealType>::Zero;
00272     for (i = 0; i < ImageDimension; ++i)
00273       {
00274       sum = NumericTraits<TRealType>::Zero;
00275       for (j = 0; j < VectorDimension; ++j)
00276         {
00277         dx =  m_DerivativeWeights[i] * m_SqrtComponentWeights[j] 
00278           * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
00279         sum += dx * dx;
00280         }
00281       accum += sum;
00282       }
00283     return vcl_sqrt(accum);
00284   }
00285 
00286   TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType &it) const
00287   {
00288     // WARNING:  ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
00289     unsigned int i, j;
00290     double Lambda[3];
00291     double CharEqn[3];
00292     double ans;
00293 
00294     vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
00295     vnl_vector_fixed<TRealType, VectorDimension>
00296       d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
00297 
00298     // Calculate the directional derivatives for each vector component using
00299     // central differences.
00300     for (i = 0; i < ImageDimension; i++)
00301       {
00302       for (j = 0; j < VectorDimension; j++)
00303         {  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00304              * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]); }
00305       }
00306 
00307     // Calculate the symmetric metric tensor g
00308     for (i = 0; i < ImageDimension; i++)
00309       {
00310       for (j = i; j < ImageDimension; j++)
00311         {
00312         g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00313         }
00314       }
00315 
00316     // Find the coefficients of the characteristic equation det(g - lambda I)=0
00317     //    CharEqn[3] = 1.0;
00318 
00319     CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
00320 
00321     CharEqn[1] =(g[0][0]*g[1][1] + g[0][0]*g[2][2] + g[1][1]*g[2][2])
00322       - (g[0][1]*g[1][0] + g[0][2]*g[2][0] + g[1][2]*g[2][1]);
00323 
00324     CharEqn[0] = g[0][0] * ( g[1][2]*g[2][1] -  g[1][1]*g[2][2]  ) +
00325       g[1][0] * (  g[2][2]*g[0][1] - g[0][2]*g[2][1] ) +
00326       g[2][0] * ( g[1][1]*g[0][2] - g[0][1]*g[1][2] );
00327 
00328     //(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])
00329     //      - (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]);
00330     
00331     // Find the eigenvalues of g
00332     int numberOfDistinctRoots =  this->CubicSolver(CharEqn, Lambda);
00333 
00334     // Define gradient magnitude as the difference between two largest
00335     // eigenvalues.  Other definitions may be appropriate here as well.
00336     if (numberOfDistinctRoots == 3) // By far the most common case
00337       {
00338       if (Lambda[0] > Lambda[1])
00339         {
00340         if ( Lambda[1] > Lambda[2] )
00341           {  ans = Lambda[0] - Lambda[1]; } // Most common, guaranteed?
00342         else
00343           {
00344           if ( Lambda[0] > Lambda[2] )
00345             { ans = Lambda[0] - Lambda[2]; }
00346           else
00347             { ans = Lambda[2] - Lambda[0]; }
00348           }
00349         }
00350       else
00351         {
00352         if ( Lambda[0] > Lambda[2] )
00353           { ans = Lambda[1] - Lambda[0]; }
00354         else
00355           {
00356           if ( Lambda[1] > Lambda[2] )
00357             { ans = Lambda[1] - Lambda[2]; }
00358           else
00359             { ans = Lambda[2] - Lambda[1]; }
00360           }            
00361         }
00362       }
00363     else if (numberOfDistinctRoots == 2)
00364       {
00365       if ( Lambda[0] > Lambda[1] )
00366         { ans = Lambda[0] - Lambda[1]; }
00367       else
00368         { ans = Lambda[1] - Lambda[0]; }
00369       }
00370     else if (numberOfDistinctRoots == 1)
00371       {
00372       ans = 0.0;
00373       }
00374     else
00375       {
00376       itkExceptionMacro( << "Undefined condition. Cubic root solver returned "
00377                          << numberOfDistinctRoots << " distinct roots." );
00378       }
00379 
00380     return ans;  
00381   }
00382   
00383   // Function is defined here because the templating confuses gcc 2.96 when defined
00384   // in .txx file.  jc 1/29/03
00385   TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
00386   {
00387     unsigned int i, j;
00388     vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
00389     vnl_vector_fixed<TRealType, VectorDimension>
00390       d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
00391     
00392     // Calculate the directional derivatives for each vector component using
00393     // central differences.
00394     for (i = 0; i < ImageDimension; i++)
00395       {
00396       for (j = 0; j < VectorDimension; j++)
00397         {  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00398              * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j] ); }
00399       }
00400     
00401     // Calculate the symmetric metric tensor g
00402     for (i = 0; i < ImageDimension; i++)
00403       {
00404       for (j = i; j < ImageDimension; j++)
00405         {
00406         g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00407         }
00408       }
00409     
00410     // Find the eigenvalues of g
00411     vnl_symmetric_eigensystem<TRealType> E(g);
00412 
00413     // Return the difference in length between the first two principle axes.
00414     // Note that other edge strength metrics may be appropriate here instead..
00415     return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
00416   }
00417   
00419   TRealType m_DerivativeWeights[itk::GetImageDimension<TInputImage>::ImageDimension];
00420 
00424   TRealType m_ComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00425   TRealType m_SqrtComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00426 
00427 private:
00428   bool m_UseImageSpacing;
00429   bool m_UsePrincipleComponents;
00430   int m_RequestedNumberOfThreads;
00431 
00432   typedef typename InputImageType::Superclass ImageBaseType;
00433 
00434   typename ImageBaseType::ConstPointer m_RealValuedInputImage;
00435   
00436   VectorGradientMagnitudeImageFilter(const Self&); //purposely not implemented
00437   void operator=(const Self&); //purposely not implemented
00438 
00439   typename ConstNeighborhoodIteratorType::RadiusType m_NeighborhoodRadius;  
00440 };
00441   
00442 } // end namespace itk
00443 
00444 #ifndef ITK_MANUAL_INSTANTIATION
00445 #include "itkVectorGradientMagnitudeImageFilter.txx"
00446 #endif
00447 
00448 #endif

Generated at Sun Jan 25 13:19:47 2004 for ITK by doxygen 1.3.3 written by Dimitri van Heesch, © 1997-2000