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 itkGaussianDerivativeOperator_hxx
19 #define itkGaussianDerivativeOperator_hxx
20
21 #include "itkGaussianDerivativeOperator.h"
22 #include "itkCompensatedSummation.h"
23 #include "itkOutputWindow.h"
24 #include "itkMacro.h"
25 #include <numeric>
26
27 namespace itk
28 {
29
30 template< typename TPixel, unsigned int VDimension, typename TAllocator >
31 GaussianDerivativeOperator< TPixel, VDimension, TAllocator >
GaussianDerivativeOperator(const Self & other)32 ::GaussianDerivativeOperator(const Self & other)
33 : NeighborhoodOperator< TPixel, VDimension, TAllocator >(other),
34 m_NormalizeAcrossScale(other.m_NormalizeAcrossScale),
35 m_Variance(other.m_Variance),
36 m_MaximumError(other.m_MaximumError),
37 m_MaximumKernelWidth(other.m_MaximumKernelWidth),
38 m_Order(other.m_Order),
39 m_Spacing(other.m_Spacing)
40 {
41 }
42
43 template< typename TPixel, unsigned int VDimension, typename TAllocator >
44 GaussianDerivativeOperator< TPixel, VDimension, TAllocator > &
45 GaussianDerivativeOperator< TPixel, VDimension, TAllocator >
operator =(const Self & other)46 ::operator=(const Self & other)
47 {
48 if(this != &other)
49 {
50 Superclass::operator=(other);
51 m_NormalizeAcrossScale = other.m_NormalizeAcrossScale;
52 m_Spacing = other.m_Spacing;
53 m_Order = other.m_Order;
54 m_Variance = other.m_Variance;
55 m_MaximumError = other.m_MaximumError;
56 m_MaximumKernelWidth = other.m_MaximumKernelWidth;
57 }
58 return *this;
59 }
60
61 template< typename TPixel, unsigned int VDimension, typename TAllocator >
62 typename GaussianDerivativeOperator< TPixel, VDimension, TAllocator >
63 ::CoefficientVector
64 GaussianDerivativeOperator< TPixel, VDimension, TAllocator >
GenerateCoefficients()65 ::GenerateCoefficients()
66 {
67
68 // compute gaussian kernel of 0-order
69 CoefficientVector coeff = this->GenerateGaussianCoefficients();
70
71 if ( m_Order == 0 )
72 {
73 return coeff;
74 }
75
76
77 // Calculate scale-space normalization factor for derivatives
78 double norm;
79 if (m_NormalizeAcrossScale && m_Order)
80 {
81 norm = std::pow(m_Variance, m_Order / 2.0);
82 }
83 else
84 {
85 norm = 1.0;
86 }
87
88 // additional normalization for spacing
89 norm /= std::pow( m_Spacing, static_cast< int >( m_Order ) );
90
91 DerivativeOperatorType derivOp;
92 derivOp.SetDirection( this->GetDirection() );
93 derivOp.SetOrder( m_Order );
94 derivOp.CreateDirectional();
95
96 // The input gaussian kernel needs to be padded with a clamped
97 // boundary condition. If N is the radius of the derivative
98 // operator, then the output kernel needs to be padded by N-1. For
99 // these values to be computed the input kernel needs to be padded
100 // by 2N-1 on both sides.
101 unsigned int N = ( derivOp.Size() - 1 ) / 2;
102
103 // copy the gaussian operator adding clamped boundary condition
104 CoefficientVector paddedCoeff( coeff.size() + 4*N - 2);
105
106 // copy the whole gaussuan operator in coeff to paddedCoef
107 // starting after the padding
108 std::copy( coeff.begin(), coeff.end(), paddedCoeff.begin() + 2*N - 1);
109
110 // padd paddedCoeff with 2*N-1 number of boundary conditions
111 std::fill( paddedCoeff.begin(), paddedCoeff.begin() + 2*N, coeff.front() );
112 std::fill( paddedCoeff.rbegin(), paddedCoeff.rbegin() + 2*N, coeff.back() );
113
114 // clear for output kernel/coeffs
115 coeff = CoefficientVector();
116
117 // Now perform convolution between derivative operators and padded gaussian
118 for ( unsigned int i = N; i < paddedCoeff.size() - N; ++i )
119 {
120 CompensatedSummation<double> conv;
121
122 // current index in derivative op
123 for ( unsigned int j = 0; j < derivOp.Size(); ++j )
124 {
125 unsigned int k = i + j - derivOp.Size() / 2;
126 conv += paddedCoeff[k] * derivOp[derivOp.Size() - 1 - j];
127 }
128
129 // normalize for scale-space and spacing
130 coeff.push_back(norm * conv.GetSum());
131 }
132
133 return coeff;
134 }
135
136 template< typename TPixel, unsigned int VDimension, typename TAllocator >
137 typename GaussianDerivativeOperator< TPixel, VDimension, TAllocator >::CoefficientVector
138 GaussianDerivativeOperator< TPixel, VDimension, TAllocator >
GenerateGaussianCoefficients() const139 ::GenerateGaussianCoefficients() const
140 {
141
142 CoefficientVector coeff;
143
144 // Use image spacing to modify variance
145 const double pixelVariance = m_Variance / ( m_Spacing * m_Spacing );
146
147 // Now create coefficients as if they were zero order coeffs
148 const double et = std::exp(-pixelVariance);
149 const double cap = 1.0 - m_MaximumError;
150 CompensatedSummation<double> sum;
151
152 // Create the kernel coefficients as a std::vector
153 coeff.push_back( et * ModifiedBesselI0(pixelVariance) );
154 sum += coeff[0];
155 coeff.push_back( et * ModifiedBesselI1(pixelVariance) );
156 sum += coeff[1] * 2.0;
157
158 for ( int i = 2; sum.GetSum() < cap; i++ )
159 {
160 coeff.push_back( et * ModifiedBesselI(i, pixelVariance) );
161 sum += coeff[i] * 2.0;
162 if ( coeff[i] < sum.GetSum()*NumericTraits<double>::epsilon() )
163 {
164 // if the coeff is less then this value then the value of cap
165 // will not change, and it's will not contribute to the operator
166 itkWarningMacro( "Kernel failed to accumulate to approximately one with current remainder "
167 << cap-sum.GetSum() << " and current coefficient " << coeff[i] << "." );
168
169 break;
170 }
171 if ( coeff.size() > m_MaximumKernelWidth )
172 {
173 itkWarningMacro("Kernel size has exceeded the specified maximum width of "
174 << m_MaximumKernelWidth << " and has been truncated to "
175 << static_cast< unsigned long >( coeff.size() ) << " elements. You can raise "
176 "the maximum width using the SetMaximumKernelWidth method.");
177 break;
178 }
179 }
180
181 // re-accumulate from smallest number to largest for maximum precision
182 sum = std::accumulate( coeff.rbegin(), coeff.rend() - 1, 0.0 );
183 sum *= 2.0;
184 sum += coeff[0]; // the first is only needed once
185
186 // Normalize the coefficients so they sum one
187 for ( auto it = coeff.begin(); it != coeff.end(); ++it )
188 {
189 *it /= sum.GetSum();
190 }
191
192 // Make symmetric
193 size_t s = coeff.size() - 1;
194 coeff.insert(coeff.begin(), s, 0);
195 std::copy( coeff.rbegin(), coeff.rbegin() + s, coeff.begin() );
196
197 return coeff;
198 }
199
200 template< typename TPixel, unsigned int VDimension, typename TAllocator >
201 double
202 GaussianDerivativeOperator< TPixel, VDimension, TAllocator >
ModifiedBesselI0(double y)203 ::ModifiedBesselI0(double y)
204 {
205 double d, accumulator;
206 double m;
207
208 if ( ( d = std::fabs(y) ) < 3.75 )
209 {
210 m = y / 3.75;
211 m *= m;
212 accumulator = 1.0 + m * ( 3.5156229 + m * ( 3.0899424 + m * ( 1.2067492
213 + m
214 * ( 0.2659732 + m * ( 0.360768e-1 + m * 0.45813e-2 ) ) ) ) );
215 }
216 else
217 {
218 m = 3.75 / d;
219 accumulator = ( std::exp(d) / std::sqrt(d) ) * ( 0.39894228 + m * ( 0.1328592e-1
220 + m
221 * ( 0.225319e-2 + m
222 * ( -0.157565e-2 + m * ( 0.916281e-2
223 +
224 m
225 * ( -0.2057706e-1
226 + m *
227 ( 0.2635537e-1 +
228 m *
229 ( -0.1647633e-1
230 +
231 m
232 *
233 0.392377e-2 ) ) ) ) ) ) ) );
234 }
235 return accumulator;
236 }
237
238 template< typename TPixel, unsigned int VDimension, typename TAllocator >
239 double
240 GaussianDerivativeOperator< TPixel, VDimension, TAllocator >
ModifiedBesselI1(double y)241 ::ModifiedBesselI1(double y)
242 {
243 double d, accumulator;
244 double m;
245
246 if ( ( d = std::fabs(y) ) < 3.75 )
247 {
248 m = y / 3.75;
249 m *= m;
250 accumulator = d * ( 0.5 + m * ( 0.87890594 + m * ( 0.51498869 + m * ( 0.15084934
251 + m
252 * ( 0.2658733e-1 + m
253 * ( 0.301532e-2 + m * 0.32411e-3 ) ) ) ) ) );
254 }
255 else
256 {
257 m = 3.75 / d;
258 accumulator = 0.2282967e-1 + m * ( -0.2895312e-1 + m * ( 0.1787654e-1
259 - m * 0.420059e-2 ) );
260 accumulator = 0.39894228 + m * ( -0.3988024e-1 + m * ( -0.362018e-2
261 + m * ( 0.163801e-2 + m * ( -0.1031555e-1 + m * accumulator ) ) ) );
262
263 accumulator *= ( std::exp(d) / std::sqrt(d) );
264 }
265
266 if ( y < 0.0 ) { return -accumulator; }
267 else { return accumulator; }
268 }
269
270 template< typename TPixel, unsigned int VDimension, typename TAllocator >
271 double
272 GaussianDerivativeOperator< TPixel, VDimension, TAllocator >
ModifiedBesselI(int n,double y)273 ::ModifiedBesselI(int n, double y)
274 {
275 constexpr double DIGITS = 10.0;
276 int j;
277 double qim, qi, qip, toy;
278 double accumulator;
279
280 if ( n < 2 )
281 {
282 throw ExceptionObject(__FILE__, __LINE__, "Order of modified bessel is > 2.", ITK_LOCATION); //
283 // placeholder
284 }
285 if ( y == 0.0 ) { return 0.0; }
286 else
287 {
288 toy = 2.0 / std::fabs(y);
289 qip = accumulator = 0.0;
290 qi = 1.0;
291 for ( j = 2 * ( n + (int)(DIGITS*std::sqrt((double)n) ) ); j > 0; j-- )
292 {
293 qim = qip + j * toy * qi;
294 qip = qi;
295 qi = qim;
296 if ( std::fabs(qi) > 1.0e10 )
297 {
298 accumulator *= 1.0e-10;
299 qi *= 1.0e-10;
300 qip *= 1.0e-10;
301 }
302 if ( j == n ) { accumulator = qip; }
303 }
304 accumulator *= ModifiedBesselI0(y) / qi;
305 if ( y < 0.0 && ( n & 1 ) ) { return -accumulator; }
306 else { return accumulator; }
307 }
308 }
309
310 template< typename TPixel, unsigned int VDimension, typename TAllocator >
311 void
312 GaussianDerivativeOperator< TPixel, VDimension, TAllocator >
PrintSelf(std::ostream & os,Indent i) const313 ::PrintSelf(std::ostream & os, Indent i) const
314 {
315 os << i << "GaussianDerivativeOperator { this=" << this
316 << ", m_NormalizeAcrossScale = " << m_NormalizeAcrossScale
317 << ", m_Order = " << m_Order
318 << ", m_Spacing = " << m_Spacing
319 << ", m_Variance = " << m_Variance
320 << ", m_MaximumError = " << m_MaximumError
321 << ", m_MaximumKernelWidth = " << m_MaximumKernelWidth
322 << "} " << std::endl;
323 Superclass::PrintSelf( os, i.GetNextIndent() );
324 }
325
326 } // end namespace itk
327
328 #endif
329