18 #ifndef __itkVectorGradientMagnitudeImageFilter_h
19 #define __itkVectorGradientMagnitudeImageFilter_h
25 #include "vnl/vnl_matrix.h"
26 #include "vnl/vnl_vector_fixed.h"
27 #include "vnl/algo/vnl_symmetric_eigensystem.h"
28 #include "vnl/vnl_math.h"
134 template<
typename TInputImage,
135 typename TRealType = float,
136 typename TOutputImage = Image< TRealType,
137 TInputImage::ImageDimension >
167 itkStaticConstMacro(ImageDimension,
unsigned int,
168 TOutputImage::ImageDimension);
171 itkStaticConstMacro(VectorDimension,
unsigned int,
172 InputPixelType::Dimension);
197 virtual void GenerateInputRequestedRegion()
204 void SetUseImageSpacingOn()
205 { this->SetUseImageSpacing(
true); }
210 void SetUseImageSpacingOff()
211 { this->SetUseImageSpacing(
false); }
215 void SetUseImageSpacing(
bool);
217 itkGetConstMacro(UseImageSpacing,
bool);
224 itkGetConstReferenceMacro(DerivativeWeights,
WeightsType);
230 itkGetConstReferenceMacro(ComponentWeights,
WeightsType);
234 itkGetConstReferenceMacro(NeighborhoodRadius,
RadiusType);
243 itkSetMacro(UsePrincipleComponents,
bool);
244 itkGetConstMacro(UsePrincipleComponents,
bool);
245 void SetUsePrincipleComponentsOn()
247 this->SetUsePrincipleComponents(
true);
251 void SetUsePrincipleComponentsOff()
253 this->SetUsePrincipleComponents(
false);
258 static int CubicSolver(
double *,
double *);
260 #ifdef ITK_USE_CONCEPT_CHECKING
277 void BeforeThreadedGenerateData();
290 void ThreadedGenerateData(
const OutputImageRegionType & outputRegionForThread,
293 void PrintSelf(std::ostream & os,
Indent indent)
const;
303 TRealType dx, sum, accum;
306 for ( i = 0; i < ImageDimension; ++i )
309 for ( j = 0; j < VectorDimension; ++j )
311 dx = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
317 return vcl_sqrt(accum);
329 vnl_vector_fixed< TRealType, VectorDimension >
330 d_phi_du[TInputImage::ImageDimension];
334 for ( i = 0; i < ImageDimension; i++ )
336 for ( j = 0; j < VectorDimension; j++ )
338 d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
344 for ( i = 0; i < ImageDimension; i++ )
346 for ( j = i; j < ImageDimension; j++ )
348 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
355 CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
357 CharEqn[1] = ( g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2] )
358 - ( g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1] );
360 CharEqn[0] = g[0][0] * ( g[1][2] * g[2][1] - g[1][1] * g[2][2] )
361 + g[1][0] * ( g[2][2] * g[0][1] - g[0][2] * g[2][1] )
362 + g[2][0] * ( g[1][1] * g[0][2] - g[0][1] * g[1][2] );
365 int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
369 if ( numberOfDistinctRoots == 3 )
371 if ( Lambda[0] > Lambda[1] )
373 if ( Lambda[1] > Lambda[2] )
375 ans = Lambda[0] - Lambda[1];
379 if ( Lambda[0] > Lambda[2] )
381 ans = Lambda[0] - Lambda[2];
385 ans = Lambda[2] - Lambda[0];
391 if ( Lambda[0] > Lambda[2] )
393 ans = Lambda[1] - Lambda[0];
397 if ( Lambda[1] > Lambda[2] )
399 ans = Lambda[1] - Lambda[2];
403 ans = Lambda[2] - Lambda[1];
408 else if ( numberOfDistinctRoots == 2 )
410 if ( Lambda[0] > Lambda[1] )
412 ans = Lambda[0] - Lambda[1];
416 ans = Lambda[1] - Lambda[0];
419 else if ( numberOfDistinctRoots == 1 )
425 itkExceptionMacro(<<
"Undefined condition. Cubic root solver returned "
426 << numberOfDistinctRoots <<
" distinct roots.");
440 vnl_vector_fixed< TRealType, VectorDimension >
441 d_phi_du[TInputImage::ImageDimension];
445 for ( i = 0; i < ImageDimension; i++ )
447 for ( j = 0; j < VectorDimension; j++ )
449 d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
455 for ( i = 0; i < ImageDimension; i++ )
457 for ( j = i; j < ImageDimension; j++ )
459 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
464 vnl_symmetric_eigensystem< TRealType > E(g);
468 return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
489 void operator=(
const Self &);
495 #ifndef ITK_MANUAL_INSTANTIATION
496 #include "itkVectorGradientMagnitudeImageFilter.hxx"