ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkVectorGradientMagnitudeImageFilter.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 #ifndef __itkVectorGradientMagnitudeImageFilter_h
00019 #define __itkVectorGradientMagnitudeImageFilter_h
00020 
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 {
00134 template< typename TInputImage,
00135           typename TRealType = float,
00136           typename TOutputImage = Image< TRealType,
00137                                          ::itk::GetImageDimension< TInputImage >::ImageDimension >
00138           >
00139 class ITK_EXPORT VectorGradientMagnitudeImageFilter:
00140   public ImageToImageFilter< TInputImage, TOutputImage >
00141 {
00142 public:
00144   typedef VectorGradientMagnitudeImageFilter              Self;
00145   typedef ImageToImageFilter< TInputImage, TOutputImage > Superclass;
00146   typedef SmartPointer< Self >                            Pointer;
00147   typedef SmartPointer< const Self >                      ConstPointer;
00148 
00150   itkNewMacro(Self);
00151 
00153   itkTypeMacro(VectorGradientMagnitudeImageFilter, ImageToImageFilter);
00154 
00157   typedef typename TOutputImage::PixelType OutputPixelType;
00158   typedef typename TInputImage::PixelType  InputPixelType;
00159 
00161   typedef TInputImage                       InputImageType;
00162   typedef TOutputImage                      OutputImageType;
00163   typedef typename InputImageType::Pointer  InputImagePointer;
00164   typedef typename OutputImageType::Pointer OutputImagePointer;
00165 
00167   itkStaticConstMacro(ImageDimension, unsigned int,
00168                       TOutputImage::ImageDimension);
00169 
00171   itkStaticConstMacro(VectorDimension, unsigned int,
00172                       InputPixelType::Dimension);
00173 
00175   typedef TRealType                                                                         RealType;
00176   typedef Vector< TRealType, ::itk::GetVectorDimension< InputPixelType >::VectorDimension > RealVectorType;
00177   typedef Image< RealVectorType, ::itk::GetImageDimension< TInputImage >::ImageDimension >  RealVectorImageType;
00178 
00181   typedef ConstNeighborhoodIterator< RealVectorImageType >   ConstNeighborhoodIteratorType;
00182   typedef typename ConstNeighborhoodIteratorType::RadiusType RadiusType;
00183 
00185   typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
00186 
00195   virtual void GenerateInputRequestedRegion()
00196   throw( InvalidRequestedRegionError );
00197 
00202   void SetUseImageSpacingOn()
00203   { this->SetUseImageSpacing(true); }
00204 
00208   void SetUseImageSpacingOff()
00209   { this->SetUseImageSpacing(false); }
00210 
00213   void SetUseImageSpacing(bool);
00214 
00215   itkGetConstMacro(UseImageSpacing, bool);
00216 
00217   typedef FixedArray< TRealType, VectorDimension > WeightsType;
00218 
00221   itkSetMacro(DerivativeWeights, WeightsType);
00222   itkGetConstReferenceMacro(DerivativeWeights, WeightsType);
00224 
00227   itkSetMacro(ComponentWeights, WeightsType);
00228   itkGetConstReferenceMacro(ComponentWeights, WeightsType);
00230 
00232   itkGetConstReferenceMacro(NeighborhoodRadius, RadiusType);
00233   itkSetMacro(NeighborhoodRadius, RadiusType);
00235 
00241   itkSetMacro(UsePrincipleComponents, bool);
00242   itkGetConstMacro(UsePrincipleComponents, bool);
00243   void SetUsePrincipleComponentsOn()
00244   {
00245     this->SetUsePrincipleComponents(true);
00246   }
00248 
00249   void SetUsePrincipleComponentsOff()
00250   {
00251     this->SetUsePrincipleComponents(false);
00252   }
00253 
00256   static int CubicSolver(double *, double *);
00257 
00258 #ifdef ITK_USE_CONCEPT_CHECKING
00259 
00260   itkConceptMacro( InputHasNumericTraitsCheck,
00261                    ( Concept::HasNumericTraits< typename InputPixelType::ValueType > ) );
00262   itkConceptMacro( RealTypeHasNumericTraitsCheck,
00263                    ( Concept::HasNumericTraits< RealType > ) );
00264 
00266 #endif
00267 protected:
00268   VectorGradientMagnitudeImageFilter();
00269   virtual ~VectorGradientMagnitudeImageFilter() {}
00270 
00274   void BeforeThreadedGenerateData();
00275 
00287   void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread,
00288                             ThreadIdType threadId);
00289 
00290   void PrintSelf(std::ostream & os, Indent indent) const;
00291 
00292   typedef typename InputImageType::Superclass ImageBaseType;
00293 
00295   itkGetConstObjectMacro(RealValuedInputImage, ImageBaseType);
00296 
00297   TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType & it) const
00298   {
00299     unsigned  i, j;
00300     TRealType dx, sum, accum;
00301 
00302     accum = NumericTraits< TRealType >::Zero;
00303     for ( i = 0; i < ImageDimension; ++i )
00304       {
00305       sum = NumericTraits< TRealType >::Zero;
00306       for ( j = 0; j < VectorDimension; ++j )
00307         {
00308         dx =  m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00309              * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
00310         sum += dx * dx;
00311         }
00312       accum += sum;
00313       }
00314     return vcl_sqrt(accum);
00315   }
00316 
00317   TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType & it) const
00318   {
00319     // WARNING:  ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
00320     unsigned int i, j;
00321     double       Lambda[3];
00322     double       CharEqn[3];
00323     double       ans;
00324 
00325     vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
00326     vnl_vector_fixed< TRealType, VectorDimension >
00327     d_phi_du[itk::GetImageDimension < TInputImage > ::ImageDimension];
00328 
00329     // Calculate the directional derivatives for each vector component using
00330     // central differences.
00331     for ( i = 0; i < ImageDimension; i++ )
00332       {
00333       for ( j = 0; j < VectorDimension; j++ )
00334         {
00335         d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00336                          * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
00337         }
00338       }
00339 
00340     // Calculate the symmetric metric tensor g
00341     for ( i = 0; i < ImageDimension; i++ )
00342       {
00343       for ( j = i; j < ImageDimension; j++ )
00344         {
00345         g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00346         }
00347       }
00348 
00349     // Find the coefficients of the characteristic equation det(g - lambda I)=0
00350     //    CharEqn[3] = 1.0;
00351 
00352     CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
00353 
00354     CharEqn[1] = ( g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2] )
00355                  - ( g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1] );
00356 
00357     CharEqn[0] = g[0][0] * ( g[1][2] * g[2][1] -  g[1][1] * g[2][2]  )
00358                  + g[1][0] * (  g[2][2] * g[0][1] - g[0][2] * g[2][1] )
00359                  + g[2][0] * ( g[1][1] * g[0][2] - g[0][1] * g[1][2] );
00360 
00361     // Find the eigenvalues of g
00362     int numberOfDistinctRoots =  this->CubicSolver(CharEqn, Lambda);
00363 
00364     // Define gradient magnitude as the difference between two largest
00365     // eigenvalues.  Other definitions may be appropriate here as well.
00366     if ( numberOfDistinctRoots == 3 ) // By far the most common case
00367       {
00368       if ( Lambda[0] > Lambda[1] )
00369         {
00370         if ( Lambda[1] > Lambda[2] )  // Most common, guaranteed?
00371           {
00372           ans = Lambda[0] - Lambda[1];
00373           }
00374         else
00375           {
00376           if ( Lambda[0] > Lambda[2] )
00377             {
00378             ans = Lambda[0] - Lambda[2];
00379             }
00380           else
00381             {
00382             ans = Lambda[2] - Lambda[0];
00383             }
00384           }
00385         }
00386       else
00387         {
00388         if ( Lambda[0] > Lambda[2] )
00389           {
00390           ans = Lambda[1] - Lambda[0];
00391           }
00392         else
00393           {
00394           if ( Lambda[1] > Lambda[2] )
00395             {
00396             ans = Lambda[1] - Lambda[2];
00397             }
00398           else
00399             {
00400             ans = Lambda[2] - Lambda[1];
00401             }
00402           }
00403         }
00404       }
00405     else if ( numberOfDistinctRoots == 2 )
00406       {
00407       if ( Lambda[0] > Lambda[1] )
00408         {
00409         ans = Lambda[0] - Lambda[1];
00410         }
00411       else
00412         {
00413         ans = Lambda[1] - Lambda[0];
00414         }
00415       }
00416     else if ( numberOfDistinctRoots == 1 )
00417       {
00418       ans = 0.0;
00419       }
00420     else
00421       {
00422       itkExceptionMacro(<< "Undefined condition. Cubic root solver returned "
00423                         << numberOfDistinctRoots << " distinct roots.");
00424       }
00425 
00426     return ans;
00427   }
00428 
00429   // Function is defined here because the templating confuses gcc 2.96 when
00430   // defined
00431   // in .hxx file.  jc 1/29/03
00432   TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType & it) const
00433   {
00434     unsigned int i, j;
00435 
00436     vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
00437     vnl_vector_fixed< TRealType, VectorDimension >
00438     d_phi_du[itk::GetImageDimension < TInputImage > ::ImageDimension];
00439 
00440     // Calculate the directional derivatives for each vector component using
00441     // central differences.
00442     for ( i = 0; i < ImageDimension; i++ )
00443       {
00444       for ( j = 0; j < VectorDimension; j++ )
00445         {
00446         d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00447                          * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
00448         }
00449       }
00450 
00451     // Calculate the symmetric metric tensor g
00452     for ( i = 0; i < ImageDimension; i++ )
00453       {
00454       for ( j = i; j < ImageDimension; j++ )
00455         {
00456         g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00457         }
00458       }
00459 
00460     // Find the eigenvalues of g
00461     vnl_symmetric_eigensystem< TRealType > E(g);
00462 
00463     // Return the difference in length between the first two principle axes.
00464     // Note that other edge strength metrics may be appropriate here instead..
00465     return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
00466   }
00467 
00469   WeightsType m_DerivativeWeights;
00470 
00474   WeightsType m_ComponentWeights;
00475   WeightsType m_SqrtComponentWeights;
00476 private:
00477   bool m_UseImageSpacing;
00478   bool m_UsePrincipleComponents;
00479 
00480   ThreadIdType m_RequestedNumberOfThreads;
00481 
00482   typename ImageBaseType::ConstPointer m_RealValuedInputImage;
00483 
00484   VectorGradientMagnitudeImageFilter(const Self &); //purposely not implemented
00485   void operator=(const Self &);                     //purposely not implemented
00486 
00487   RadiusType m_NeighborhoodRadius;
00488 };
00489 } // end namespace itk
00490 
00491 #ifndef ITK_MANUAL_INSTANTIATION
00492 #include "itkVectorGradientMagnitudeImageFilter.hxx"
00493 #endif
00494 
00495 #endif
00496