00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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
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
00326
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
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
00344
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
00356 int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
00357
00358
00359
00360 if (numberOfDistinctRoots == 3)
00361 {
00362 if (Lambda[0] > Lambda[1])
00363 {
00364 if ( Lambda[1] > Lambda[2] )
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
00424
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
00433
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
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
00453 vnl_symmetric_eigensystem<TRealType> E(g);
00454
00455
00456
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&);
00477 void operator=(const Self&);
00478
00479 RadiusType m_NeighborhoodRadius;
00480 };
00481
00482 }
00483
00484 #ifndef ITK_MANUAL_INSTANTIATION
00485 #include "itkVectorGradientMagnitudeImageFilter.txx"
00486 #endif
00487
00488 #endif
00489