| 1 |
|
/*========================================================================= |
| 2 |
|
|
| 3 |
|
Program: Insight Segmentation & Registration Toolkit |
| 4 |
|
Module: $RCSfile: itkGaussianOperator.txx.html,v $ |
| 5 |
|
Language: C++ |
| 6 |
|
Date: $Date: 2006/01/17 19:15:36 $ |
| 7 |
|
Version: $Revision: 1.4 $ |
| 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 |
DEF |
#ifndef _itkGaussianOperator_txx |
| 18 |
DEF |
#define _itkGaussianOperator_txx |
| 19 |
|
#include "itkGaussianOperator.h" |
| 20 |
|
#include "itkOutputWindow.h" |
| 21 |
|
#include "itkMacro.h" |
| 22 |
|
namespace itk |
| 23 |
|
{ |
| 24 |
|
|
| 25 |
|
template<class TPixel,unsigned int VDimension, class TAllocator> |
| 26 |
|
typename GaussianOperator<TPixel,VDimension, TAllocator> |
| 27 |
|
::CoefficientVector |
| 28 |
|
GaussianOperator<TPixel,VDimension, TAllocator> |
| 29 |
|
::GenerateCoefficients() |
| 30 |
|
{ |
| 31 |
|
CoefficientVector coeff; |
| 32 |
|
double sum; |
| 33 |
|
int i; |
| 34 |
|
int j; |
| 35 |
|
typename CoefficientVector::iterator it; |
| 36 |
|
|
| 37 |
|
const double et = ::exp(-m_Variance); |
| 38 |
|
const double cap = 1.0 - m_MaximumError; |
| 39 |
|
|
| 40 |
|
// Create the kernel coefficients as a std::vector |
| 41 |
|
sum = 0.0; |
| 42 |
|
coeff.push_back(et * ModifiedBesselI0(m_Variance)); |
| 43 |
|
sum += coeff[0]; |
| 44 |
|
coeff.push_back(et * ModifiedBesselI1(m_Variance)); |
| 45 |
|
sum += coeff[1]* 2.0; |
| 46 |
|
|
| 47 |
|
for (i=2; sum < cap; i++) |
| 48 |
|
{ |
| 49 |
|
coeff.push_back(et* ModifiedBesselI(i, m_Variance)); |
| 50 |
|
sum += coeff[i] *2.0; |
| 51 |
|
if (coeff[i] <= 0.0) break; // failsafe |
| 52 |
|
if (coeff.size() > m_MaximumKernelWidth ) |
| 53 |
|
{ |
| 54 |
LEN |
itkWarningMacro("Kernel size has exceeded the specified maximum width of " << m_MaximumKernelWidth << " and has been truncated to " << static_cast<unsigned long>( coeff.size() ) << " elements. You can raise the maximum width using the SetMaximumKernelWidth method."); |
| 55 |
|
break; |
| 56 |
|
} |
| 57 |
|
} |
| 58 |
|
// Normalize the coefficients so that their sum is one. |
| 59 |
|
for (it = coeff.begin(); it < coeff.end(); ++it) |
| 60 |
|
{ |
| 61 |
|
*it /= sum; |
| 62 |
|
} |
| 63 |
|
|
| 64 |
|
// Make symmetric |
| 65 |
|
j = static_cast<int>( coeff.size() ) - 1; |
| 66 |
|
coeff.insert(coeff.begin(), j, 0); |
| 67 |
|
for (i=0, it = coeff.end()-1; i < j; --it, ++i) |
| 68 |
|
{ |
| 69 |
|
coeff[i] = *it; |
| 70 |
|
} |
| 71 |
|
|
| 72 |
|
return coeff; |
| 73 |
|
} |
| 74 |
|
|
| 75 |
|
template<class TPixel,unsigned int VDimension, class TAllocator> |
| 76 |
|
double |
| 77 |
|
GaussianOperator<TPixel,VDimension, TAllocator> |
| 78 |
|
::ModifiedBesselI0(double y) |
| 79 |
|
{ |
| 80 |
|
double d, accumulator; |
| 81 |
|
double m; |
| 82 |
|
|
| 83 |
|
if ((d=fabs(y)) < 3.75) |
| 84 |
|
{ |
| 85 |
|
m=y/3.75; |
| 86 |
|
m*=m; |
| 87 |
|
accumulator = 1.0 + m *(3.5156229+m*(3.0899424+m*(1.2067492 |
| 88 |
LEN |
+ m*(0.2659732+m*(0.360768e-1 +m*0.45813e-2))))); |
| 89 |
|
} |
| 90 |
|
else |
| 91 |
|
{ |
| 92 |
|
m=3.5/d; |
| 93 |
|
accumulator =(::exp(d)/::sqrt(d))*(0.39894228+m*(0.1328592e-1 |
| 94 |
LEN |
+m*(0.225319e-2+m*(-0.157565e-2+m*(0.916281e-2 |
| 95 |
LEN |
+m*(-0.2057706e-1+m*(0.2635537e-1+m*(-0.1647633e-1 |
| 96 |
LEN |
+m*0.392377e-2)))))))); |
| 97 |
|
} |
| 98 |
|
return accumulator; |
| 99 |
|
} |
| 100 |
|
|
| 101 |
|
|
| 102 |
|
template<class TPixel,unsigned int VDimension, class TAllocator> |
| 103 |
|
double |
| 104 |
|
GaussianOperator<TPixel,VDimension, TAllocator> |
| 105 |
|
::ModifiedBesselI1(double y) |
| 106 |
|
{ |
| 107 |
|
double d, accumulator; |
| 108 |
|
double m; |
| 109 |
|
|
| 110 |
|
if ((d=fabs(y)) < 3.75) |
| 111 |
|
{ |
| 112 |
|
m=y/3.75; |
| 113 |
|
m*=m; |
| 114 |
|
accumulator = d*(0.5+m*(0.87890594+m*(0.51498869+m*(0.15084934 |
| 115 |
LEN |
+m*(0.2658733e-1+m*(0.301532e-2+m*0.32411e-3)))))); |
| 116 |
|
} |
| 117 |
|
else |
| 118 |
|
{ |
| 119 |
|
m=3.75/d; |
| 120 |
|
accumulator = 0.2282967e-1+m*(-0.2895312e-1+m*(0.1787654e-1 |
| 121 |
|
-m*0.420059e-2)); |
| 122 |
|
accumulator = 0.39894228+m*(-0.3988024e-1+m*(-0.362018e-2 |
| 123 |
LEN |
+m*(0.163801e-2+m*(-0.1031555e-1+m*accumulator)))); |
| 124 |
|
|
| 125 |
|
accumulator *= (::exp(d)/::sqrt(d)); |
| 126 |
|
} |
| 127 |
|
|
| 128 |
|
if (y<0.0) return -accumulator; |
| 129 |
|
else return accumulator; |
| 130 |
|
|
| 131 |
|
} |
| 132 |
|
|
| 133 |
|
template<class TPixel,unsigned int VDimension, class TAllocator> |
| 134 |
|
double |
| 135 |
|
GaussianOperator<TPixel,VDimension, TAllocator> |
| 136 |
|
::ModifiedBesselI(int n, double y) |
| 137 |
|
{ |
| 138 |
|
const double ACCURACY = 40.0; |
| 139 |
|
int j; |
| 140 |
|
double qim, qi, qip, toy; |
| 141 |
|
double accumulator; |
| 142 |
|
|
| 143 |
|
if (n<2) |
| 144 |
|
{ |
| 145 |
LEN |
throw ExceptionObject(__FILE__, __LINE__, "Order of modified bessel is > 2.", ITK_LOCATION); // placeholder |
| 146 |
|
} |
| 147 |
|
if (y==0.0) return 0.0; |
| 148 |
|
else |
| 149 |
|
{ |
| 150 |
|
toy=2.0/fabs(y); |
| 151 |
|
qip=accumulator=0.0; |
| 152 |
|
qi=1.0; |
| 153 |
SEM |
for (j=2*(n+(int)::sqrt(ACCURACY*n)); j>0 ; j--) |
| 154 |
|
{ |
| 155 |
|
qim=qip+j*toy*qi; |
| 156 |
|
qip=qi; |
| 157 |
|
qi=qim; |
| 158 |
|
if (fabs(qi) > 1.0e10) |
| 159 |
|
{ |
| 160 |
|
accumulator*=1.0e-10; |
| 161 |
|
qi *=1.0e-10; |
| 162 |
|
qip*=1.0e-10; |
| 163 |
|
} |
| 164 |
|
if (j==n) accumulator=qip; |
| 165 |
|
} |
| 166 |
|
accumulator *= ModifiedBesselI0(y)/qi; |
| 167 |
|
if (y<0.0 && (n&1)) return -accumulator; |
| 168 |
|
else return accumulator; |
| 169 |
|
} |
| 170 |
|
} |
| 171 |
|
|
| 172 |
|
}// end namespace itk |
| 173 |
|
|
| 174 |
|
#endif |
| 175 |
|
|
