ITK  5.2.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_ASSIGN(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 
167  static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;
168 
170  static constexpr unsigned int VectorDimension = InputPixelType::Dimension;
171 
173  using RealType = TRealType;
176 
181 
183  using OutputImageRegionType = typename Superclass::OutputImageRegionType;
184 
193  void
194  GenerateInputRequestedRegion() override;
195 
200  void
202  {
203  this->SetUseImageSpacing(true);
204  }
205 
209  void
211  {
212  this->SetUseImageSpacing(false);
213  }
214 
217  void
218  SetUseImageSpacing(bool);
219 
220  itkGetConstMacro(UseImageSpacing, bool);
221 
223 
226  itkSetMacro(DerivativeWeights, WeightsType);
227  itkGetConstReferenceMacro(DerivativeWeights, WeightsType);
229 
232  itkSetMacro(ComponentWeights, WeightsType);
233  itkGetConstReferenceMacro(ComponentWeights, WeightsType);
235 
241  itkSetMacro(UsePrincipleComponents, bool);
242  itkGetConstMacro(UsePrincipleComponents, bool);
243  void
245  {
246  this->SetUsePrincipleComponents(true);
247  }
249 
250  void
252  {
253  this->SetUsePrincipleComponents(false);
254  }
255 
258  static int
259  CubicSolver(double *, double *);
260 
261 #ifdef ITK_USE_CONCEPT_CHECKING
262  // Begin concept checking
264  itkConceptMacro(RealTypeHasNumericTraitsCheck, (Concept::HasNumericTraits<RealType>));
265  // End concept checking
266 #endif
267 
268 protected:
270  ~VectorGradientMagnitudeImageFilter() override = default;
271 
275  void
276  BeforeThreadedGenerateData() override;
277 
289  void
290  DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override;
291 
292 
293  void
294  PrintSelf(std::ostream & os, Indent indent) const override;
295 
296  using ImageBaseType = typename InputImageType::Superclass;
297 
299  itkGetConstObjectMacro(RealValuedInputImage, RealVectorImageType);
300 
301  TRealType
303  {
304  unsigned i, j;
305  TRealType dx, sum, accum;
306 
308  for (i = 0; i < ImageDimension; ++i)
309  {
311  for (j = 0; j < VectorDimension; ++j)
312  {
313  dx = m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
314  sum += dx * dx;
315  }
316  accum += sum;
317  }
318  return std::sqrt(accum);
319  }
320 
321  TRealType
323  {
324  // WARNING: ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
325  unsigned int i, j;
326  double Lambda[3];
327  double CharEqn[3];
328  double ans;
329 
330  vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
331  vnl_vector_fixed<TRealType, VectorDimension> d_phi_du[TInputImage::ImageDimension];
332 
333  // Calculate the directional derivatives for each vector component using
334  // central differences.
335  for (i = 0; i < ImageDimension; i++)
336  {
337  for (j = 0; j < VectorDimension; j++)
338  {
339  d_phi_du[i][j] =
340  m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
341  }
342  }
343 
344  // Calculate the symmetric metric tensor g
345  for (i = 0; i < ImageDimension; i++)
346  {
347  for (j = i; j < ImageDimension; j++)
348  {
349  g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
350  }
351  }
352 
353  // Find the coefficients of the characteristic equation det(g - lambda I)=0
354  // CharEqn[3] = 1.0;
355 
356  CharEqn[2] = -(g[0][0] + g[1][1] + g[2][2]);
357 
358  CharEqn[1] = (g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2]) -
359  (g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1]);
360 
361  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]) +
362  g[2][0] * (g[1][1] * g[0][2] - g[0][1] * g[1][2]);
363 
364  // Find the eigenvalues of g
365  int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
366 
367  // Define gradient magnitude as the difference between two largest
368  // eigenvalues. Other definitions may be appropriate here as well.
369  if (numberOfDistinctRoots == 3) // By far the most common case
370  {
371  if (Lambda[0] > Lambda[1])
372  {
373  if (Lambda[1] > Lambda[2]) // Most common, guaranteed?
374  {
375  ans = Lambda[0] - Lambda[1];
376  }
377  else
378  {
379  if (Lambda[0] > Lambda[2])
380  {
381  ans = Lambda[0] - Lambda[2];
382  }
383  else
384  {
385  ans = Lambda[2] - Lambda[0];
386  }
387  }
388  }
389  else
390  {
391  if (Lambda[0] > Lambda[2])
392  {
393  ans = Lambda[1] - Lambda[0];
394  }
395  else
396  {
397  if (Lambda[1] > Lambda[2])
398  {
399  ans = Lambda[1] - Lambda[2];
400  }
401  else
402  {
403  ans = Lambda[2] - Lambda[1];
404  }
405  }
406  }
407  }
408  else if (numberOfDistinctRoots == 2)
409  {
410  if (Lambda[0] > Lambda[1])
411  {
412  ans = Lambda[0] - Lambda[1];
413  }
414  else
415  {
416  ans = Lambda[1] - Lambda[0];
417  }
418  }
419  else if (numberOfDistinctRoots == 1)
420  {
421  ans = 0.0;
422  }
423  else
424  {
425  itkExceptionMacro(<< "Undefined condition. Cubic root solver returned " << numberOfDistinctRoots
426  << " distinct roots.");
427  }
428 
429  return ans;
430  }
431 
432  // Function is defined here because the templating confuses gcc 2.96 when
433  // defined
434  // in .hxx file. jc 1/29/03
435  TRealType
437  {
438  unsigned int i, j;
439 
440  vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
441  vnl_vector_fixed<TRealType, VectorDimension> d_phi_du[TInputImage::ImageDimension];
442 
443  // Calculate the directional derivatives for each vector component using
444  // central differences.
445  for (i = 0; i < ImageDimension; i++)
446  {
447  for (j = 0; j < VectorDimension; j++)
448  {
449  d_phi_du[i][j] =
450  m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
451  }
452  }
453 
454  // Calculate the symmetric metric tensor g
455  for (i = 0; i < ImageDimension; i++)
456  {
457  for (j = i; j < ImageDimension; j++)
458  {
459  g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
460  }
461  }
462 
463  // Find the eigenvalues of g
464  vnl_symmetric_eigensystem<TRealType> E(g);
465 
466  // Return the difference in length between the first two principle axes.
467  // Note that other edge strength metrics may be appropriate here instead..
468  return (E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2));
469  }
470 
473 
479 
480 private:
483 
485 
487 };
488 } // end namespace itk
489 
490 #ifndef ITK_MANUAL_INSTANTIATION
491 # include "itkVectorGradientMagnitudeImageFilter.hxx"
492 #endif
493 
494 #endif
itk::VectorGradientMagnitudeImageFilter::m_DerivativeWeights
WeightsType m_DerivativeWeights
Definition: itkVectorGradientMagnitudeImageFilter.h:472
itk::VectorGradientMagnitudeImageFilter::RadiusType
typename ConstNeighborhoodIteratorType::RadiusType RadiusType
Definition: itkVectorGradientMagnitudeImageFilter.h:180
itk::VectorGradientMagnitudeImageFilter::m_RequestedNumberOfThreads
ThreadIdType m_RequestedNumberOfThreads
Definition: itkVectorGradientMagnitudeImageFilter.h:484
itk::ConstNeighborhoodIterator::RadiusType
typename Superclass::RadiusType RadiusType
Definition: itkConstNeighborhoodIterator.h:70
itk::Concept::HasNumericTraits
Definition: itkConceptChecking.h:712
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::InputPixelType
typename TInputImage::PixelType InputPixelType
Definition: itkVectorGradientMagnitudeImageFilter.h:158
itkImage.h
itk::VectorGradientMagnitudeImageFilter::SetUsePrincipleComponentsOff
void SetUsePrincipleComponentsOff()
Definition: itkVectorGradientMagnitudeImageFilter.h:251
itk::SmartPointer< Self >
itk::Indent
Control indentation during Print() invocation.
Definition: itkIndent.h:49
itk::ConstNeighborhoodIterator::GetPrevious
ITK_ITERATOR_VIRTUAL PixelType GetPrevious(const unsigned axis, NeighborIndexType i) const ITK_ITERATOR_FINAL
Definition: itkConstNeighborhoodIterator.h:260
itk::VectorGradientMagnitudeImageFilter::SetUseImageSpacingOn
void SetUseImageSpacingOn()
Definition: itkVectorGradientMagnitudeImageFilter.h:201
itk::VectorGradientMagnitudeImageFilter::m_UsePrincipleComponents
bool m_UsePrincipleComponents
Definition: itkVectorGradientMagnitudeImageFilter.h:482
itk::VectorGradientMagnitudeImageFilter::m_SqrtComponentWeights
WeightsType m_SqrtComponentWeights
Definition: itkVectorGradientMagnitudeImageFilter.h:478
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:436
itk::ImageToImageFilter::InputImagePointer
typename InputImageType::Pointer InputImagePointer
Definition: itkImageToImageFilter.h:130
itk::VectorGradientMagnitudeImageFilter::SetUsePrincipleComponentsOn
void SetUsePrincipleComponentsOn()
Definition: itkVectorGradientMagnitudeImageFilter.h:244
itk::ImageToImageFilter::InputImageType
TInputImage InputImageType
Definition: itkImageToImageFilter.h:129
itkImageToImageFilter.h
itk::VectorGradientMagnitudeImageFilter::EvaluateAtNeighborhood3D
TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType &it) const
Definition: itkVectorGradientMagnitudeImageFilter.h:322
itk::VectorGradientMagnitudeImageFilter::RealType
TRealType RealType
Definition: itkVectorGradientMagnitudeImageFilter.h:173
itk::FixedArray< TRealType, VectorDimension >
itk::VectorGradientMagnitudeImageFilter::SetUseImageSpacingOff
void SetUseImageSpacingOff()
Definition: itkVectorGradientMagnitudeImageFilter.h:210
itk::VectorGradientMagnitudeImageFilter::ImageBaseType
typename InputImageType::Superclass ImageBaseType
Definition: itkVectorGradientMagnitudeImageFilter.h:296
itk::VectorGradientMagnitudeImageFilter::m_UseImageSpacing
bool m_UseImageSpacing
Definition: itkVectorGradientMagnitudeImageFilter.h:481
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:64
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:50
itk::VectorGradientMagnitudeImageFilter::m_RealValuedInputImage
RealVectorImageType::ConstPointer m_RealValuedInputImage
Definition: itkVectorGradientMagnitudeImageFilter.h:486
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::VectorGradientMagnitudeImageFilter::m_ComponentWeights
WeightsType m_ComponentWeights
Definition: itkVectorGradientMagnitudeImageFilter.h:477
itk::ConstNeighborhoodIterator::GetNext
ITK_ITERATOR_VIRTUAL PixelType GetNext(const unsigned axis, NeighborIndexType i) const ITK_ITERATOR_FINAL
Definition: itkConstNeighborhoodIterator.h:242
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itk::VectorGradientMagnitudeImageFilter::NonPCEvaluateAtNeighborhood
TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
Definition: itkVectorGradientMagnitudeImageFilter.h:302
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itkMath.h
itk::ImageSource::OutputImageType
TOutputImage OutputImageType
Definition: itkImageSource.h:90