ITK: itkVectorGradientMagnitudeImageFilter.h Source 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