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: 2006/04/03 15:07:52 $
00007   Version:   $Revision: 1.12 $
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     //(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])
00356     //      - (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]);
00357     
00358     // Find the eigenvalues of g
00359     int numberOfDistinctRoots =  this->CubicSolver(CharEqn, Lambda);
00360 
00361     // Define gradient magnitude as the difference between two largest
00362     // eigenvalues.  Other definitions may be appropriate here as well.
00363     if (numberOfDistinctRoots == 3) // By far the most common case
00364       {
00365       if (Lambda[0] > Lambda[1])
00366         {
00367         if ( Lambda[1] > Lambda[2] )
00368           {  ans = Lambda[0] - Lambda[1]; } // Most common, guaranteed?
00369         else
00370           {
00371           if ( Lambda[0] > Lambda[2] )
00372             { ans = Lambda[0] - Lambda[2]; }
00373           else
00374             { ans = Lambda[2] - Lambda[0]; }
00375           }
00376         }
00377       else
00378         {
00379         if ( Lambda[0] > Lambda[2] )
00380           { ans = Lambda[1] - Lambda[0]; }
00381         else
00382           {
00383           if ( Lambda[1] > Lambda[2] )
00384             { ans = Lambda[1] - Lambda[2]; }
00385           else
00386             { ans = Lambda[2] - Lambda[1]; }
00387           }            
00388         }
00389       }
00390     else if (numberOfDistinctRoots == 2)
00391       {
00392       if ( Lambda[0] > Lambda[1] )
00393         { ans = Lambda[0] - Lambda[1]; }
00394       else
00395         { ans = Lambda[1] - Lambda[0]; }
00396       }
00397     else if (numberOfDistinctRoots == 1)
00398       {
00399       ans = 0.0;
00400       }
00401     else
00402       {
00403       itkExceptionMacro( << "Undefined condition. Cubic root solver returned "
00404                          << numberOfDistinctRoots << " distinct roots." );
00405       }
00406 
00407     return ans;  
00408   }
00409   
00410   // Function is defined here because the templating confuses gcc 2.96 when defined
00411   // in .txx file.  jc 1/29/03
00412   TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
00413   {
00414     unsigned int i, j;
00415     vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
00416     vnl_vector_fixed<TRealType, VectorDimension>
00417       d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
00418     
00419     // Calculate the directional derivatives for each vector component using
00420     // central differences.
00421     for (i = 0; i < ImageDimension; i++)
00422       {
00423       for (j = 0; j < VectorDimension; j++)
00424         {  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00425              * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j] ); }
00426       }
00427     
00428     // Calculate the symmetric metric tensor g
00429     for (i = 0; i < ImageDimension; i++)
00430       {
00431       for (j = i; j < ImageDimension; j++)
00432         {
00433         g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00434         }
00435       }
00436     
00437     // Find the eigenvalues of g
00438     vnl_symmetric_eigensystem<TRealType> E(g);
00439 
00440     // Return the difference in length between the first two principle axes.
00441     // Note that other edge strength metrics may be appropriate here instead..
00442     return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
00443   }
00444   
00446   TRealType m_DerivativeWeights[itk::GetImageDimension<TInputImage>::ImageDimension];
00447 
00451   TRealType m_ComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00452   TRealType m_SqrtComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00453 
00454 private:
00455   bool m_UseImageSpacing;
00456   bool m_UsePrincipleComponents;
00457   int m_RequestedNumberOfThreads;
00458 
00459 
00460   typename ImageBaseType::ConstPointer m_RealValuedInputImage;
00461   
00462   VectorGradientMagnitudeImageFilter(const Self&); //purposely not implemented
00463   void operator=(const Self&); //purposely not implemented
00464 
00465   RadiusType    m_NeighborhoodRadius;  
00466 };
00467   
00468 } // end namespace itk
00469 
00470 #ifndef ITK_MANUAL_INSTANTIATION
00471 #include "itkVectorGradientMagnitudeImageFilter.txx"
00472 #endif
00473 
00474 #endif
00475 

Generated at Mon Mar 12 03:19:11 2007 for ITK by doxygen 1.5.1 written by Dimitri van Heesch, © 1997-2000