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 = TRealType{};
327 for (i = 0; i < ImageDimension; ++i)
328 {
329 sum = TRealType{};
330 for (j = 0; j < VectorDimension; ++j)
331 {
332 dx = m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
333 sum += dx * dx;
334 }
335 accum += sum;
336 }
337 return std::sqrt(accum);
338 }
339
340 TRealType
342 {
343 // WARNING: ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
344 unsigned int i, j;
345 double Lambda[3];
346 double CharEqn[3];
347 double ans;
348
349 vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
350 vnl_vector_fixed<TRealType, VectorDimension> d_phi_du[TInputImage::ImageDimension];
351
352 // Calculate the directional derivatives for each vector component using
353 // central differences.
354 for (i = 0; i < ImageDimension; ++i)
355 {
356 for (j = 0; j < VectorDimension; ++j)
357 {
358 d_phi_du[i][j] =
359 m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
360 }
361 }
362
363 // Calculate the symmetric metric tensor g
364 for (i = 0; i < ImageDimension; ++i)
365 {
366 for (j = i; j < ImageDimension; ++j)
367 {
368 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
369 }
370 }
371
372 // Find the coefficients of the characteristic equation det(g - lambda I)=0
373 // CharEqn[3] = 1.0;
374
375 CharEqn[2] = -(g[0][0] + g[1][1] + g[2][2]);
376
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]);
379
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]);
382
383 // Find the eigenvalues of g
384 int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
385
386 // Define gradient magnitude as the difference between two largest
387 // eigenvalues. Other definitions may be appropriate here as well.
388 if (numberOfDistinctRoots == 3) // By far the most common case
389 {
390 if (Lambda[0] > Lambda[1])
391 {
392 if (Lambda[1] > Lambda[2]) // Most common, guaranteed?
393 {
394 ans = Lambda[0] - Lambda[1];
395 }
396 else
397 {
398 if (Lambda[0] > Lambda[2])
399 {
400 ans = Lambda[0] - Lambda[2];
401 }
402 else
403 {
404 ans = Lambda[2] - Lambda[0];
405 }
406 }
407 }
408 else
409 {
410 if (Lambda[0] > Lambda[2])
411 {
412 ans = Lambda[1] - Lambda[0];
413 }
414 else
415 {
416 if (Lambda[1] > Lambda[2])
417 {
418 ans = Lambda[1] - Lambda[2];
419 }
420 else
421 {
422 ans = Lambda[2] - Lambda[1];
423 }
424 }
425 }
426 }
427 else if (numberOfDistinctRoots == 2)
428 {
429 if (Lambda[0] > Lambda[1])
430 {
431 ans = Lambda[0] - Lambda[1];
432 }
433 else
434 {
435 ans = Lambda[1] - Lambda[0];
436 }
437 }
438 else if (numberOfDistinctRoots == 1)
439 {
440 ans = 0.0;
441 }
442 else
443 {
444 itkExceptionMacro("Undefined condition. Cubic root solver returned " << numberOfDistinctRoots
445 << " distinct roots.");
446 }
447
448 return ans;
449 }
450
451 // Function is defined here because the templating confuses gcc 2.96 when
452 // defined
453 // in .hxx file. jc 1/29/03
454 TRealType
456 {
457 unsigned int i, j;
458
459 vnl_matrix<TRealType> g(ImageDimension, ImageDimension);
460 vnl_vector_fixed<TRealType, VectorDimension> d_phi_du[TInputImage::ImageDimension];
461
462 // Calculate the directional derivatives for each vector component using
463 // central differences.
464 for (i = 0; i < ImageDimension; ++i)
465 {
466 for (j = 0; j < VectorDimension; ++j)
467 {
468 d_phi_du[i][j] =
469 m_DerivativeWeights[i] * m_SqrtComponentWeights[j] * 0.5 * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
470 }
471 }
472
473 // Calculate the symmetric metric tensor g
474 for (i = 0; i < ImageDimension; ++i)
475 {
476 for (j = i; j < ImageDimension; ++j)
477 {
478 g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
479 }
480 }
481
482 // Find the eigenvalues of g
483 vnl_symmetric_eigensystem<TRealType> E(g);
484
485 // Return the difference in length between the first two principle axes.
486 // Note that other edge strength metrics may be appropriate here instead..
487 return (E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2));
488 }
489
491 DerivativeWeightsType m_DerivativeWeights = DerivativeWeightsType::Filled(1);
492
496 ComponentWeightsType m_ComponentWeights = ComponentWeightsType::Filled(1);
497 ComponentWeightsType m_SqrtComponentWeights = ComponentWeightsType::Filled(1);
500private:
501 bool m_UseImageSpacing{ true };
502 bool m_UsePrincipleComponents{};
503
504 ThreadIdType m_RequestedNumberOfWorkUnits{};
505
506 typename RealVectorImageType::ConstPointer m_RealValuedInputImage{};
507};
508} // end namespace itk
509
510#ifndef ITK_MANUAL_INSTANTIATION
511# include "itkVectorGradientMagnitudeImageFilter.hxx"
512#endif
513
514#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