ITK  5.3.0
Insight Toolkit
itkVectorGradientMagnitudeImageFilter.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright NumFOCUS
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef itkVectorGradientMagnitudeImageFilter_h
19 #define itkVectorGradientMagnitudeImageFilter_h
20 
22 #include "itkImageToImageFilter.h"
23 #include "itkImage.h"
24 #include "itkVector.h"
25 #include "vnl/vnl_matrix.h"
26 #include "vnl/vnl_vector_fixed.h"
27 #include "vnl/algo/vnl_symmetric_eigensystem.h"
28 #include "itkMath.h"
29 
30 namespace itk
31 {
104 
135 template <typename TInputImage,
136  typename TRealType = float,
137  typename TOutputImage = Image<TRealType, TInputImage::ImageDimension>>
138 class ITK_TEMPLATE_EXPORT VectorGradientMagnitudeImageFilter : public ImageToImageFilter<TInputImage, TOutputImage>
139 {
140 public:
141  ITK_DISALLOW_COPY_AND_MOVE(VectorGradientMagnitudeImageFilter);
142 
148 
150  itkNewMacro(Self);
151 
154 
157  using OutputPixelType = typename TOutputImage::PixelType;
158  using InputPixelType = typename TInputImage::PixelType;
159 
161  using InputImageType = TInputImage;
162  using OutputImageType = TOutputImage;
163  using InputImagePointer = typename InputImageType::Pointer;
164  using OutputImagePointer = typename OutputImageType::Pointer;
165 
168  static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
169  static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;
170 
172  static constexpr unsigned int VectorDimension = InputPixelType::Dimension;
173 
175  using RealType = TRealType;
178 
183 
185  using OutputImageRegionType = typename Superclass::OutputImageRegionType;
186 
195  void
196  GenerateInputRequestedRegion() override;
197 
203  void
204  SetUseImageSpacing(bool);
205  itkGetConstMacro(UseImageSpacing, bool);
206  itkBooleanMacro(UseImageSpacing);
208 
209 #if !defined(ITK_FUTURE_LEGACY_REMOVE)
210 
215  void
216  SetUseImageSpacingOn()
217  {
218  this->SetUseImageSpacing(true);
219  }
220 
225  void
226  SetUseImageSpacingOff()
227  {
228  this->SetUseImageSpacing(false);
229  }
230 #endif
231 
234 #if !defined(ITK_LEGACY_REMOVE)
235  using WeightsType [[deprecated("Use DerivativeWeightsType or ComponentWeightsType instead.")]] = ComponentWeightsType;
236 #endif
237 
240  itkSetMacro(DerivativeWeights, DerivativeWeightsType);
241  itkGetConstReferenceMacro(DerivativeWeights, DerivativeWeightsType);
243 
246  itkSetMacro(ComponentWeights, ComponentWeightsType);
247  itkGetConstReferenceMacro(ComponentWeights, ComponentWeightsType);
249 
255  itkSetMacro(UsePrincipleComponents, bool);
256  itkGetConstMacro(UsePrincipleComponents, bool);
257  itkBooleanMacro(UsePrincipleComponents);
259 
260 #if !defined(ITK_FUTURE_LEGACY_REMOVE)
261 
262  void
263  SetUsePrincipleComponentsOn()
264  {
265  this->SetUsePrincipleComponents(true);
266  }
267 
269  void
270  SetUsePrincipleComponentsOff()
271  {
272  this->SetUsePrincipleComponents(false);
273  }
274 #endif
275 
278  static int
279  CubicSolver(const double *, double *);
280 
281 #ifdef ITK_USE_CONCEPT_CHECKING
282  // Begin concept checking
283  itkConceptMacro(SameDimensionCheck, (Concept::SameDimension<InputImageDimension, ImageDimension>));
284  itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits<typename InputPixelType::ValueType>));
285  itkConceptMacro(RealTypeHasNumericTraitsCheck, (Concept::HasNumericTraits<RealType>));
286  // End concept checking
287 #endif
288 
289 protected:
290  VectorGradientMagnitudeImageFilter();
291  ~VectorGradientMagnitudeImageFilter() override = default;
292 
296  void
297  BeforeThreadedGenerateData() override;
298 
310  void
311  DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override;
312 
313 
314  void
315  PrintSelf(std::ostream & os, Indent indent) const override;
316 
317  using ImageBaseType = typename InputImageType::Superclass;
318 
320  itkGetConstObjectMacro(RealValuedInputImage, RealVectorImageType);
321 
322  TRealType
324  {
325  unsigned i, j;
326  TRealType dx, sum, accum;
327 
329  for (i = 0; i < ImageDimension; ++i)
330  {
332  for (j = 0; j < VectorDimension; ++j)
333  {
334  dx = m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
335  sum += dx * dx;
336  }
337  accum += sum;
338  }
339  return std::sqrt(accum);
340  }
341 
342  TRealType
344  {
345  // WARNING: ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
346  unsigned int i, j;
347  double Lambda[3];
348  double CharEqn[3];
349  double ans;
350 
351  vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
352  vnl_vector_fixed<TRealType, VectorDimension> d_phi_du[TInputImage::ImageDimension];
353 
354  // Calculate the directional derivatives for each vector component using
355  // central differences.
356  for (i = 0; i < ImageDimension; ++i)
357  {
358  for (j = 0; j < VectorDimension; ++j)
359  {
360  d_phi_du[i][j] =
361  m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
362  }
363  }
364 
365  // Calculate the symmetric metric tensor g
366  for (i = 0; i < ImageDimension; ++i)
367  {
368  for (j = i; j < ImageDimension; ++j)
369  {
370  g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
371  }
372  }
373 
374  // Find the coefficients of the characteristic equation det(g - lambda I)=0
375  // CharEqn[3] = 1.0;
376 
377  CharEqn[2] = -(g[0][0] + g[1][1] + g[2][2]);
378 
379  CharEqn[1] = (g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2]) -
380  (g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1]);
381 
382  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]) +
383  g[2][0] * (g[1][1] * g[0][2] - g[0][1] * g[1][2]);
384 
385  // Find the eigenvalues of g
386  int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
387 
388  // Define gradient magnitude as the difference between two largest
389  // eigenvalues. Other definitions may be appropriate here as well.
390  if (numberOfDistinctRoots == 3) // By far the most common case
391  {
392  if (Lambda[0] > Lambda[1])
393  {
394  if (Lambda[1] > Lambda[2]) // Most common, guaranteed?
395  {
396  ans = Lambda[0] - Lambda[1];
397  }
398  else
399  {
400  if (Lambda[0] > Lambda[2])
401  {
402  ans = Lambda[0] - Lambda[2];
403  }
404  else
405  {
406  ans = Lambda[2] - Lambda[0];
407  }
408  }
409  }
410  else
411  {
412  if (Lambda[0] > Lambda[2])
413  {
414  ans = Lambda[1] - Lambda[0];
415  }
416  else
417  {
418  if (Lambda[1] > Lambda[2])
419  {
420  ans = Lambda[1] - Lambda[2];
421  }
422  else
423  {
424  ans = Lambda[2] - Lambda[1];
425  }
426  }
427  }
428  }
429  else if (numberOfDistinctRoots == 2)
430  {
431  if (Lambda[0] > Lambda[1])
432  {
433  ans = Lambda[0] - Lambda[1];
434  }
435  else
436  {
437  ans = Lambda[1] - Lambda[0];
438  }
439  }
440  else if (numberOfDistinctRoots == 1)
441  {
442  ans = 0.0;
443  }
444  else
445  {
446  itkExceptionMacro(<< "Undefined condition. Cubic root solver returned " << numberOfDistinctRoots
447  << " distinct roots.");
448  }
449 
450  return ans;
451  }
452 
453  // Function is defined here because the templating confuses gcc 2.96 when
454  // defined
455  // in .hxx file. jc 1/29/03
456  TRealType
458  {
459  unsigned int i, j;
460 
461  vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
462  vnl_vector_fixed<TRealType, VectorDimension> d_phi_du[TInputImage::ImageDimension];
463 
464  // Calculate the directional derivatives for each vector component using
465  // central differences.
466  for (i = 0; i < ImageDimension; ++i)
467  {
468  for (j = 0; j < VectorDimension; ++j)
469  {
470  d_phi_du[i][j] =
471  m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
472  }
473  }
474 
475  // Calculate the symmetric metric tensor g
476  for (i = 0; i < ImageDimension; ++i)
477  {
478  for (j = i; j < ImageDimension; ++j)
479  {
480  g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
481  }
482  }
483 
484  // Find the eigenvalues of g
485  vnl_symmetric_eigensystem<TRealType> E(g);
486 
487  // Return the difference in length between the first two principle axes.
488  // Note that other edge strength metrics may be appropriate here instead..
489  return (E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2));
490  }
491 
493  DerivativeWeightsType m_DerivativeWeights = DerivativeWeightsType::Filled(1);
494 
498  ComponentWeightsType m_ComponentWeights = ComponentWeightsType::Filled(1);
499  ComponentWeightsType m_SqrtComponentWeights = ComponentWeightsType::Filled(1);
501 
502 private:
505 
507 
509 };
510 } // end namespace itk
511 
512 #ifndef ITK_MANUAL_INSTANTIATION
513 # include "itkVectorGradientMagnitudeImageFilter.hxx"
514 #endif
515 
516 #endif
itk::VectorGradientMagnitudeImageFilter::RadiusType
typename ConstNeighborhoodIteratorType::RadiusType RadiusType
Definition: itkVectorGradientMagnitudeImageFilter.h:182
itk::ConstNeighborhoodIterator::RadiusType
typename Superclass::RadiusType RadiusType
Definition: itkConstNeighborhoodIterator.h:71
itk::ImageSource::OutputImagePointer
typename OutputImageType::Pointer OutputImagePointer
Definition: itkImageSource.h:91
itk::VectorGradientMagnitudeImageFilter
Computes a scalar, gradient magnitude image from a multiple channel (pixels are vectors) input.
Definition: itkVectorGradientMagnitudeImageFilter.h:138
itkNeighborhoodIterator.h
itk::Vector
A templated class holding a n-Dimensional vector.
Definition: itkVector.h:62
itk::VectorGradientMagnitudeImageFilter::m_RequestedNumberOfWorkUnits
ThreadIdType m_RequestedNumberOfWorkUnits
Definition: itkVectorGradientMagnitudeImageFilter.h:506
itk::VectorGradientMagnitudeImageFilter::InputPixelType
typename TInputImage::PixelType InputPixelType
Definition: itkVectorGradientMagnitudeImageFilter.h:158
itkImage.h
itk::SmartPointer< Self >
itk::ConstNeighborhoodIterator::GetPrevious
ITK_ITERATOR_VIRTUAL PixelType GetPrevious(const unsigned axis, NeighborIndexType i) const ITK_ITERATOR_FINAL
Definition: itkConstNeighborhoodIterator.h:261
itk::VectorGradientMagnitudeImageFilter::m_UsePrincipleComponents
bool m_UsePrincipleComponents
Definition: itkVectorGradientMagnitudeImageFilter.h:504
itk::ThreadIdType
unsigned int ThreadIdType
Definition: itkIntTypes.h:99
itk::ImageToImageFilter
Base class for filters that take an image as input and produce an image as output.
Definition: itkImageToImageFilter.h:108
itk::ImageSource
Base class for all process objects that output image data.
Definition: itkImageSource.h:67
itk::VectorGradientMagnitudeImageFilter::EvaluateAtNeighborhood
TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
Definition: itkVectorGradientMagnitudeImageFilter.h:457
itk::ImageToImageFilter::InputImagePointer
typename InputImageType::Pointer InputImagePointer
Definition: itkImageToImageFilter.h:130
itk::ImageToImageFilter::InputImageType
TInputImage InputImageType
Definition: itkImageToImageFilter.h:129
itkImageToImageFilter.h
itk::VectorGradientMagnitudeImageFilter::EvaluateAtNeighborhood3D
TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType &it) const
Definition: itkVectorGradientMagnitudeImageFilter.h:343
itk::VectorGradientMagnitudeImageFilter::RealType
TRealType RealType
Definition: itkVectorGradientMagnitudeImageFilter.h:175
itk::FixedArray< TRealType, VectorDimension >
itk::VectorGradientMagnitudeImageFilter::ImageBaseType
typename InputImageType::Superclass ImageBaseType
Definition: itkVectorGradientMagnitudeImageFilter.h:317
itk::VectorGradientMagnitudeImageFilter::m_UseImageSpacing
bool m_UseImageSpacing
Definition: itkVectorGradientMagnitudeImageFilter.h:503
itk::ImageSource::OutputImageRegionType
typename OutputImageType::RegionType OutputImageRegionType
Definition: itkImageSource.h:92
itk::NumericTraits::ZeroValue
static T ZeroValue()
Definition: itkNumericTraits.h:148
itkConceptMacro
#define itkConceptMacro(name, concept)
Definition: itkConceptChecking.h:65
itk::VectorGradientMagnitudeImageFilter::OutputPixelType
typename TOutputImage::PixelType OutputPixelType
Definition: itkVectorGradientMagnitudeImageFilter.h:157
itk
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
Definition: itkAnnulusOperator.h:24
itk::ConstNeighborhoodIterator
Const version of NeighborhoodIterator, defining iteration of a local N-dimensional neighborhood of pi...
Definition: itkConstNeighborhoodIterator.h:51
itk::VectorGradientMagnitudeImageFilter::m_RealValuedInputImage
RealVectorImageType::ConstPointer m_RealValuedInputImage
Definition: itkVectorGradientMagnitudeImageFilter.h:508
itk::ProcessObject
The base class for all process objects (source, filters, mappers) in the Insight data processing pipe...
Definition: itkProcessObject.h:138
itkVector.h
itk::ConstNeighborhoodIterator::GetNext
ITK_ITERATOR_VIRTUAL PixelType GetNext(const unsigned axis, NeighborIndexType i) const ITK_ITERATOR_FINAL
Definition: itkConstNeighborhoodIterator.h:243
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:88
itk::VectorGradientMagnitudeImageFilter::NonPCEvaluateAtNeighborhood
TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
Definition: itkVectorGradientMagnitudeImageFilter.h:323
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itkMath.h
itk::ImageSource::OutputImageType
TOutputImage OutputImageType
Definition: itkImageSource.h:90