ITK  4.2.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,
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;
178 
183 
185  typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
186 
195  virtual void GenerateInputRequestedRegion()
197 
202  void SetUseImageSpacingOn()
203  { this->SetUseImageSpacing(true); }
204 
208  void SetUseImageSpacingOff()
209  { this->SetUseImageSpacing(false); }
210 
213  void SetUseImageSpacing(bool);
214 
215  itkGetConstMacro(UseImageSpacing, bool);
216 
218 
221  itkSetMacro(DerivativeWeights, WeightsType);
222  itkGetConstReferenceMacro(DerivativeWeights, WeightsType);
224 
227  itkSetMacro(ComponentWeights, WeightsType);
228  itkGetConstReferenceMacro(ComponentWeights, WeightsType);
230 
232  itkGetConstReferenceMacro(NeighborhoodRadius, RadiusType);
233  itkSetMacro(NeighborhoodRadius, RadiusType);
235 
241  itkSetMacro(UsePrincipleComponents, bool);
242  itkGetConstMacro(UsePrincipleComponents, bool);
243  void SetUsePrincipleComponentsOn()
244  {
245  this->SetUsePrincipleComponents(true);
246  }
248 
249  void SetUsePrincipleComponentsOff()
250  {
251  this->SetUsePrincipleComponents(false);
252  }
253 
256  static int CubicSolver(double *, double *);
257 
258 #ifdef ITK_USE_CONCEPT_CHECKING
259 
260  itkConceptMacro( InputHasNumericTraitsCheck,
262  itkConceptMacro( RealTypeHasNumericTraitsCheck,
264 
266 #endif
267 protected:
270 
274  void BeforeThreadedGenerateData();
275 
287  void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread,
288  ThreadIdType threadId);
289 
290  void PrintSelf(std::ostream & os, Indent indent) const;
291 
292  typedef typename InputImageType::Superclass ImageBaseType;
293 
295  itkGetConstObjectMacro(RealValuedInputImage, ImageBaseType);
296 
297  TRealType NonPCEvaluateAtNeighborhood(const ConstNeighborhoodIteratorType & it) const
298  {
299  unsigned i, j;
300  TRealType dx, sum, accum;
301 
303  for ( i = 0; i < ImageDimension; ++i )
304  {
306  for ( j = 0; j < VectorDimension; ++j )
307  {
308  dx = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
309  * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
310  sum += dx * dx;
311  }
312  accum += sum;
313  }
314  return vcl_sqrt(accum);
315  }
316 
317  TRealType EvaluateAtNeighborhood3D(const ConstNeighborhoodIteratorType & it) const
318  {
319  // WARNING: ONLY CALL THIS METHOD WHEN PROCESSING A 3D IMAGE
320  unsigned int i, j;
321  double Lambda[3];
322  double CharEqn[3];
323  double ans;
324 
325  vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
326  vnl_vector_fixed< TRealType, VectorDimension >
327  d_phi_du[itk::GetImageDimension < TInputImage > ::ImageDimension];
328 
329  // Calculate the directional derivatives for each vector component using
330  // central differences.
331  for ( i = 0; i < ImageDimension; i++ )
332  {
333  for ( j = 0; j < VectorDimension; j++ )
334  {
335  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
336  * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
337  }
338  }
339 
340  // Calculate the symmetric metric tensor g
341  for ( i = 0; i < ImageDimension; i++ )
342  {
343  for ( j = i; j < ImageDimension; j++ )
344  {
345  g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
346  }
347  }
348 
349  // Find the coefficients of the characteristic equation det(g - lambda I)=0
350  // CharEqn[3] = 1.0;
351 
352  CharEqn[2] = -( g[0][0] + g[1][1] + g[2][2] );
353 
354  CharEqn[1] = ( g[0][0] * g[1][1] + g[0][0] * g[2][2] + g[1][1] * g[2][2] )
355  - ( g[0][1] * g[1][0] + g[0][2] * g[2][0] + g[1][2] * g[2][1] );
356 
357  CharEqn[0] = g[0][0] * ( g[1][2] * g[2][1] - g[1][1] * g[2][2] )
358  + g[1][0] * ( g[2][2] * g[0][1] - g[0][2] * g[2][1] )
359  + g[2][0] * ( g[1][1] * g[0][2] - g[0][1] * g[1][2] );
360 
361  // Find the eigenvalues of g
362  int numberOfDistinctRoots = this->CubicSolver(CharEqn, Lambda);
363 
364  // Define gradient magnitude as the difference between two largest
365  // eigenvalues. Other definitions may be appropriate here as well.
366  if ( numberOfDistinctRoots == 3 ) // By far the most common case
367  {
368  if ( Lambda[0] > Lambda[1] )
369  {
370  if ( Lambda[1] > Lambda[2] ) // Most common, guaranteed?
371  {
372  ans = Lambda[0] - Lambda[1];
373  }
374  else
375  {
376  if ( Lambda[0] > Lambda[2] )
377  {
378  ans = Lambda[0] - Lambda[2];
379  }
380  else
381  {
382  ans = Lambda[2] - Lambda[0];
383  }
384  }
385  }
386  else
387  {
388  if ( Lambda[0] > Lambda[2] )
389  {
390  ans = Lambda[1] - Lambda[0];
391  }
392  else
393  {
394  if ( Lambda[1] > Lambda[2] )
395  {
396  ans = Lambda[1] - Lambda[2];
397  }
398  else
399  {
400  ans = Lambda[2] - Lambda[1];
401  }
402  }
403  }
404  }
405  else if ( numberOfDistinctRoots == 2 )
406  {
407  if ( Lambda[0] > Lambda[1] )
408  {
409  ans = Lambda[0] - Lambda[1];
410  }
411  else
412  {
413  ans = Lambda[1] - Lambda[0];
414  }
415  }
416  else if ( numberOfDistinctRoots == 1 )
417  {
418  ans = 0.0;
419  }
420  else
421  {
422  itkExceptionMacro(<< "Undefined condition. Cubic root solver returned "
423  << numberOfDistinctRoots << " distinct roots.");
424  }
425 
426  return ans;
427  }
428 
429  // Function is defined here because the templating confuses gcc 2.96 when
430  // defined
431  // in .hxx file. jc 1/29/03
432  TRealType EvaluateAtNeighborhood(const ConstNeighborhoodIteratorType & it) const
433  {
434  unsigned int i, j;
435 
436  vnl_matrix< TRealType > g(ImageDimension, ImageDimension);
437  vnl_vector_fixed< TRealType, VectorDimension >
438  d_phi_du[itk::GetImageDimension < TInputImage > ::ImageDimension];
439 
440  // Calculate the directional derivatives for each vector component using
441  // central differences.
442  for ( i = 0; i < ImageDimension; i++ )
443  {
444  for ( j = 0; j < VectorDimension; j++ )
445  {
446  d_phi_du[i][j] = m_DerivativeWeights[i] * m_SqrtComponentWeights[j]
447  * 0.5 * ( it.GetNext(i)[j] - it.GetPrevious(i)[j] );
448  }
449  }
450 
451  // Calculate the symmetric metric tensor g
452  for ( i = 0; i < ImageDimension; i++ )
453  {
454  for ( j = i; j < ImageDimension; j++ )
455  {
456  g[j][i] = g[i][j] = dot_product(d_phi_du[i], d_phi_du[j]);
457  }
458  }
459 
460  // Find the eigenvalues of g
461  vnl_symmetric_eigensystem< TRealType > E(g);
462 
463  // Return the difference in length between the first two principle axes.
464  // Note that other edge strength metrics may be appropriate here instead..
465  return ( E.get_eigenvalue(ImageDimension - 1) - E.get_eigenvalue(ImageDimension - 2) );
466  }
467 
470 
476 private:
479 
481 
482  typename ImageBaseType::ConstPointer m_RealValuedInputImage;
483 
484  VectorGradientMagnitudeImageFilter(const Self &); //purposely not implemented
485  void operator=(const Self &); //purposely not implemented
486 
488 };
489 } // end namespace itk
490 
491 #ifndef ITK_MANUAL_INSTANTIATION
492 #include "itkVectorGradientMagnitudeImageFilter.hxx"
493 #endif
494 
495 #endif
496