ITK  4.9.0
Insight Segmentation and Registration Toolkit
itkVectorGradientMagnitudeImageFilter.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
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 "vnl/vnl_math.h"
29 
30 namespace itk
31 {
134 template< typename TInputImage,
135  typename TRealType = float,
136  typename TOutputImage = Image< TRealType,
137  TInputImage::ImageDimension >
138  >
140  public ImageToImageFilter< TInputImage, TOutputImage >
141 {
142 public:
148 
150  itkNewMacro(Self);
151 
154 
157  typedef typename TOutputImage::PixelType OutputPixelType;
158  typedef typename TInputImage::PixelType InputPixelType;
159 
161  typedef TInputImage InputImageType;
162  typedef TOutputImage OutputImageType;
163  typedef typename InputImageType::Pointer InputImagePointer;
164  typedef typename OutputImageType::Pointer OutputImagePointer;
165 
167  itkStaticConstMacro(ImageDimension, unsigned int,
168  TOutputImage::ImageDimension);
169 
171  itkStaticConstMacro(VectorDimension, unsigned int,
172  InputPixelType::Dimension);
173 
175  typedef TRealType RealType;
180 
185 
188 
197  virtual void GenerateInputRequestedRegion() ITK_OVERRIDE;
198 
204  { this->SetUseImageSpacing(true); }
205 
210  { this->SetUseImageSpacing(false); }
211 
214  void SetUseImageSpacing(bool);
215 
216  itkGetConstMacro(UseImageSpacing, bool);
217 
219 
222  itkSetMacro(DerivativeWeights, WeightsType);
223  itkGetConstReferenceMacro(DerivativeWeights, WeightsType);
225 
228  itkSetMacro(ComponentWeights, WeightsType);
229  itkGetConstReferenceMacro(ComponentWeights, WeightsType);
231 
232 #ifndef ITK_LEGACY_REMOVE
233 
234  virtual const RadiusType GetNeighborhoodRadius() const
235  {
236  RadiusType r1;
237  r1.Fill(1);
238  return r1;
239  }
241 
243  virtual void SetNeighborhoodRadius(const RadiusType) {}
244 #endif
245 
251  itkSetMacro(UsePrincipleComponents, bool);
252  itkGetConstMacro(UsePrincipleComponents, bool);
254  {
255  this->SetUsePrincipleComponents(true);
256  }
258 
260  {
261  this->SetUsePrincipleComponents(false);
262  }
263 
266  static int CubicSolver(double *, double *);
267 
268 #ifdef ITK_USE_CONCEPT_CHECKING
269  // Begin concept checking
270  itkConceptMacro( InputHasNumericTraitsCheck,
272  itkConceptMacro( RealTypeHasNumericTraitsCheck,
274  // End concept checking
275 #endif
276 
277 protected:
280 
284  void BeforeThreadedGenerateData() ITK_OVERRIDE;
285 
297  void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread,
298  ThreadIdType threadId) ITK_OVERRIDE;
299 
300  void PrintSelf(std::ostream & os, Indent indent) const ITK_OVERRIDE;
301 
303 
305  itkGetConstObjectMacro(RealValuedInputImage, ImageBaseType);
306 
308  {
309  unsigned i, j;
310  TRealType dx, sum, accum;
311 
313  for ( i = 0; i < ImageDimension; ++i )
314  {
316  for ( j = 0; j < VectorDimension; ++j )
317  {
319  * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
320  sum += dx * dx;
321  }
322  accum += sum;
323  }
324  return std::sqrt(accum);
325  }
326 
328  {
329  // WARNING: ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
330  unsigned int i, j;
331  double Lambda[3];
332  double CharEqn[3];
333  double ans;
334 
335  vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
336  vnl_vector_fixed< TRealType, VectorDimension >
337  d_phi_du[TInputImage::ImageDimension];
338 
339  // Calculate the directional derivatives for each vector component using
340  // central differences.
341  for ( i = 0; i < ImageDimension; i++ )
342  {
343  for ( j = 0; j < VectorDimension; j++ )
344  {
345  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
346  * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
347  }
348  }
349 
350  // Calculate the symmetric metric tensor g
351  for ( i = 0; i < ImageDimension; i++ )
352  {
353  for ( j = i; j < ImageDimension; j++ )
354  {
355  g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
356  }
357  }
358 
359  // Find the coefficients of the characteristic equation det(g - lambda I)=0
360  // CharEqn[3] = 1.0;
361 
362  CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
363 
364  CharEqn[1] = ( g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2] )
365  - ( g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1] );
366 
367  CharEqn[0] = g[0][0] * ( g[1][2] * g[2][1] - g[1][1] * g[2][2] )
368  + g[1][0] * ( g[2][2] * g[0][1] - g[0][2] * g[2][1] )
369  + g[2][0] * ( g[1][1] * g[0][2] - g[0][1] * g[1][2] );
370 
371  // Find the eigenvalues of g
372  int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
373 
374  // Define gradient magnitude as the difference between two largest
375  // eigenvalues. Other definitions may be appropriate here as well.
376  if ( numberOfDistinctRoots == 3 ) // By far the most common case
377  {
378  if ( Lambda[0] > Lambda[1] )
379  {
380  if ( Lambda[1] > Lambda[2] ) // Most common, guaranteed?
381  {
382  ans = Lambda[0] - Lambda[1];
383  }
384  else
385  {
386  if ( Lambda[0] > Lambda[2] )
387  {
388  ans = Lambda[0] - Lambda[2];
389  }
390  else
391  {
392  ans = Lambda[2] - Lambda[0];
393  }
394  }
395  }
396  else
397  {
398  if ( Lambda[0] > Lambda[2] )
399  {
400  ans = Lambda[1] - Lambda[0];
401  }
402  else
403  {
404  if ( Lambda[1] > Lambda[2] )
405  {
406  ans = Lambda[1] - Lambda[2];
407  }
408  else
409  {
410  ans = Lambda[2] - Lambda[1];
411  }
412  }
413  }
414  }
415  else if ( numberOfDistinctRoots == 2 )
416  {
417  if ( Lambda[0] > Lambda[1] )
418  {
419  ans = Lambda[0] - Lambda[1];
420  }
421  else
422  {
423  ans = Lambda[1] - Lambda[0];
424  }
425  }
426  else if ( numberOfDistinctRoots == 1 )
427  {
428  ans = 0.0;
429  }
430  else
431  {
432  itkExceptionMacro(<< "Undefined condition. Cubic root solver returned "
433  << numberOfDistinctRoots << " distinct roots.");
434  }
435 
436  return ans;
437  }
438 
439  // Function is defined here because the templating confuses gcc 2.96 when
440  // defined
441  // in .hxx file. jc 1/29/03
443  {
444  unsigned int i, j;
445 
446  vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
447  vnl_vector_fixed< TRealType, VectorDimension >
448  d_phi_du[TInputImage::ImageDimension];
449 
450  // Calculate the directional derivatives for each vector component using
451  // central differences.
452  for ( i = 0; i < ImageDimension; i++ )
453  {
454  for ( j = 0; j < VectorDimension; j++ )
455  {
456  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
457  * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
458  }
459  }
460 
461  // Calculate the symmetric metric tensor g
462  for ( i = 0; i < ImageDimension; i++ )
463  {
464  for ( j = i; j < ImageDimension; j++ )
465  {
466  g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
467  }
468  }
469 
470  // Find the eigenvalues of g
471  vnl_symmetric_eigensystem< TRealType > E(g);
472 
473  // Return the difference in length between the first two principle axes.
474  // Note that other edge strength metrics may be appropriate here instead..
475  return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
476  }
477 
480 
486 
487 private:
490 
492 
493  typename ImageBaseType::ConstPointer m_RealValuedInputImage;
494 
495  VectorGradientMagnitudeImageFilter(const Self &) ITK_DELETE_FUNCTION;
496  void operator=(const Self &) ITK_DELETE_FUNCTION;
497 };
498 } // end namespace itk
499 
500 #ifndef ITK_MANUAL_INSTANTIATION
501 #include "itkVectorGradientMagnitudeImageFilter.hxx"
502 #endif
503 
504 #endif
ConstNeighborhoodIterator< RealVectorImageType > ConstNeighborhoodIteratorType
void PrintSelf(std::ostream &os, Indent indent) const override
FixedArray< TRealType, VectorDimension > WeightsType
virtual void GenerateInputRequestedRegion() override
TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType &it) const
virtual void SetUsePrincipleComponents(bool _arg)
The base class for all process objects (source, filters, mappers) in the Insight data processing pipe...
ImageToImageFilter< TInputImage, TOutputImage > Superclass
Base class for all process objects that output image data.
Const version of NeighborhoodIterator, defining iteration of a local N-dimensional neighborhood of pi...
ConstNeighborhoodIteratorType::RadiusType RadiusType
virtual PixelType GetPrevious(const unsigned axis, NeighborIndexType i) const
Computes a scalar, gradient magnitude image from a multiple channel (pixels are vectors) input...
virtual PixelType GetNext(const unsigned axis, NeighborIndexType i) const
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId) override
typedef(Concept::HasNumericTraits< typename InputPixelType::ValueType >) InputHasNumericTraitsCheck
A templated class holding a n-Dimensional vector.
Definition: itkVector.h:62
static int CubicSolver(double *, double *)
Vector< TRealType, InputPixelType::Dimension > RealVectorType
TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const
unsigned int ThreadIdType
Definition: itkIntTypes.h:159
Base class for filters that take an image as input and produce an image as output.
Control indentation during Print() invocation.
Definition: itkIndent.h:49
Superclass::OutputImageRegionType OutputImageRegionType
#define itkConceptMacro(name, concept)
Image< RealVectorType, TInputImage::ImageDimension > RealVectorImageType
Templated n-dimensional image class.
Definition: itkImage.h:75
TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType &it) const