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,
167 itkStaticConstMacro(ImageDimension,
unsigned int,
168 TOutputImage::ImageDimension);
171 itkStaticConstMacro(VectorDimension,
unsigned int,
172 InputPixelType::Dimension);
195 virtual void GenerateInputRequestedRegion()
202 void SetUseImageSpacingOn()
203 { this->SetUseImageSpacing(
true); }
208 void SetUseImageSpacingOff()
209 { this->SetUseImageSpacing(
false); }
213 void SetUseImageSpacing(
bool);
215 itkGetConstMacro(UseImageSpacing,
bool);
222 itkGetConstReferenceMacro(DerivativeWeights,
WeightsType);
228 itkGetConstReferenceMacro(ComponentWeights,
WeightsType);
232 itkGetConstReferenceMacro(NeighborhoodRadius,
RadiusType);
241 itkSetMacro(UsePrincipleComponents,
bool);
242 itkGetConstMacro(UsePrincipleComponents,
bool);
243 void SetUsePrincipleComponentsOn()
245 this->SetUsePrincipleComponents(
true);
249 void SetUsePrincipleComponentsOff()
251 this->SetUsePrincipleComponents(
false);
256 static int CubicSolver(
double *,
double *);
258 #ifdef ITK_USE_CONCEPT_CHECKING
274 void BeforeThreadedGenerateData();
287 void ThreadedGenerateData(
const OutputImageRegionType & outputRegionForThread,
290 void PrintSelf(std::ostream & os,
Indent indent)
const;
300 TRealType dx, sum, accum;
303 for ( i = 0; i < ImageDimension; ++i )
306 for ( j = 0; j < VectorDimension; ++j )
308 dx = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
314 return vcl_sqrt(accum);
326 vnl_vector_fixed< TRealType, VectorDimension >
331 for ( i = 0; i < ImageDimension; i++ )
333 for ( j = 0; j < VectorDimension; j++ )
335 d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
341 for ( i = 0; i < ImageDimension; i++ )
343 for ( j = i; j < ImageDimension; j++ )
345 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
352 CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
354 CharEqn[1] = ( g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2] )
355 - ( g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1] );
357 CharEqn[0] = g[0][0] * ( g[1][2] * g[2][1] - g[1][1] * g[2][2] )
358 + g[1][0] * ( g[2][2] * g[0][1] - g[0][2] * g[2][1] )
359 + g[2][0] * ( g[1][1] * g[0][2] - g[0][1] * g[1][2] );
362 int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
366 if ( numberOfDistinctRoots == 3 )
368 if ( Lambda[0] > Lambda[1] )
370 if ( Lambda[1] > Lambda[2] )
372 ans = Lambda[0] - Lambda[1];
376 if ( Lambda[0] > Lambda[2] )
378 ans = Lambda[0] - Lambda[2];
382 ans = Lambda[2] - Lambda[0];
388 if ( Lambda[0] > Lambda[2] )
390 ans = Lambda[1] - Lambda[0];
394 if ( Lambda[1] > Lambda[2] )
396 ans = Lambda[1] - Lambda[2];
400 ans = Lambda[2] - Lambda[1];
405 else if ( numberOfDistinctRoots == 2 )
407 if ( Lambda[0] > Lambda[1] )
409 ans = Lambda[0] - Lambda[1];
413 ans = Lambda[1] - Lambda[0];
416 else if ( numberOfDistinctRoots == 1 )
422 itkExceptionMacro(<<
"Undefined condition. Cubic root solver returned "
423 << numberOfDistinctRoots <<
" distinct roots.");
437 vnl_vector_fixed< TRealType, VectorDimension >
442 for ( i = 0; i < ImageDimension; i++ )
444 for ( j = 0; j < VectorDimension; j++ )
446 d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
452 for ( i = 0; i < ImageDimension; i++ )
454 for ( j = i; j < ImageDimension; j++ )
456 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
461 vnl_symmetric_eigensystem< TRealType > E(g);
465 return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
485 void operator=(
const Self &);
491 #ifndef ITK_MANUAL_INSTANTIATION
492 #include "itkVectorGradientMagnitudeImageFilter.hxx"