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: 2004/04/15 15:42:47 $
00007   Version:   $Revision: 1.10 $
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   typedef typename ConstNeighborhoodIteratorType::RadiusType RadiusType;  
00182   
00184   typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
00185 
00194   virtual void GenerateInputRequestedRegion() throw(InvalidRequestedRegionError);
00195 
00199   void SetUseImageSpacingOn()
00200   { this->SetUseImageSpacing(true); }
00201 
00205   void SetUseImageSpacingOff()
00206   { this->SetUseImageSpacing(false); }
00207 
00210   void SetUseImageSpacing(bool);                         
00211   itkGetMacro(UseImageSpacing, bool);
00212 
00215   void SetDerivativeWeights(TRealType data[]);
00216   itkGetVectorMacro(DerivativeWeights, const TRealType, itk::GetImageDimension<TInputImage>::ImageDimension);
00217 
00220   itkSetVectorMacro(ComponentWeights, TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
00221   itkGetVectorMacro(ComponentWeights, const TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
00222   
00228   itkSetMacro(UsePrincipleComponents, bool);
00229   itkGetMacro(UsePrincipleComponents, bool);
00230   void SetUsePrincipleComponentsOn()
00231   {
00232     this->SetUsePrincipleComponents(true);
00233   }
00234   void SetUsePrincipleComponentsOff()
00235   {
00236     this->SetUsePrincipleComponents(false);
00237   }
00238 
00241   static int CubicSolver(double *, double *);
00242 
00243 protected:
00244   VectorGradientMagnitudeImageFilter();
00245   virtual ~VectorGradientMagnitudeImageFilter() {}
00246 
00250   void BeforeThreadedGenerateData ();
00251 
00263   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
00264                             int threadId );
00265 
00266   void PrintSelf(std::ostream& os, Indent indent) const;
00267 
00268   typedef typename InputImageType::Superclass ImageBaseType;
00269 
00271   itkGetConstObjectMacro( RealValuedInputImage, ImageBaseType );
00272 
00274   itkGetConstReferenceMacro( NeighborhoodRadius, RadiusType );
00275   itkSetMacro( NeighborhoodRadius, RadiusType );
00276   
00277 
00278   TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
00279   {
00280     unsigned i, j;
00281     TRealType dx, sum, accum;
00282     accum = NumericTraits<TRealType>::Zero;
00283     for (i = 0; i < ImageDimension; ++i)
00284       {
00285       sum = NumericTraits<TRealType>::Zero;
00286       for (j = 0; j < VectorDimension; ++j)
00287         {
00288         dx =  m_DerivativeWeights[i] * m_SqrtComponentWeights[j] 
00289           * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
00290         sum += dx * dx;
00291         }
00292       accum += sum;
00293       }
00294     return vcl_sqrt(accum);
00295   }
00296 
00297   TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType &it) const
00298   {
00299     // WARNING:  ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
00300     unsigned int i, j;
00301     double Lambda[3];
00302     double CharEqn[3];
00303     double ans;
00304 
00305     vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
00306     vnl_vector_fixed<TRealType, VectorDimension>
00307       d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
00308 
00309     // Calculate the directional derivatives for each vector component using
00310     // central differences.
00311     for (i = 0; i < ImageDimension; i++)
00312       {
00313       for (j = 0; j < VectorDimension; j++)
00314         {  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00315              * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]); }
00316       }
00317 
00318     // Calculate the symmetric metric tensor g
00319     for (i = 0; i < ImageDimension; i++)
00320       {
00321       for (j = i; j < ImageDimension; j++)
00322         {
00323         g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00324         }
00325       }
00326 
00327     // Find the coefficients of the characteristic equation det(g - lambda I)=0
00328     //    CharEqn[3] = 1.0;
00329 
00330     CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
00331 
00332     CharEqn[1] =(g[0][0]*g[1][1] + g[0][0]*g[2][2] + g[1][1]*g[2][2])
00333       - (g[0][1]*g[1][0] + g[0][2]*g[2][0] + g[1][2]*g[2][1]);
00334 
00335     CharEqn[0] = g[0][0] * ( g[1][2]*g[2][1] -  g[1][1]*g[2][2]  ) +
00336       g[1][0] * (  g[2][2]*g[0][1] - g[0][2]*g[2][1] ) +
00337       g[2][0] * ( g[1][1]*g[0][2] - g[0][1]*g[1][2] );
00338 
00339     //(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])
00340     //      - (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]);
00341     
00342     // Find the eigenvalues of g
00343     int numberOfDistinctRoots =  this->CubicSolver(CharEqn, Lambda);
00344 
00345     // Define gradient magnitude as the difference between two largest
00346     // eigenvalues.  Other definitions may be appropriate here as well.
00347     if (numberOfDistinctRoots == 3) // By far the most common case
00348       {
00349       if (Lambda[0] > Lambda[1])
00350         {
00351         if ( Lambda[1] > Lambda[2] )
00352           {  ans = Lambda[0] - Lambda[1]; } // Most common, guaranteed?
00353         else
00354           {
00355           if ( Lambda[0] > Lambda[2] )
00356             { ans = Lambda[0] - Lambda[2]; }
00357           else
00358             { ans = Lambda[2] - Lambda[0]; }
00359           }
00360         }
00361       else
00362         {
00363         if ( Lambda[0] > Lambda[2] )
00364           { ans = Lambda[1] - Lambda[0]; }
00365         else
00366           {
00367           if ( Lambda[1] > Lambda[2] )
00368             { ans = Lambda[1] - Lambda[2]; }
00369           else
00370             { ans = Lambda[2] - Lambda[1]; }
00371           }            
00372         }
00373       }
00374     else if (numberOfDistinctRoots == 2)
00375       {
00376       if ( Lambda[0] > Lambda[1] )
00377         { ans = Lambda[0] - Lambda[1]; }
00378       else
00379         { ans = Lambda[1] - Lambda[0]; }
00380       }
00381     else if (numberOfDistinctRoots == 1)
00382       {
00383       ans = 0.0;
00384       }
00385     else
00386       {
00387       itkExceptionMacro( << "Undefined condition. Cubic root solver returned "
00388                          << numberOfDistinctRoots << " distinct roots." );
00389       }
00390 
00391     return ans;  
00392   }
00393   
00394   // Function is defined here because the templating confuses gcc 2.96 when defined
00395   // in .txx file.  jc 1/29/03
00396   TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
00397   {
00398     unsigned int i, j;
00399     vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
00400     vnl_vector_fixed<TRealType, VectorDimension>
00401       d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
00402     
00403     // Calculate the directional derivatives for each vector component using
00404     // central differences.
00405     for (i = 0; i < ImageDimension; i++)
00406       {
00407       for (j = 0; j < VectorDimension; j++)
00408         {  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00409              * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j] ); }
00410       }
00411     
00412     // Calculate the symmetric metric tensor g
00413     for (i = 0; i < ImageDimension; i++)
00414       {
00415       for (j = i; j < ImageDimension; j++)
00416         {
00417         g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00418         }
00419       }
00420     
00421     // Find the eigenvalues of g
00422     vnl_symmetric_eigensystem<TRealType> E(g);
00423 
00424     // Return the difference in length between the first two principle axes.
00425     // Note that other edge strength metrics may be appropriate here instead..
00426     return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
00427   }
00428   
00430   TRealType m_DerivativeWeights[itk::GetImageDimension<TInputImage>::ImageDimension];
00431 
00435   TRealType m_ComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00436   TRealType m_SqrtComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00437 
00438 private:
00439   bool m_UseImageSpacing;
00440   bool m_UsePrincipleComponents;
00441   int m_RequestedNumberOfThreads;
00442 
00443 
00444   typename ImageBaseType::ConstPointer m_RealValuedInputImage;
00445   
00446   VectorGradientMagnitudeImageFilter(const Self&); //purposely not implemented
00447   void operator=(const Self&); //purposely not implemented
00448 
00449   RadiusType    m_NeighborhoodRadius;  
00450 };
00451   
00452 } // end namespace itk
00453 
00454 #ifndef ITK_MANUAL_INSTANTIATION
00455 #include "itkVectorGradientMagnitudeImageFilter.txx"
00456 #endif
00457 
00458 #endif

Generated at Wed Mar 30 00:13:58 2005 for ITK by doxygen 1.3.9.1 written by Dimitri van Heesch, © 1997-2000