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/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 typedef typename ConstNeighborhoodIteratorType::RadiusType RadiusType;
00182
00184 typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
00185
00194 virtual void GenerateInputRequestedRegion() throw(InvalidRequestedRegionError);
00195
00199 void SetUseImageSpacingOn()
00200 { this->SetUseImageSpacing(true); }
00201
00205 void SetUseImageSpacingOff()
00206 { this->SetUseImageSpacing(false); }
00207
00210 void SetUseImageSpacing(bool);
00211 itkGetMacro(UseImageSpacing, bool);
00212
00215 void SetDerivativeWeights(TRealType data[]);
00216 itkGetVectorMacro(DerivativeWeights, const TRealType, itk::GetImageDimension<TInputImage>::ImageDimension);
00217
00220 itkSetVectorMacro(ComponentWeights, TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
00221 itkGetVectorMacro(ComponentWeights, const TRealType, itk::GetVectorDimension<InputPixelType>::VectorDimension);
00222
00228 itkSetMacro(UsePrincipleComponents, bool);
00229 itkGetMacro(UsePrincipleComponents, bool);
00230 void SetUsePrincipleComponentsOn()
00231 {
00232 this->SetUsePrincipleComponents(true);
00233 }
00234 void SetUsePrincipleComponentsOff()
00235 {
00236 this->SetUsePrincipleComponents(false);
00237 }
00238
00241 static int CubicSolver(double *, double *);
00242
00243 protected:
00244 VectorGradientMagnitudeImageFilter();
00245 virtual ~VectorGradientMagnitudeImageFilter() {}
00246
00250 void BeforeThreadedGenerateData ();
00251
00263 void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
00264 int threadId );
00265
00266 void PrintSelf(std::ostream& os, Indent indent) const;
00267
00268 typedef typename InputImageType::Superclass ImageBaseType;
00269
00271 itkGetConstObjectMacro( RealValuedInputImage, ImageBaseType );
00272
00274 itkGetConstReferenceMacro( NeighborhoodRadius, RadiusType );
00275 itkSetMacro( NeighborhoodRadius, RadiusType );
00276
00277
00278 TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
00279 {
00280 unsigned i, j;
00281 TRealType dx, sum, accum;
00282 accum = NumericTraits<TRealType>::Zero;
00283 for (i = 0; i < ImageDimension; ++i)
00284 {
00285 sum = NumericTraits<TRealType>::Zero;
00286 for (j = 0; j < VectorDimension; ++j)
00287 {
00288 dx = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00289 * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
00290 sum += dx * dx;
00291 }
00292 accum += sum;
00293 }
00294 return vcl_sqrt(accum);
00295 }
00296
00297 TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType &it) const
00298 {
00299
00300 unsigned int i, j;
00301 double Lambda[3];
00302 double CharEqn[3];
00303 double ans;
00304
00305 vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
00306 vnl_vector_fixed<TRealType, VectorDimension>
00307 d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
00308
00309
00310
00311 for (i = 0; i < ImageDimension; i++)
00312 {
00313 for (j = 0; j < VectorDimension; j++)
00314 { d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00315 * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]); }
00316 }
00317
00318
00319 for (i = 0; i < ImageDimension; i++)
00320 {
00321 for (j = i; j < ImageDimension; j++)
00322 {
00323 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00324 }
00325 }
00326
00327
00328
00329
00330 CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
00331
00332 CharEqn[1] =(g[0][0]*g[1][1] + g[0][0]*g[2][2] + g[1][1]*g[2][2])
00333 - (g[0][1]*g[1][0] + g[0][2]*g[2][0] + g[1][2]*g[2][1]);
00334
00335 CharEqn[0] = g[0][0] * ( g[1][2]*g[2][1] - g[1][1]*g[2][2] ) +
00336 g[1][0] * ( g[2][2]*g[0][1] - g[0][2]*g[2][1] ) +
00337 g[2][0] * ( g[1][1]*g[0][2] - g[0][1]*g[1][2] );
00338
00339
00340
00341
00342
00343 int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
00344
00345
00346
00347 if (numberOfDistinctRoots == 3)
00348 {
00349 if (Lambda[0] > Lambda[1])
00350 {
00351 if ( Lambda[1] > Lambda[2] )
00352 { ans = Lambda[0] - Lambda[1]; }
00353 else
00354 {
00355 if ( Lambda[0] > Lambda[2] )
00356 { ans = Lambda[0] - Lambda[2]; }
00357 else
00358 { ans = Lambda[2] - Lambda[0]; }
00359 }
00360 }
00361 else
00362 {
00363 if ( Lambda[0] > Lambda[2] )
00364 { ans = Lambda[1] - Lambda[0]; }
00365 else
00366 {
00367 if ( Lambda[1] > Lambda[2] )
00368 { ans = Lambda[1] - Lambda[2]; }
00369 else
00370 { ans = Lambda[2] - Lambda[1]; }
00371 }
00372 }
00373 }
00374 else if (numberOfDistinctRoots == 2)
00375 {
00376 if ( Lambda[0] > Lambda[1] )
00377 { ans = Lambda[0] - Lambda[1]; }
00378 else
00379 { ans = Lambda[1] - Lambda[0]; }
00380 }
00381 else if (numberOfDistinctRoots == 1)
00382 {
00383 ans = 0.0;
00384 }
00385 else
00386 {
00387 itkExceptionMacro( << "Undefined condition. Cubic root solver returned "
00388 << numberOfDistinctRoots << " distinct roots." );
00389 }
00390
00391 return ans;
00392 }
00393
00394
00395
00396 TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
00397 {
00398 unsigned int i, j;
00399 vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
00400 vnl_vector_fixed<TRealType, VectorDimension>
00401 d_phi_du[itk::GetImageDimension<TInputImage>::ImageDimension];
00402
00403
00404
00405 for (i = 0; i < ImageDimension; i++)
00406 {
00407 for (j = 0; j < VectorDimension; j++)
00408 { d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
00409 * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j] ); }
00410 }
00411
00412
00413 for (i = 0; i < ImageDimension; i++)
00414 {
00415 for (j = i; j < ImageDimension; j++)
00416 {
00417 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
00418 }
00419 }
00420
00421
00422 vnl_symmetric_eigensystem<TRealType> E(g);
00423
00424
00425
00426 return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
00427 }
00428
00430 TRealType m_DerivativeWeights[itk::GetImageDimension<TInputImage>::ImageDimension];
00431
00435 TRealType m_ComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00436 TRealType m_SqrtComponentWeights[itk::GetVectorDimension<InputPixelType>::VectorDimension];
00437
00438 private:
00439 bool m_UseImageSpacing;
00440 bool m_UsePrincipleComponents;
00441 int m_RequestedNumberOfThreads;
00442
00443
00444 typename ImageBaseType::ConstPointer m_RealValuedInputImage;
00445
00446 VectorGradientMagnitudeImageFilter(const Self&);
00447 void operator=(const Self&);
00448
00449 RadiusType m_NeighborhoodRadius;
00450 };
00451
00452 }
00453
00454 #ifndef ITK_MANUAL_INSTANTIATION
00455 #include "itkVectorGradientMagnitudeImageFilter.txx"
00456 #endif
00457
00458 #endif