Orfeo Toolbox  3.20
itkGaussianDerivativeOperator.txx
Go to the documentation of this file.
1 /*=========================================================================
2
3  Program: Insight Segmentation & Registration Toolkit
4  Module: \$RCSfile: itkGaussianDerivativeOperator.txx,v \$
5  Language: C++
6  Date: \$Date: 2009-04-06 16:49:36 \$
7  Version: \$Revision: 1.5 \$
8
11
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15
16 =========================================================================*/
17
18 #ifndef __itkGaussianDerivativeOperator_txx
19 #define __itkGaussianDerivativeOperator_txx
20
22 #include "itkOutputWindow.h"
23 #include "itkMacro.h"
24
25 namespace itk
26 {
27
28 template<class TPixel,unsigned int VDimension, class TAllocator>
29 typename GaussianDerivativeOperator<TPixel,VDimension, TAllocator>
30 ::CoefficientVector
33 {
34  CoefficientVector coeff;
35
36  // Use image spacing to modify variance
37  m_Variance /= ( m_Spacing * m_Spacing );
38
39  // Calculate normalization factor for derivatives when necessary
40  double norm = m_NormalizeAcrossScale && m_Order ? vcl_pow( m_Variance, m_Order/2.0 ) : 1.0;
41
42  if( !this->GetUseDerivativeOperator() )
43  {
44  // Coefficient of the polynomial that multiplies the gaussian
45  // Gaussian derivatives always take the form
46  // G'(n)(x,t) = P(n)(x,t) * G(x,t)
47  // where P(n)(x,sigma) is a polynomial of the same order as the derivative.
48  // For first order derivatives the polynomial simply corresponds to multiplying
49  // the gaussian by -x/t.
50
51  std::vector<int> polyCoeffs;
52
53  if( m_Order == 1 ) // -x/t
54  {
55  polyCoeffs.push_back(0);
56  polyCoeffs.push_back(-1);
57  }
58  else if( m_Order == 2 ) // ( x^2-t )/t^2
59  {
60  polyCoeffs.push_back(-1);
61  polyCoeffs.push_back(0);
62  polyCoeffs.push_back(1);
63  }
64  else if( m_Order == 3 ) // (-x^3+3xt)/t^3
65  {
66  polyCoeffs.push_back(0);
67  polyCoeffs.push_back(3);
68  polyCoeffs.push_back(0);
69  polyCoeffs.push_back(-1);
70  }
71  else if( m_Order > 3 ) // recursively calculate derivative of polynomial
72  {
73  polyCoeffs.push_back(0);
74  polyCoeffs.push_back(3);
75  polyCoeffs.push_back(0);
76  polyCoeffs.push_back(-1);
77
78  unsigned int i,j;
79  for( i = 4; i <= m_Order; ++i )
80  {
81  if( i%2 == 0 ) // even order
82  {
83  for( j = 1; j < polyCoeffs.size(); j += 2 )
84  {
85  polyCoeffs[j-1] += j*polyCoeffs[j];
86  if( j < polyCoeffs.size()-1 )
87  polyCoeffs[j+1] -= polyCoeffs[j];
88  polyCoeffs[j] = 0;
89  }
90  polyCoeffs.push_back(1); // add highest order new element
91  }
92  else // odd order
93  {
94  polyCoeffs[1] = -polyCoeffs[0];
95  polyCoeffs[0] = 0;
96  for( j = 2; j < polyCoeffs.size(); j += 2 )
97  {
98  polyCoeffs[j-1] += j*polyCoeffs[j];
99  polyCoeffs[j+1] -= polyCoeffs[j];
100  polyCoeffs[j] = 0;
101  }
102  polyCoeffs.push_back(-1); // add highest order new element
103  }
104  }
105  }
106
107  // Now create coefficients as if they were zero order coeffs
108
109  double sum;
110  int i;
111  int j;
112  typename CoefficientVector::iterator it;
113
114  const double et = vcl_exp(-m_Variance);
115  const double cap = 1.0 - m_MaximumError;
116
117  // Create the kernel coefficients as a std::vector
118  sum = 0.0;
119  coeff.push_back(et * ModifiedBesselI0(m_Variance));
120  sum += coeff[0];
121  coeff.push_back(et * ModifiedBesselI1(m_Variance));
122  sum += coeff[1] * 2.0;
123
124  for (i = 2; sum < cap; i++ )
125  {
126  coeff.push_back(et* ModifiedBesselI(i, m_Variance));
127  sum += coeff[i] * 2.0;
128  if (coeff[i] <= 0.0) break; // failsafe
129  if (coeff.size() > m_MaximumKernelWidth )
130  {
131  itkWarningMacro("Kernel size has exceeded the specified maximum width of "
132  << m_MaximumKernelWidth << " and has been truncated to " <<
133  static_cast<unsigned long>( coeff.size() ) << " elements. You can raise "
134  "the maximum width using the SetMaximumKernelWidth method.");
135  break;
136  }
137  }
138
139  // Normalize the coefficients so they sum one
140  for (it = coeff.begin(); it < coeff.end(); ++it)
141  {
142  *it /= sum;
143  }
144
145  // Make symmetric
146  j = static_cast<int>( coeff.size() ) - 1;
147  coeff.insert(coeff.begin(), j, 0);
148  for (i=0, it = coeff.end()-1; i < j; --it, ++i)
149  {
150  coeff[i] = *it;
151  }
152
153  if( m_Order == 0 )
154  return coeff;
155
156  // Now multiply modify the coefficients taking into account the derivative polynomial
157  // and order
158  unsigned int k;
159  for( i=-(int)coeff.size()/2, it = coeff.begin(); i<=(int)coeff.size()/2; ++i )
160  {
161  sum = 0.0;
162  if( m_Order%2 == 0 ) // even
163  {
164  for( j = 0, k = (int)m_Order/2; j < (int)polyCoeffs.size(); j += 2, k-- )
165  {
166  sum += polyCoeffs[j] * vcl_pow( m_Variance,(double)k ) * vcl_pow( i*m_Spacing,(double)j );
167  }
168  }
169  else // odd
170  {
171  for( j = 1, k = (int)(m_Order-1)/2; j < (int)polyCoeffs.size(); j += 2, k-- )
172  {
173  sum += polyCoeffs[j] * vcl_pow( m_Variance,(double)k ) * vcl_pow( i*m_Spacing,(double)j );
174  }
175  }
176  sum *= norm / vcl_pow( m_Variance, static_cast<int>( m_Order ) );
177  (*it) *= sum;
178  ++it;
179  }
180  }
181
182  else // m_UseDerivativeOperator = true
183  {
184  GaussianOperatorType gaussOp;
185  gaussOp.SetDirection( this->GetDirection() );
186  gaussOp.SetMaximumKernelWidth( m_MaximumKernelWidth );
187  gaussOp.SetMaximumError( m_MaximumError );
188  gaussOp.SetVariance( m_Variance / ( m_Spacing * m_Spacing ) );
189  gaussOp.CreateDirectional();
190
191  DerivativeOperatorType derivOp;
192  derivOp.SetDirection( this->GetDirection() );
193  derivOp.SetOrder( m_Order );
194  derivOp.ScaleCoefficients( 1.0 / m_Spacing );
195  derivOp.CreateDirectional();
196
197  // Now perform convolution between both operators
198
199  int i,j,k;
200  double conv;
201
202  for( i=0; i<(int)gaussOp.Size(); ++i ) // current index in gaussian op
203  {
204  conv = 0.0;
205  for( j=0; j<(int)derivOp.Size(); ++j ) // current index in derivative op
206  {
207  k = i + j - derivOp.Size()/2;
208  if( k >= 0 && k < (int)gaussOp.Size() )
209  conv += gaussOp[k] * derivOp[derivOp.Size()-1-j];
210  }
211  coeff.push_back( norm * conv );
212  }
213  }
214
215  return coeff;
216 }
217
218 template<class TPixel,unsigned int VDimension, class TAllocator>
219 double
222 {
223  double d, accumulator;
224  double m;
225
226  if ((d=vcl_fabs(y)) < 3.75)
227  {
228  m=y/3.75;
229  m *= m;
230  accumulator = 1.0 + m *(3.5156229+m*(3.0899424+m*(1.2067492
231  + m*(0.2659732+m*(0.360768e-1 +m*0.45813e-2)))));
232  }
233  else
234  {
235  m=3.5/d;
236  accumulator =(vcl_exp(d)/vcl_sqrt(d))*(0.39894228+m*(0.1328592e-1
237  +m*(0.225319e-2+m*(-0.157565e-2+m*(0.916281e-2
238  +m*(-0.2057706e-1+m*(0.2635537e-1+m*(-0.1647633e-1
239  +m*0.392377e-2))))))));
240  }
241  return accumulator;
242 }
243
244
245 template<class TPixel,unsigned int VDimension, class TAllocator>
246 double
249 {
250  double d, accumulator;
251  double m;
252
253  if ((d=vcl_fabs(y)) < 3.75)
254  {
255  m=y/3.75;
256  m *= m;
257  accumulator = d*(0.5+m*(0.87890594+m*(0.51498869+m*(0.15084934
258  +m*(0.2658733e-1+m*(0.301532e-2+m*0.32411e-3))))));
259  }
260  else
261  {
262  m=3.75/d;
263  accumulator = 0.2282967e-1+m*(-0.2895312e-1+m*(0.1787654e-1
264  -m*0.420059e-2));
265  accumulator = 0.39894228+m*(-0.3988024e-1+m*(-0.362018e-2
266  +m*(0.163801e-2+m*(-0.1031555e-1+m*accumulator))));
267
268  accumulator *= (vcl_exp(d)/vcl_sqrt(d));
269  }
270
271  if (y<0.0) return -accumulator;
272  else return accumulator;
273 }
274
275
276 template<class TPixel,unsigned int VDimension, class TAllocator>
277 double
279 ::ModifiedBesselI(int n, double y)
280 {
281  const double ACCURACY = 40.0;
282  int j;
283  double qim, qi, qip, toy;
284  double accumulator;
285
286  if (n<2)
287  {
288  throw ExceptionObject(__FILE__, __LINE__, "Order of modified bessel is > 2.", ITK_LOCATION); // placeholder
289  }
290  if (y==0.0) return 0.0;
291  else
292  {
293  toy=2.0/vcl_fabs(y);
294  qip=accumulator=0.0;
295  qi=1.0;
296  for (j=2*(n+(int)vcl_sqrt(ACCURACY*n)); j>0; j--)
297  {
298  qim=qip+j*toy*qi;
299  qip=qi;
300  qi=qim;
301  if (vcl_fabs(qi) > 1.0e10)
302  {
303  accumulator *= 1.0e-10;
304  qi *= 1.0e-10;
305  qip *= 1.0e-10;
306  }
307  if (j==n) accumulator=qip;
308  }
309  accumulator *= ModifiedBesselI0(y)/qi;
310  if (y<0.0 && (n&1)) return -accumulator;
311  else return accumulator;
312  }
313 }
314
315 }// end namespace itk
316
317 #endif

Generated at Sat Nov 9 2013 14:30:48 for Orfeo Toolbox with doxygen 1.8.3.1