ITK  4.3.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 
187  typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
188 
197  virtual void GenerateInputRequestedRegion()
199 
204  void SetUseImageSpacingOn()
205  { this->SetUseImageSpacing(true); }
206 
210  void SetUseImageSpacingOff()
211  { this->SetUseImageSpacing(false); }
212 
215  void SetUseImageSpacing(bool);
216 
217  itkGetConstMacro(UseImageSpacing, bool);
218 
220 
223  itkSetMacro(DerivativeWeights, WeightsType);
224  itkGetConstReferenceMacro(DerivativeWeights, WeightsType);
226 
229  itkSetMacro(ComponentWeights, WeightsType);
230  itkGetConstReferenceMacro(ComponentWeights, WeightsType);
232 
234  itkGetConstReferenceMacro(NeighborhoodRadius, RadiusType);
235  itkSetMacro(NeighborhoodRadius, RadiusType);
237 
243  itkSetMacro(UsePrincipleComponents, bool);
244  itkGetConstMacro(UsePrincipleComponents, bool);
245  void SetUsePrincipleComponentsOn()
246  {
247  this->SetUsePrincipleComponents(true);
248  }
250 
251  void SetUsePrincipleComponentsOff()
252  {
253  this->SetUsePrincipleComponents(false);
254  }
255 
258  static int CubicSolver(double *, double *);
259 
260 #ifdef ITK_USE_CONCEPT_CHECKING
261 
262  itkConceptMacro( InputHasNumericTraitsCheck,
264  itkConceptMacro( RealTypeHasNumericTraitsCheck,
266 
268 #endif
269 
270 protected:
273 
277  void BeforeThreadedGenerateData();
278 
290  void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread,
291  ThreadIdType threadId);
292 
293  void PrintSelf(std::ostream & os, Indent indent) const;
294 
295  typedef typename InputImageType::Superclass ImageBaseType;
296 
298  itkGetConstObjectMacro(RealValuedInputImage, ImageBaseType);
299 
300  TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType & it) const
301  {
302  unsigned i, j;
303  TRealType dx, sum, accum;
304 
306  for ( i = 0; i < ImageDimension; ++i )
307  {
309  for ( j = 0; j < VectorDimension; ++j )
310  {
311  dx = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
312  * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
313  sum += dx * dx;
314  }
315  accum += sum;
316  }
317  return vcl_sqrt(accum);
318  }
319 
320  TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType & it) const
321  {
322  // WARNING: ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
323  unsigned int i, j;
324  double Lambda[3];
325  double CharEqn[3];
326  double ans;
327 
328  vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
329  vnl_vector_fixed< TRealType, VectorDimension >
330  d_phi_du[TInputImage::ImageDimension];
331 
332  // Calculate the directional derivatives for each vector component using
333  // central differences.
334  for ( i = 0; i < ImageDimension; i++ )
335  {
336  for ( j = 0; j < VectorDimension; j++ )
337  {
338  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
339  * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
340  }
341  }
342 
343  // Calculate the symmetric metric tensor g
344  for ( i = 0; i < ImageDimension; i++ )
345  {
346  for ( j = i; j < ImageDimension; j++ )
347  {
348  g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
349  }
350  }
351 
352  // Find the coefficients of the characteristic equation det(g - lambda I)=0
353  // CharEqn[3] = 1.0;
354 
355  CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
356 
357  CharEqn[1] = ( g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2] )
358  - ( g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1] );
359 
360  CharEqn[0] = g[0][0] * ( g[1][2] * g[2][1] - g[1][1] * g[2][2] )
361  + 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 "
426  << numberOfDistinctRoots << " 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 EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType & it) const
436  {
437  unsigned int i, j;
438 
439  vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
440  vnl_vector_fixed< TRealType, VectorDimension >
441  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] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
450  * 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 
486  typename ImageBaseType::ConstPointer m_RealValuedInputImage;
487 
488  VectorGradientMagnitudeImageFilter(const Self &); //purposely not implemented
489  void operator=(const Self &); //purposely not implemented
490 
492 };
493 } // end namespace itk
494 
495 #ifndef ITK_MANUAL_INSTANTIATION
496 #include "itkVectorGradientMagnitudeImageFilter.hxx"
497 #endif
498 
499 #endif
500