ITK  6.0.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 * https://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
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
30namespace itk
31{
105 * \par Constraints
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
109 *
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.
113 *
114 * \par Performance
115 * This filter will automatically multithread if run with
116 * SetUsePrincipleComponents=Off or on 3D data in UsePrincipleComponents=On
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
120 *
121 * \par References
122 *
123 * [1] G. Sapiro and D. Ringach, "Anisotropic Diffusion of Multivalued Images
124 * with Application to Color Filtering," IEEE Transactions on Image Processing,
125 * Vol 5, No. 11 pp. 1582-1586, 1996
126
127 * \ingroup GradientFilters
128 *
129 * \sa Image
130 * \sa Neighborhood
133 * \ingroup ITKImageGradient
134 */
135template <typename TInputImage,
136 typename TRealType = float,
137 typename TOutputImage = Image<TRealType, TInputImage::ImageDimension>>
138class ITK_TEMPLATE_EXPORT VectorGradientMagnitudeImageFilter : public ImageToImageFilter<TInputImage, TOutputImage>
139{
140public:
142
148
151
154
157 using OutputPixelType = typename TOutputImage::PixelType;
158 using InputPixelType = typename TInputImage::PixelType;
159
161 using InputImageType = TInputImage;
162 using OutputImageType = TOutputImage;
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 typename Superclass::OutputImageRegionType;
186
195 void
197
203 void
205 itkGetConstMacro(UseImageSpacing, bool);
206 itkBooleanMacro(UseImageSpacing);
209#if !defined(ITK_FUTURE_LEGACY_REMOVE)
215 void
217 {
218 this->SetUseImageSpacing(true);
219 }
220
225 void
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
260#if !defined(ITK_FUTURE_LEGACY_REMOVE)
262 void
264 {
265 this->SetUsePrincipleComponents(true);
266 }
267
269 void
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
285 itkConceptMacro(RealTypeHasNumericTraitsCheck, (Concept::HasNumericTraits<RealType>));
286 // End concept checking
287#endif
288
289protected:
292
296 void
298
310 void
311 DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override;
312
313
314 void
315 PrintSelf(std::ostream & os, Indent indent) const override;
316
318
321
322 TRealType
324 {
325 unsigned int i, j;
326 TRealType dx, sum, accum;
327
328 accum = TRealType{};
329 for (i = 0; i < ImageDimension; ++i)
330 {
331 sum = TRealType{};
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);
502private:
503 bool m_UseImageSpacing{ true };
504 bool m_UsePrincipleComponents{};
505
506 ThreadIdType m_RequestedNumberOfWorkUnits{};
507
508 typename RealVectorImageType::ConstPointer m_RealValuedInputImage{};
509};
510} // end namespace itk
511
512#ifndef ITK_MANUAL_INSTANTIATION
513# include "itkVectorGradientMagnitudeImageFilter.hxx"
514#endif
515
516#endif
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.
Definition: itkFixedArray.h:54
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.
Definition: itkImage.h:89
Control indentation during Print() invocation.
Definition: itkIndent.h:50
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.
itkGetConstMacro(UsePrincipleComponents, bool)
itkGetConstReferenceMacro(DerivativeWeights, DerivativeWeightsType)
~VectorGradientMagnitudeImageFilter() override=default
itkBooleanMacro(UsePrincipleComponents)
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 >))
itkGetConstReferenceMacro(ComponentWeights, ComponentWeightsType)
itkConceptMacro(RealTypeHasNumericTraitsCheck,(Concept::HasNumericTraits< RealType >))
static int CubicSolver(const double *, double *)
itkGetConstMacro(UseImageSpacing, bool)
typename ConstNeighborhoodIteratorType::RadiusType RadiusType
itkSetMacro(DerivativeWeights, DerivativeWeightsType)
TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
itkSetMacro(UsePrincipleComponents, bool)
itkOverrideGetNameOfClassMacro(VectorGradientMagnitudeImageFilter)
itkConceptMacro(SameDimensionCheck,(Concept::SameDimension< InputImageDimension, ImageDimension >))
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) override
TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
itkSetMacro(ComponentWeights, ComponentWeightsType)
A templated class holding a n-Dimensional vector.
Definition: itkVector.h:63
BinaryGeneratorImageFilter< TInputImage1, TInputImage2, TOutputImage > Superclass
SmartPointer< Self > Pointer
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
Definition: itkIntTypes.h:102
*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