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"
106 * The filter
requires an image with at least two dimensions and a vector
107 * length of at least 2. The theory supports extension to scalar
images, but
108 * the implementation of the
itk vector classes
do not
110 * The
template parameter TRealType must be floating
point (
float or
double) or
111 * a user-defined
"real" numerical type with arithmetic operations defined
112 * sufficient to compute derivatives.
115 * This filter will automatically multithread
if run with
117 * mode. Unfortunately the ND eigen solver used is not thread safe (a special
118 * 3D solver is), so it cannot multithread
for data other than 3D in
123 * [1] G. Sapiro and D. Ringach,
"Anisotropic Diffusion of Multivalued Images
127 * \ingroup GradientFilters
133 * \ingroup ITKImageGradient
135template <
typename TInputImage,
136 typename TRealType = float,
168 static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
169 static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;
185 using typename Superclass::OutputImageRegionType;
209#if !defined(ITK_FUTURE_LEGACY_REMOVE)
218 this->SetUseImageSpacing(
true);
228 this->SetUseImageSpacing(
false);
234#if !defined(ITK_LEGACY_REMOVE)
260#if !defined(ITK_FUTURE_LEGACY_REMOVE)
265 this->SetUsePrincipleComponents(
true);
272 this->SetUsePrincipleComponents(
false);
281#ifdef ITK_USE_CONCEPT_CHECKING
326 TRealType dx, sum, accum = TRealType{};
327 for (i = 0; i < ImageDimension; ++i)
330 for (j = 0; j < VectorDimension; ++j)
332 dx = m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.
GetNext(i)[j] - it.
GetPrevious(i)[j]);
337 return std::sqrt(accum);
349 vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
350 vnl_vector_fixed<TRealType, VectorDimension> d_phi_du[TInputImage::ImageDimension];
354 for (i = 0; i < ImageDimension; ++i)
356 for (j = 0; j < VectorDimension; ++j)
359 m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.
GetNext(i)[j] - it.
GetPrevious(i)[j]);
364 for (i = 0; i < ImageDimension; ++i)
366 for (j = i; j < ImageDimension; ++j)
368 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
375 CharEqn[2] = -(g[0][0] + g[1][1] + g[2][2]);
377 CharEqn[1] = (g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2]) -
378 (g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1]);
380 CharEqn[0] = g[0][0] * (g[1][2] * g[2][1] - g[1][1] * g[2][2]) + g[1][0] * (g[2][2] * g[0][1] - g[0][2] * g[2][1]) +
381 g[2][0] * (g[1][1] * g[0][2] - g[0][1] * g[1][2]);
384 int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
388 if (numberOfDistinctRoots == 3)
390 if (Lambda[0] > Lambda[1])
392 if (Lambda[1] > Lambda[2])
394 ans = Lambda[0] - Lambda[1];
398 if (Lambda[0] > Lambda[2])
400 ans = Lambda[0] - Lambda[2];
404 ans = Lambda[2] - Lambda[0];
410 if (Lambda[0] > Lambda[2])
412 ans = Lambda[1] - Lambda[0];
416 if (Lambda[1] > Lambda[2])
418 ans = Lambda[1] - Lambda[2];
422 ans = Lambda[2] - Lambda[1];
427 else if (numberOfDistinctRoots == 2)
429 if (Lambda[0] > Lambda[1])
431 ans = Lambda[0] - Lambda[1];
435 ans = Lambda[1] - Lambda[0];
438 else if (numberOfDistinctRoots == 1)
444 itkExceptionMacro(
"Undefined condition. Cubic root solver returned " << numberOfDistinctRoots
445 <<
" distinct roots.");
459 vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
460 vnl_vector_fixed<TRealType, VectorDimension> d_phi_du[TInputImage::ImageDimension];
464 for (i = 0; i < ImageDimension; ++i)
466 for (j = 0; j < VectorDimension; ++j)
469 m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.
GetNext(i)[j] - it.
GetPrevious(i)[j]);
474 for (i = 0; i < ImageDimension; ++i)
476 for (j = i; j < ImageDimension; ++j)
478 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
483 vnl_symmetric_eigensystem<TRealType> E(g);
487 return (E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2));
501 bool m_UseImageSpacing{
true };
502 bool m_UsePrincipleComponents{};
510#ifndef ITK_MANUAL_INSTANTIATION
511# include "itkVectorGradientMagnitudeImageFilter.hxx"
Const version of NeighborhoodIterator, defining iteration of a local N-dimensional neighborhood of pi...
PixelType GetPrevious(const unsigned int axis, NeighborIndexType i) const
PixelType GetNext(const unsigned int axis, NeighborIndexType i) const
Simulate a standard C array with copy semantics.
typename OutputImageType::RegionType OutputImageRegionType
Base class for filters that take an image as input and produce an image as output.
Templated n-dimensional image class.
Control indentation during Print() invocation.
Light weight base class for most itk classes.
Defines iteration of a local N-dimensional neighborhood of pixels across an itk::Image.
Virtual class that defines a common interface to all neighborhood operator subtypes.
A light-weight container object for storing an N-dimensional neighborhood of values.
Computes a scalar, gradient magnitude image from a multiple channel (pixels are vectors) input.
void SetUseImageSpacing(bool)
itkGetConstMacro(UsePrincipleComponents, bool)
itkGetConstReferenceMacro(DerivativeWeights, DerivativeWeightsType)
TOutputImage OutputImageType
void SetUsePrincipleComponentsOff()
~VectorGradientMagnitudeImageFilter() override=default
typename OutputImageType::Pointer OutputImagePointer
itkBooleanMacro(UsePrincipleComponents)
void SetUseImageSpacingOn()
void PrintSelf(std::ostream &os, Indent indent) const override
itkGetConstObjectMacro(RealValuedInputImage, RealVectorImageType)
ITK_DISALLOW_COPY_AND_MOVE(VectorGradientMagnitudeImageFilter)
TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType &it) const
itkConceptMacro(InputHasNumericTraitsCheck,(Concept::HasNumericTraits< typename InputPixelType::ValueType >))
VectorGradientMagnitudeImageFilter()
itkGetConstReferenceMacro(ComponentWeights, ComponentWeightsType)
TInputImage InputImageType
void SetUseImageSpacingOff()
itkConceptMacro(RealTypeHasNumericTraitsCheck,(Concept::HasNumericTraits< RealType >))
static int CubicSolver(const double *, double *)
typename TInputImage::PixelType InputPixelType
void GenerateInputRequestedRegion() override
itkGetConstMacro(UseImageSpacing, bool)
typename TOutputImage::PixelType OutputPixelType
typename ConstNeighborhoodIteratorType::RadiusType RadiusType
itkSetMacro(DerivativeWeights, DerivativeWeightsType)
typename InputImageType::Pointer InputImagePointer
TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
itkSetMacro(UsePrincipleComponents, bool)
itkOverrideGetNameOfClassMacro(VectorGradientMagnitudeImageFilter)
void BeforeThreadedGenerateData() override
void SetUsePrincipleComponentsOn()
itkConceptMacro(SameDimensionCheck,(Concept::SameDimension< InputImageDimension, ImageDimension >))
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) override
TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
itkSetMacro(ComponentWeights, ComponentWeightsType)
typename InputImageType::Superclass ImageBaseType
itkBooleanMacro(UseImageSpacing)
A templated class holding a n-Dimensional vector.
BinaryGeneratorImageFilter< TInputImage1, TInputImage2, TOutputImage > Superclass
SmartPointer< Self > Pointer
constexpr unsigned int Dimension
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
*par Constraints *The filter image with at least two dimensions and a vector *length of at least The theory supports extension to scalar but *the implementation of the itk vector classes do not **The template parameter TRealType must be floating so it cannot multithread for data other than in Anisotropic Diffusion of Multivalued Images *with Application to Color IEEE Transactions on Image No pp
*par Constraints *The filter image with at least two dimensions and a vector *length of at least The theory supports extension to scalar but *the implementation of the itk vector classes do not **The template parameter TRealType must be floating so it cannot multithread for data other than in Anisotropic Diffusion of Multivalued Images *with Application to Color IEEE Transactions on Image * Vol
*par Constraints *The filter image with at least two dimensions and a vector *length of at least The theory supports extension to scalar but *the implementation of the itk vector classes do not **The template parameter TRealType must be floating point(float or double) or *a user-defined "real" numerical type with arithmetic operations defined *sufficient to compute derivatives. **\par Performance *This filter will automatically multithread if run with *SetUsePrincipleComponents
*par Constraints *The filter image with at least two dimensions and a vector *length of at least The theory supports extension to scalar but *the implementation of the itk vector classes do not **The template parameter TRealType must be floating so it cannot multithread for data other than in * UsePrincipleComponents
unsigned int ThreadIdType
*par Constraints *The filter image with at least two dimensions and a vector *length of at least The theory supports extension to scalar but *the implementation of the itk vector classes do not **The template parameter TRealType must be floating so it cannot multithread for data other than in Anisotropic Diffusion of Multivalued Images *with Application to Color IEEE Transactions on Image Processing
*par Constraints *The filter image with at least two dimensions and a vector *length of at least The theory supports extension to scalar images