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 /*=========================================================================
19 *
20 * Portions of this file are subject to the VTK Toolkit Version 3 copyright.
21 *
22 * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
23 *
24 * For complete copyright, license and disclaimer of warranty information
25 * please refer to the NOTICE file at the top of the ITK source tree.
26 *
27 *=========================================================================*/
28 #ifndef itkBSplineCenteredResampleImageFilterBase_hxx
29 #define itkBSplineCenteredResampleImageFilterBase_hxx
30 #include "itkBSplineCenteredResampleImageFilterBase.h"
31 #include "itkImageLinearIteratorWithIndex.h"
32 #include "itkProgressReporter.h"
33 namespace itk
34 {
35 /**
36 * Standard "PrintSelf" method
37 */
38 template< typename TInputImage, typename TOutputImage >
39 void
40 BSplineCenteredResampleImageFilterBase< TInputImage, TOutputImage >
PrintSelf(std::ostream & os,Indent indent) const41 ::PrintSelf(
42 std::ostream & os,
43 Indent indent) const
44 {
45 Superclass::PrintSelf(os, indent);
46 }
47
48 template< typename TInputImage, typename TOutputImage >
49 void BSplineCenteredResampleImageFilterBase< TInputImage, TOutputImage >
InitializePyramidSplineFilter(int SplineOrder)50 ::InitializePyramidSplineFilter(int SplineOrder)
51 {
52 switch ( SplineOrder )
53 {
54 case 1:
55 this->m_GSize = 21;
56 this->m_HSize = 2;
57 this->m_G.resize(this->m_GSize);
58 this->m_H.resize(this->m_HSize);
59 this->m_G[0] = 1.;
60 this->m_G[1] = 0.333333;
61 this->m_G[2] = -0.333333;
62 this->m_G[3] = -0.111111;
63 this->m_G[4] = 0.111111;
64 this->m_G[5] = 0.037037;
65 this->m_G[6] = -0.037037;
66 this->m_G[7] = -0.0123457;
67 this->m_G[8] = 0.0123457;
68 this->m_G[9] = 0.00411523;
69 this->m_G[10] = -0.00411523;
70 this->m_G[11] = -0.00137174;
71 this->m_G[12] = 0.00137174;
72 this->m_G[13] = 0.000457247;
73 this->m_G[14] = -0.000457247;
74 this->m_G[15] = -0.000152416;
75 this->m_G[16] = 0.000152416;
76 this->m_G[17] = 0.0000508053;
77 this->m_G[18] = -0.0000508053;
78 this->m_G[19] = -0.0000169351;
79 this->m_G[20] = 0.0000169351;
80 this->m_H[0] = 1;
81 this->m_H[1] = 0.5;
82 break;
83
84 case 2:
85 this->m_GSize = 21;
86 this->m_HSize = 11;
87 this->m_G.resize(this->m_GSize);
88 this->m_H.resize(this->m_HSize);
89 this->m_G[0] = 0.738417;
90 this->m_G[1] = 0.307916;
91 this->m_G[2] = -0.171064;
92 this->m_G[3] = -0.0799199;
93 this->m_G[4] = 0.0735791;
94 this->m_G[5] = 0.03108;
95 this->m_G[6] = -0.0307862;
96 this->m_G[7] = -0.0128561;
97 this->m_G[8] = 0.0128425;
98 this->m_G[9] = 0.00535611;
99 this->m_G[10] = -0.00535548;
100 this->m_G[11] = -0.00223325;
101 this->m_G[12] = 0.00223322;
102 this->m_G[13] = 0.000931242;
103 this->m_G[14] = -0.00093124;
104 this->m_G[15] = -0.000388322;
105 this->m_G[16] = 0.000388322;
106 this->m_G[17] = 0.000161928;
107 this->m_G[18] = -0.000161928;
108 this->m_G[19] = -0.0000675233;
109 this->m_G[20] = 0.0000675233;
110 this->m_H[0] = 1.20711;
111 this->m_H[1] = 0.585786;
112 this->m_H[2] = -0.12132;
113 this->m_H[3] = -0.100505;
114 this->m_H[4] = 0.0208153;
115 this->m_H[5] = 0.0172439;
116 this->m_H[6] = -0.00357134;
117 this->m_H[7] = -0.00295859;
118 this->m_H[8] = 0.000612745;
119 this->m_H[9] = 0.000507614;
120 this->m_H[10] = -0.00010513;
121 break;
122
123 case 3:
124 this->m_GSize = 21;
125 this->m_HSize = 16;
126 this->m_G.resize(this->m_GSize);
127 this->m_H.resize(this->m_HSize);
128 this->m_G[0] = 0.708792;
129 this->m_G[1] = 0.328616;
130 this->m_G[2] = -0.165157;
131 this->m_G[3] = -0.114448;
132 this->m_G[4] = 0.0944036;
133 this->m_G[5] = 0.0543881;
134 this->m_G[6] = -0.05193;
135 this->m_G[7] = -0.0284868;
136 this->m_G[8] = 0.0281854;
137 this->m_G[9] = 0.0152877;
138 this->m_G[10] = -0.0152508;
139 this->m_G[11] = -0.00825077;
140 this->m_G[12] = 0.00824629;
141 this->m_G[13] = 0.00445865;
142 this->m_G[14] = -0.0044582;
143 this->m_G[15] = -0.00241009;
144 this->m_G[16] = 0.00241022;
145 this->m_G[17] = 0.00130278;
146 this->m_G[18] = -0.00130313;
147 this->m_G[19] = -0.000704109;
148 this->m_G[20] = 0.000704784;
149 this->m_H[0] = 1.13726;
150 this->m_H[1] = 0.625601;
151 this->m_H[2] = -0.0870191;
152 this->m_H[3] = -0.159256;
153 this->m_H[4] = 0.0233167;
154 this->m_H[5] = 0.0426725;
155 this->m_H[6] = -0.00624769;
156 this->m_H[7] = -0.0114341;
157 this->m_H[8] = 0.00167406;
158 this->m_H[9] = 0.00306375;
159 this->m_H[10] = -0.000448564;
160 this->m_H[11] = -0.000820929;
161 this->m_H[12] = 0.000120192;
162 this->m_H[13] = 0.000219967;
163 this->m_H[14] = -0.0000322054;
164 this->m_H[15] = -0.00005894;
165 break;
166
167 case 4:
168 this->m_GSize = 21;
169 this->m_HSize = 20;
170 this->m_G.resize(this->m_GSize);
171 this->m_H.resize(this->m_HSize);
172 this->m_G[0] = 0.673072;
173 this->m_G[1] = 0.331218;
174 this->m_G[2] = -0.139359;
175 this->m_G[3] = -0.12051;
176 this->m_G[4] = 0.086389;
177 this->m_G[5] = 0.0611801;
178 this->m_G[6] = -0.0542989;
179 this->m_G[7] = -0.034777;
180 this->m_G[8] = 0.033388;
181 this->m_G[9] = 0.0206275;
182 this->m_G[10] = -0.0203475;
183 this->m_G[11] = -0.0124183;
184 this->m_G[12] = 0.0123625;
185 this->m_G[13] = 0.00751369;
186 this->m_G[14] = -0.00750374;
187 this->m_G[15] = -0.00455348;
188 this->m_G[16] = 0.00455363;
189 this->m_G[17] = 0.00276047;
190 this->m_G[18] = -0.00276406;
191 this->m_G[19] = -0.00167279;
192 this->m_G[20] = 0.00167938;
193 this->m_H[0] = 1.14324;
194 this->m_H[1] = 0.643609;
195 this->m_H[2] = -0.0937888;
196 this->m_H[3] = -0.194993;
197 this->m_H[4] = 0.030127;
198 this->m_H[5] = 0.0699433;
199 this->m_H[6] = -0.0108345;
200 this->m_H[7] = -0.0252663;
201 this->m_H[8] = 0.00391424;
202 this->m_H[9] = 0.00912967;
203 this->m_H[10] = -0.00141437;
204 this->m_H[11] = -0.00329892;
205 this->m_H[12] = 0.000511068;
206 this->m_H[13] = 0.00119204;
207 this->m_H[14] = -0.00018467;
208 this->m_H[15] = -0.000430732;
209 this->m_H[16] = 0.0000667289;
210 this->m_H[17] = 0.000155641;
211 this->m_H[18] = -0.0000241119;
212 this->m_H[19] = -0.0000562395;
213 break;
214 default:
215 // Throw an execption for unsupported splines.
216 ExceptionObject err(__FILE__, __LINE__);
217 err.SetLocation(ITK_LOCATION);
218 err.SetDescription(
219 "SplineOrder for Centered pyramid filter must be 0 through 4. Requested spline order has not been implemented.");
220 throw err;
221 break;
222 }
223 }
224
225 /** Reduce1DImage - reduces the vector of data (in) by a
226 * factor of 2 and writes the results to the location specified
227 * by the Iterator (out). inTraverseSize is the size of the in vector.
228 */
229 template< typename TInputImage, typename TOutputImage >
230 void BSplineCenteredResampleImageFilterBase< TInputImage, TOutputImage >
Reduce1DImage(const std::vector<double> & in,OutputImageIterator & out,unsigned int inTraverseSize,ProgressReporter & progress)231 ::Reduce1DImage(const std::vector< double > & in, OutputImageIterator & out,
232 unsigned int inTraverseSize, ProgressReporter & progress)
233 {
234 IndexValueType i1, i2;
235
236 SizeValueType outK, inK;
237 SizeValueType outTraverseSize = inTraverseSize / 2;
238
239 inTraverseSize = outTraverseSize * 2; // ensures that an even number is used.
240 SizeValueType inModK; // number for modulus math of in
241 inModK = 2L * inTraverseSize;
242
243 // TODO: Need to allocate this once as a scratch variable instead of each
244 // time through.
245 std::vector< double > temp;
246 temp.resize(inTraverseSize);
247
248 for ( inK = 0; inK < inTraverseSize; inK++ )
249 {
250 temp[inK] = in[inK] * this->m_G[0];
251
252 for ( int i = 1; i < this->m_GSize; i++ )
253 {
254 // Calculate indices for left and right of symmetrical filter.
255 i1 = inK - i;
256 i2 = inK + i;
257 // reflect at boundaries if necessary
258 if ( i1 < 0 )
259 {
260 i1 = ( 2L * inTraverseSize - 1L - i1 ) % inModK;
261 if ( i1 >= (IndexValueType)inTraverseSize )
262 {
263 i1 = inModK - i1 - 1L;
264 }
265 }
266 if ( i2 >= (IndexValueType)inTraverseSize ) // originally (i2 > (inTraverseSize - 1)
267 // )
268 {
269 i2 = i2 % inModK;
270 if ( i2 >= (IndexValueType)inTraverseSize )
271 {
272 i2 = inModK - i2 - 1L;
273 }
274 }
275 temp[inK] += this->m_G[i] * ( in[i1] + in[i2] );
276 }
277 }
278
279 for ( outK = 0; outK < outTraverseSize; outK++ )
280 {
281 i1 = 2 * outK;
282 double outVal = ( temp[i1] + temp[i1 + 1] ) / 2.0;
283 out.Set( static_cast< OutputImagePixelType >( outVal ) );
284 ++out;
285 progress.CompletedPixel();
286 }
287 }
288
289 /* Expand1DImage - expands the vector of data (in) by a
290 * factor of 2 and writes the results to the location specified
291 * by the Iterator (out). inTraverseSize is the size of the in vector.
292 */
293 template< typename TInputImage, typename TOutputImage >
294 void BSplineCenteredResampleImageFilterBase< TInputImage, TOutputImage >
Expand1DImage(const std::vector<double> & in,OutputImageIterator & out,unsigned int inTraverseSize,ProgressReporter & progress)295 ::Expand1DImage(const std::vector< double > & in, OutputImageIterator & out,
296 unsigned int inTraverseSize, ProgressReporter & progress)
297 {
298 IndexValueType i1, i2;
299
300 IndexValueType inK;
301 SizeValueType outTraverseSize = inTraverseSize * 2;
302 //inTraverseSize = outTraverseSize/2; // ensures that an even number is used.
303 IndexValueType inModK; // number for modulus math of in
304
305 inModK = outTraverseSize;
306 IndexValueType k0 = ( this->m_HSize / 2 ) * 2 - 1L;
307
308 double outVal, outVal2;
309
310 for ( inK = 0; inK < (IndexValueType)inTraverseSize; inK++ )
311 {
312 //outK = inK * 2L;
313 outVal = in[inK] * this->m_H[0];
314 for ( int k = 2; k < this->m_HSize; k += 2 )
315 {
316 i1 = inK - k / 2L;
317 i2 = inK + k / 2L;
318 if ( i1 < 0 )
319 { // provide
320 // correct border
321 // condition
322 i1 = ( 2L * (IndexValueType)inTraverseSize - 1L - i1 ) % inModK; // pseudo mirror
323 // image
324 if ( i1 >= (IndexValueType)inTraverseSize )
325 {
326 i1 = outTraverseSize - i1 - 1L;
327 }
328 }
329 if ( i2 >= (IndexValueType)inTraverseSize )
330 {
331 i2 = i2 % inModK;
332 if ( i2 >= (IndexValueType)inTraverseSize )
333 {
334 i2 = outTraverseSize - i2 - 1L;
335 }
336 }
337 outVal += this->m_H[k] * ( in[i1] + in[i2] );
338 }
339 out.Set( static_cast< OutputImagePixelType >( outVal ) );
340 ++out;
341 outVal2 = 0;
342 for ( IndexValueType k = -k0; k < this->m_HSize; k += 2L )
343 {
344 IndexValueType kk = std::abs( static_cast< int >( k ) );
345 i1 = inK + ( k + 1L ) / 2L;
346 if ( i1 < 0L )
347 {
348 i1 = ( 2 * inTraverseSize - 1 - i1 ) % inModK;
349 if ( i1 > (IndexValueType)inTraverseSize - 1L )
350 {
351 i1 = outTraverseSize - i1 - 1L;
352 }
353 }
354 if ( i1 >= (IndexValueType)inTraverseSize )
355 {
356 i1 = i1 % inModK;
357 if ( i1 >= (IndexValueType)inTraverseSize )
358 {
359 i1 = outTraverseSize - i1 - 1;
360 }
361 }
362 outVal2 += this->m_H[kk] * in[i1];
363 }
364 out.Set( static_cast< OutputImagePixelType >( outVal2 ) );
365 ++out;
366 }
367
368 // Now apply the Haar[-x]
369 --out;
370 for ( IndexValueType j = outTraverseSize - 1; j > 0L; j-- )
371 {
372 // out[j] = (out[j] + out[j-1]/2.0;
373 outVal = out.Get();
374 --out;
375 outVal2 = out.Get();
376 outVal = ( outVal + outVal2 ) / 2.0;
377 ++out;
378 out.Set( static_cast< OutputImagePixelType >( outVal ) );
379 --out;
380 progress.CompletedPixel();
381 }
382 // out[0] /= 2.0;
383 out.Set( static_cast< OutputImagePixelType >( out.Get() / 2.0 ) );
384
385 // TODO: Temporary fix for itkImageLinearIteratorWithIndex::NextLine() not
386 // setting m_Position correctly
387 out.GoToEndOfLine();
388 }
389 } // namespace itk
390
391 #endif
392