ITK  4.6.0
Insight Segmentation and Registration Toolkit
itkDescoteauxSheetnessImageFilter.h
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: itkDescoteauxSheetnessImageFilter.h
5  Language: C++
6  Date: $Date$
7  Version: $Revision$
8 
9  Copyright (c) Insight Software Consortium. All rights reserved.
10  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 #ifndef __itkDescoteauxSheetnessImageFilter_h
18 #define __itkDescoteauxSheetnessImageFilter_h
19 
21 #include "vnl/vnl_math.h"
22 
23 namespace itk
24 {
25 
40 namespace Function {
41 
42 template< class TInput, class TOutput>
43 class Sheetness
44 {
45 public:
47  {
48  m_Alpha = 0.5; // suggested value in the paper
49  m_Gamma = 0.5; // suggested value in the paper;
50  m_C = 1.0;
51  m_DetectBrightSheets = true;
52  }
54  bool operator!=( const Sheetness & ) const
55  {
56  return false;
57  }
58  bool operator==( const Sheetness & other ) const
59  {
60  return !(*this != other);
61  }
62  inline TOutput operator()( const TInput & A )
63  {
64  double sheetness = 0.0;
65 
66  double a1 = static_cast<double>( A[0] );
67  double a2 = static_cast<double>( A[1] );
68  double a3 = static_cast<double>( A[2] );
69 
70  double l1 = vnl_math_abs( a1 );
71  double l2 = vnl_math_abs( a2 );
72  double l3 = vnl_math_abs( a3 );
73 
74  //
75  // Sort the values by their absolute value.
76  // At the end of the sorting we should have
77  //
78  // l1 <= l2 <= l3
79  //
80  if( l2 > l3 )
81  {
82  double tmpl = l3;
83  l3 = l2;
84  l2 = tmpl;
85  double tmpa = a3;
86  a3 = a2;
87  a2 = tmpa;
88  }
89 
90  if( l1 > l2 )
91  {
92  double tmp = l1;
93  l1 = l2;
94  l2 = tmp;
95  double tmpa = a1;
96  a1 = a2;
97  a2 = tmpa;
98  }
99 
100  if( l2 > l3 )
101  {
102  double tmp = l3;
103  l3 = l2;
104  l2 = tmp;
105  double tmpa = a3;
106  a3 = a2;
107  a2 = tmpa;
108  }
109 
110  if( this->m_DetectBrightSheets )
111  {
112  if( a3 > 0.0 )
113  {
114  return static_cast<TOutput>( sheetness );
115  }
116  }
117  else
118  {
119  if( a3 < 0.0 )
120  {
121  return static_cast<TOutput>( sheetness );
122  }
123  }
124 
125 
126  //
127  // Avoid divisions by zero (or close to zero)
128  //
129  if( static_cast<double>( l3 ) < vnl_math::eps )
130  {
131  return static_cast<TOutput>( sheetness );
132  }
133 
134  const double Rs = l2 / l3;
135  const double Rb = vnl_math_abs( l3 + l3 - l2 - l1 ) / l3;
136  const double Rn = vcl_sqrt( l3*l3 + l2*l2 + l1*l1 );
137 
138  sheetness = vcl_exp( - ( Rs * Rs ) / ( 2.0 * m_Alpha * m_Alpha ) );
139  sheetness *= ( 1.0 - vcl_exp( - ( Rb * Rb ) / ( 2.0 * m_Gamma * m_Gamma ) ) );
140  sheetness *= ( 1.0 - vcl_exp( - ( Rn * Rn ) / ( 2.0 * m_C * m_C ) ) );
141 
142  return static_cast<TOutput>( sheetness );
143  }
144  void SetAlpha( double value )
145  {
146  this->m_Alpha = value;
147  }
148  void SetGamma( double value )
149  {
150  this->m_Gamma = value;
151  }
152  void SetC( double value )
153  {
154  this->m_C = value;
155  }
156  void SetDetectBrightSheets( bool value )
157  {
158  this->m_DetectBrightSheets = value;
159  }
160  void SetDetectDarkSheets( bool value )
161  {
162  this->m_DetectBrightSheets = !value;
163  }
164 
165 private:
166  double m_Alpha;
167  double m_Gamma;
168  double m_C;
170 };
171 }
172 
173 template <class TInputImage, class TOutputImage>
175  public
176 UnaryFunctorImageFilter<TInputImage,TOutputImage,
177  Function::Sheetness< typename TInputImage::PixelType,
178  typename TOutputImage::PixelType> >
179 {
180 public:
183  typedef UnaryFunctorImageFilter<
184  TInputImage,TOutputImage,
186  typename TInputImage::PixelType,
187  typename TOutputImage::PixelType> > Superclass;
190 
192  itkNewMacro(Self);
193 
195  itkTypeMacro(DescoteauxSheetnessImageFilter,
197 
199  void SetSheetnessNormalization( double value )
200  {
201  this->GetFunctor().SetAlpha( value );
202  }
203 
205  void SetBloobinessNormalization( double value )
206  {
207  this->GetFunctor().SetGamma( value );
208  }
209 
211  void SetNoiseNormalization( double value )
212  {
213  this->GetFunctor().SetC( value );
214  }
215  void SetDetectBrightSheets( bool value )
216  {
217  this->GetFunctor().SetDetectBrightSheets( value );
218  }
219  void SetDetectDarkSheets( bool value )
220  {
221  this->GetFunctor().SetDetectDarkSheets( value );
222  }
224 
225 #ifdef ITK_USE_CONCEPT_CHECKING
226 
227  typedef typename TInputImage::PixelType InputPixelType;
228  itkConceptMacro(BracketOperatorsCheck,
230  itkConceptMacro(DoubleConvertibleToOutputCheck,
232 
234 #endif
235 
236 protected:
239 
240 private:
241  DescoteauxSheetnessImageFilter(const Self&); //purposely not implemented
242  void operator=(const Self&); //purposely not implemented
243 
244 };
245 
246 } // end namespace itk
247 
248 
249 #endif
UnaryFunctorImageFilter< TInputImage, TOutputImage, Function::Sheetness< typename TInputImage::PixelType, typename TOutputImage::PixelType > > Superclass
Base class for all process objects that output image data.
Computes a measure of Sheetness from the Hessian Eigenvalues.
Implements pixel-wise generic operation on one image.
bool operator==(const Sheetness &other) const
#define itkConceptMacro(name, concept)
bool operator!=(const Sheetness &) const