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: 2008-10-17 16:30:55 $
00007   Version:   $Revision: 1.13 $
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/vnl_vector_fixed.h"
00027 #include "vnl/algo/vnl_symmetric_eigensystem.h"
00028 #include "vnl/vnl_math.h"
00029 
00030 namespace itk
00031 {
00133 template < typename TInputImage,
00134            typename TRealType = float,
00135            typename TOutputImage = Image< TRealType,
00136                                           ::itk::GetImageDimension<TInputImage>::ImageDimension >
00137 >
00138 class ITK_EXPORT VectorGradientMagnitudeImageFilter :
00139     public ImageToImageFilter< TInputImage, TOutputImage >
00140 {
00141 public:
00143   typedef VectorGradientMagnitudeImageFilter              Self;
00144   typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
00145   typedef SmartPointer<Self>                              Pointer;
00146   typedef SmartPointer<const Self>                        ConstPointer;
00147 
00149   itkNewMacro(Self);
00150 
00152   itkTypeMacro(VectorGradientMagnitudeImageFilter, ImageToImageFilter);
00153 
00156   typedef typename TOutputImage::PixelType OutputPixelType;
00157   typedef typename TInputImage::PixelType  InputPixelType;
00158 
00160   typedef TInputImage                       InputImageType;
00161   typedef TOutputImage                      OutputImageType;
00162   typedef typename InputImageType::Pointer  InputImagePointer;
00163   typedef typename OutputImageType::Pointer OutputImagePointer;
00164 
00166   itkStaticConstMacro(ImageDimension, unsigned int,
00167                       TOutputImage::ImageDimension);
00168 
00170   itkStaticConstMacro(VectorDimension, unsigned int,
00171                       InputPixelType::Dimension);
00172 
00174   typedef TRealType RealType;
00175   typedef Vector<TRealType, ::itk::GetVectorDimension<InputPixelType>::VectorDimension> RealVectorType;
00176   typedef Image<RealVectorType, ::itk::GetImageDimension<TInputImage>::ImageDimension>  RealVectorImageType;
00177 
00178 
00181   typedef ConstNeighborhoodIterator<RealVectorImageType> ConstNeighborhoodIteratorType;
00182   typedef typename ConstNeighborhoodIteratorType::RadiusType RadiusType;  
00183 
00185   typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
00186 
00195   virtual void GenerateInputRequestedRegion() throw(InvalidRequestedRegionError);
00196 
00200   void SetUseImageSpacingOn()
00201     { this->SetUseImageSpacing(true); }
00202 
00206   void SetUseImageSpacingOff()
00207     { this->SetUseImageSpacing(false); }
00208 
00211   void SetUseImageSpacing(bool);
00212   itkGetMacro(UseImageSpacing, bool);
00214 
00217   void SetDerivativeWeights(TRealType data[]);
00218   itkGetVectorMacro(DerivativeWeights, const TRealType, itk::GetImageDimension<TInputImage>::ImageDimension);
00220 
00223   itkSetVectorMacro(ComponentWeights, TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
00224   itkGetVectorMacro(ComponentWeights, const TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
00226 
00232   itkSetMacro(UsePrincipleComponents, bool);
00233   itkGetMacro(UsePrincipleComponents, bool);
00234   void SetUsePrincipleComponentsOn()
00235     {
00236     this->SetUsePrincipleComponents(true);
00237     }
00238   void SetUsePrincipleComponentsOff()
00239     {
00240     this->SetUsePrincipleComponents(false);
00241     }
00243 
00246   static int CubicSolver(double *, double *);
00247 
00248 #ifdef ITK_USE_CONCEPT_CHECKING
00249 
00250   itkConceptMacro(InputHasNumericTraitsCheck,
00251                   (Concept::HasNumericTraits<typename InputPixelType::ValueType>));
00252   itkConceptMacro(RealTypeHasNumericTraitsCheck,
00253                   (Concept::HasNumericTraits<RealType>));
00254 
00256 #endif
00257 
00258 protected:
00259   VectorGradientMagnitudeImageFilter();
00260   virtual ~VectorGradientMagnitudeImageFilter() {}
00261 
00265   void BeforeThreadedGenerateData ();
00266 
00278   void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
00279                             int threadId );
00280 
00281   void PrintSelf(std::ostream& os, Indent indent) const;
00282 
00283   typedef typename InputImageType::Superclass ImageBaseType;
00284 
00286   itkGetConstObjectMacro( RealValuedInputImage, ImageBaseType );
00287 
00289   itkGetConstReferenceMacro( NeighborhoodRadius, RadiusType );
00290   itkSetMacro( NeighborhoodRadius, RadiusType );
00292 
00293 
00294   TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
00295   {
00296     unsigned i, j;
00297     TRealType dx, sum, accum;
00298     accum = NumericTraits<TRealType>::Zero;
00299     for (i = 0; i < ImageDimension; ++i)
00300       {
00301       sum = NumericTraits<TRealType>::Zero;
00302       for (j = 0; j < VectorDimension; ++j)
00303         {
00304         dx =  m_DerivativeWeights[i] * m_SqrtComponentWeights[j] 
00305           * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
00306         sum += dx * dx;
00307         }
00308       accum += sum;
00309       }
00310     return vcl_sqrt(accum);
00311   }
00312 
00313   TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType &it) const
00314   {
00315     // WARNING:  ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
00316     unsigned int i, j;
00317     double Lambda[3];
00318     double CharEqn[3];
00319     double ans;
00320 
00321     vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
00322     vnl_vector_fixed<TRealType, VectorDimension>
00323       d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
00324 
00325     // Calculate the directional derivatives for each vector component using
00326     // central differences.
00327     for (i = 0; i < ImageDimension; i++)
00328       {
00329       for (j = 0; j < VectorDimension; j++)
00330         {  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00331              * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]); }
00332       }
00333 
00334     // Calculate the symmetric metric tensor g
00335     for (i = 0; i < ImageDimension; i++)
00336       {
00337       for (j = i; j < ImageDimension; j++)
00338         {
00339         g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00340         }
00341       }
00342 
00343     // Find the coefficients of the characteristic equation det(g - lambda I)=0
00344     //    CharEqn[3] = 1.0;
00345 
00346     CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
00347 
00348     CharEqn[1] =(g[0][0]*g[1][1] + g[0][0]*g[2][2] + g[1][1]*g[2][2])
00349       - (g[0][1]*g[1][0] + g[0][2]*g[2][0] + g[1][2]*g[2][1]);
00350 
00351     CharEqn[0] = g[0][0] * ( g[1][2]*g[2][1] -  g[1][1]*g[2][2]  ) +
00352       g[1][0] * (  g[2][2]*g[0][1] - g[0][2]*g[2][1] ) +
00353       g[2][0] * ( g[1][1]*g[0][2] - g[0][1]*g[1][2] );
00354 
00355     // Find the eigenvalues of g
00356     int numberOfDistinctRoots =  this->CubicSolver(CharEqn, Lambda);
00357 
00358     // Define gradient magnitude as the difference between two largest
00359     // eigenvalues.  Other definitions may be appropriate here as well.
00360     if (numberOfDistinctRoots == 3) // By far the most common case
00361       {
00362       if (Lambda[0] > Lambda[1])
00363         {
00364         if ( Lambda[1] > Lambda[2] )  // Most common, guaranteed?
00365           {
00366           ans = Lambda[0] - Lambda[1];
00367           }
00368         else
00369           {
00370           if ( Lambda[0] > Lambda[2] )
00371             {
00372             ans = Lambda[0] - Lambda[2];
00373             }
00374           else
00375             {
00376             ans = Lambda[2] - Lambda[0];
00377             }
00378           }
00379         }
00380       else
00381         {
00382         if ( Lambda[0] > Lambda[2] )
00383           {
00384           ans = Lambda[1] - Lambda[0];
00385           }
00386         else
00387           {
00388           if ( Lambda[1] > Lambda[2] )
00389             {
00390             ans = Lambda[1] - Lambda[2];
00391             }
00392           else
00393             {
00394             ans = Lambda[2] - Lambda[1];
00395             }
00396           }
00397         }
00398       }
00399     else if (numberOfDistinctRoots == 2)
00400       {
00401       if ( Lambda[0] > Lambda[1] )
00402         {
00403         ans = Lambda[0] - Lambda[1];
00404         }
00405       else
00406         {
00407         ans = Lambda[1] - Lambda[0];
00408         }
00409       }
00410     else if (numberOfDistinctRoots == 1)
00411       {
00412       ans = 0.0;
00413       }
00414     else
00415       {
00416       itkExceptionMacro( << "Undefined condition. Cubic root solver returned "
00417                          << numberOfDistinctRoots << " distinct roots." );
00418       }
00419 
00420     return ans;  
00421   }
00422   
00423   // Function is defined here because the templating confuses gcc 2.96 when defined
00424   // in .txx file.  jc 1/29/03
00425   TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
00426     {
00427     unsigned int i, j;
00428     vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
00429     vnl_vector_fixed<TRealType, VectorDimension>
00430       d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
00431     
00432     // Calculate the directional derivatives for each vector component using
00433     // central differences.
00434     for (i = 0; i < ImageDimension; i++)
00435       {
00436       for (j = 0; j < VectorDimension; j++)
00437         {
00438         d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00439           * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j] );
00440         }
00441       }
00442     
00443     // Calculate the symmetric metric tensor g
00444     for (i = 0; i < ImageDimension; i++)
00445       {
00446       for (j = i; j < ImageDimension; j++)
00447         {
00448         g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00449         }
00450       }
00451     
00452     // Find the eigenvalues of g
00453     vnl_symmetric_eigensystem<TRealType> E(g);
00454 
00455     // Return the difference in length between the first two principle axes.
00456     // Note that other edge strength metrics may be appropriate here instead..
00457     return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
00458     }
00459   
00461   TRealType m_DerivativeWeights[itk::GetImageDimension<TInputImage>::ImageDimension];
00462 
00466   TRealType m_ComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00467   TRealType m_SqrtComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00468 
00469 private:
00470   bool m_UseImageSpacing;
00471   bool m_UsePrincipleComponents;
00472   int  m_RequestedNumberOfThreads;
00473 
00474   typename ImageBaseType::ConstPointer m_RealValuedInputImage;
00475   
00476   VectorGradientMagnitudeImageFilter(const Self&); //purposely not implemented
00477   void operator=(const Self&); //purposely not implemented
00478 
00479   RadiusType    m_NeighborhoodRadius;  
00480 };
00481   
00482 } // end namespace itk
00483 
00484 #ifndef ITK_MANUAL_INSTANTIATION
00485 #include "itkVectorGradientMagnitudeImageFilter.txx"
00486 #endif
00487 
00488 #endif
00489 

Generated at Sat Feb 28 13:45:26 2009 for ITK by doxygen 1.5.6 written by Dimitri van Heesch, © 1997-2000