- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

interpolate.h
説明を見る。
1 //
2 // Copyright (c) 2003-2011, MIST Project, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28
33 #ifndef __INCLUDE_MIST_INTERPOLATE__
34 #define __INCLUDE_MIST_INTERPOLATE__
35
36
37 #ifndef __INCLUDE_MIST_H__
38 #include "mist.h"
39 #endif
40
41 #ifndef __INCLUDE_MIST_LIMITS__
42 #include "limits.h"
43 #endif
44
45
46 #ifndef __INCLUDE_MIST_THREAD__
47 #include "thread.h"
48 #endif
49
50
51 #include <cmath>
52
53 // mist名前空間の始まり
54 _MIST_BEGIN
55
56 // 補間関数で利用するユーティリティー関数群
57 namespace __interpolate_utility__
58 {
59 // 浮動小数の場合は四捨五入しない
60 template < bool b >
61 struct _round_
62 {
63 template < class T >
64 static void convert( double in, T &out )
65 {
66 out = static_cast< T >( in );
67 }
68
69 template < class T >
70 static void convert( const rgb< double > &in, rgb< T > &out )
71 {
72 out.r = static_cast< T >( in.r );
73 out.g = static_cast< T >( in.g );
74 out.b = static_cast< T >( in.b );
75 }
76
77 template < class T >
78 static void convert( const rgb< double > &in, rgba< T > &out )
79 {
80 out.r = static_cast< T >( in.r );
81 out.g = static_cast< T >( in.g );
82 out.b = static_cast< T >( in.b );
83 }
84
85 template < class T >
86 static void convert( const bgr< double > &in, bgr< T > &out )
87 {
88 out.r = static_cast< T >( in.r );
89 out.g = static_cast< T >( in.g );
90 out.b = static_cast< T >( in.b );
91 }
92
93 template < class T >
94 static void convert( const bgr< double > &in, bgra< T > &out )
95 {
96 out.r = static_cast< T >( in.r );
97 out.g = static_cast< T >( in.g );
98 out.b = static_cast< T >( in.b );
99 }
100 };
101
102 // 整数の場合は四捨五入する
103 template < >
104 struct _round_< true >
105 {
106 template < class T >
107 static void convert( double in, T &out )
108 {
109 out = static_cast< T >( in + 0.5 );
110 }
111
112 template < class T >
113 static void convert( const rgb< double > &in, rgb< T > &out )
114 {
115 out.r = static_cast< T >( in.r + 0.5 );
116 out.g = static_cast< T >( in.g + 0.5 );
117 out.b = static_cast< T >( in.b + 0.5 );
118 }
119
120 template < class T >
121 static void convert( const rgb< double > &in, rgba< T > &out )
122 {
123 out.r = static_cast< T >( in.r + 0.5 );
124 out.g = static_cast< T >( in.g + 0.5 );
125 out.b = static_cast< T >( in.b + 0.5 );
126 }
127
128 template < class T >
129 static void convert( const bgr< double > &in, bgr< T > &out )
130 {
131 out.r = static_cast< T >( in.r + 0.5 );
132 out.g = static_cast< T >( in.g + 0.5 );
133 out.b = static_cast< T >( in.b + 0.5 );
134 }
135
136 template < class T >
137 static void convert( const bgr< double > &in, bgra< T > &out )
138 {
139 out.r = static_cast< T >( in.r + 0.5 );
140 out.g = static_cast< T >( in.g + 0.5 );
141 out.b = static_cast< T >( in.b + 0.5 );
142 }
143 };
144
145 template < class T >
146 inline void round( const double in, T &out )
147 {
148 _round_< is_integer< T >::value >::convert( in, out );
149 }
150
151 template < class T >
152 inline void round( const rgb< double > &in, rgb< T > &out )
153 {
154 _round_< is_integer< T >::value >::convert( in, out );
155 }
156
157 template < class T >
158 inline void round( const rgb< double > &in, rgba< T > &out )
159 {
160 _round_< is_integer< T >::value >::convert( in, out );
161 }
162
163 template < class T >
164 inline void round( const bgr< double > &in, bgr< T > &out )
165 {
166 _round_< is_integer< T >::value >::convert( in, out );
167 }
168
169 template < class T >
170 inline void round( const bgr< double > &in, bgra< T > &out )
171 {
172 _round_< is_integer< T >::value >::convert( in, out );
173 }
174 }
175
176 // 最近傍型補間
177 namespace __nearest__
178 {
179 template < class Array1, class Array2 >
180 void interpolate( const Array1 &in, Array2 &out,
181 typename Array1::size_type thread_idx, typename Array1::size_type thread_numx,
182 typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
183 typename Array1::size_type thread_idz, typename Array1::size_type thread_numz )
184 {
185 typedef typename Array1::size_type size_type;
186 typedef typename Array1::value_type value_type;
187 typedef typename Array2::value_type out_value_type;
188
189 size_type i, j, k;
190 size_type iw = in.width( );
191 size_type ih = in.height( );
192 size_type id = in.depth( );
193 size_type ow = out.width( );
194 size_type oh = out.height( );
195 size_type od = out.depth( );
196
197 double sx = static_cast< double >( iw ) / static_cast< double >( ow );
198 double sy = static_cast< double >( ih ) / static_cast< double >( oh );
199 double sz = static_cast< double >( id ) / static_cast< double >( od );
200 size_type x, y, z;
201
202 for( k = thread_idz ; k < od ; k += thread_numz )
203 {
204 z = static_cast< size_type >( sz * k + 0.5 );
205 z = z < id ? z : id - 1;
206 for( j = thread_idy ; j < oh ; j += thread_numy )
207 {
208 y = static_cast< size_type >( sy * j + 0.5 );
209 y = y < ih ? y : ih - 1;
210 for( i = thread_idx ; i < ow ; i += thread_numx )
211 {
212 x = static_cast< size_type >( sx * i + 0.5 );
213 x = x < iw ? x : iw - 1;
214
215 __interpolate_utility__::round( in( x, y, z ), out( i, j, k ) );
216 }
217 }
218 }
219 }
220 }
221
222
223 // 平均値型補間
224 namespace __mean__
225 {
226 template < bool b >
227 struct _mean_
228 {
229 template < class T, class Allocator >
230 static double mean___( const array< T, Allocator > &in,
231 typename array< T, Allocator >::size_type i1,
232 typename array< T, Allocator >::size_type i2,
233 typename array< T, Allocator >::size_type j1,
234 typename array< T, Allocator >::size_type j2,
235 typename array< T, Allocator >::size_type k1,
236 typename array< T, Allocator >::size_type k2,
237 double xs,
238 double xe,
239 double ys,
240 double ye,
241 double zs,
242 double ze )
243 {
244 typedef typename array< T, Allocator >::value_type value_type;
245 typedef typename array< T, Allocator >::size_type size_type;
246 double pix = 0.0;
247
248 {
249 double a = ( i1 + 1 - xs ) * 0.5;
250 pix += ( in[ i1 ] + in[ i1 + 1 ] ) * a;
251 }
252
253 for( size_type i = i1 + 1 ; i < i2 - 1 ; i++ )
254 {
255 pix += ( in[ i ] + in[ i + 1 ] ) * 0.5;
256 }
257
258 {
259 double a = ( xe - i2 + 1 ) * 0.5;
260 pix += ( in[ i2 - 1 ] + in[ i2 ] ) * a;
261 }
262
263 double min = type_limits< value_type >::minimum( );
264 double max = type_limits< value_type >::maximum( );
265 pix /= ( xe - xs ) * ( ye - ys ) * ( ze - zs );
266 pix = pix > min ? pix : min;
267 pix = pix < max ? pix : max;
268 return( pix );
269 }
270
271 template < class T, class Allocator >
272 static double mean___( const array2< T, Allocator > &in,
273 typename array2< T, Allocator >::size_type i1,
274 typename array2< T, Allocator >::size_type i2,
275 typename array2< T, Allocator >::size_type j1,
276 typename array2< T, Allocator >::size_type j2,
277 typename array2< T, Allocator >::size_type /* k1 */,
278 typename array2< T, Allocator >::size_type /* k2 */,
279 double xs,
280 double xe,
281 double ys,
282 double ye,
283 double /* zs */,
284 double /* ze */ )
285 {
286 typedef typename array2< T, Allocator >::value_type value_type;
287 typedef typename array2< T, Allocator >::size_type size_type;
288 typedef typename array2< T, Allocator >::const_pointer const_pointer;
289 double pix = 0.0, sum = 0.0;
290 double yy;
291
292 if( i2 - i1 < 3 )
293 {
294 double xx[ 2 ] = { ( i1 + 1 - xs ) * 0.25, ( xe - i2 + 1 ) * 0.25 };
295 for( size_type j = j1 ; j < j2 ; j++ )
296 {
297 if( j == j1 )
298 {
299 yy = j1 + 1 - ys;
300 }
301 else if( j == j2 - 1 )
302 {
303 yy = ye - j;
304 }
305 else
306 {
307 yy = 1.0;
308 }
309
310 double aa[ 2 ] = { xx[ 0 ] * yy, xx[ 1 ] * yy };
311
312 const_pointer p0 = &in( 0, j );
313 const_pointer p1 = &in( 0, j + 1 );
314
315 pix += ( p0[ i1 ] + p1[ i1 ] + p0[ i1 + 1 ] + p1[ i1 + 1 ] ) * aa[ 0 ];
316 pix += ( p0[ i2 - 1 ] + p1[ i2 - 1 ] + p0[ i2 ] + p1[ i2 ] ) * aa[ 1 ];
317 sum += ( aa[ 0 ] + aa[ 1 ] ) * 4.0;
318 }
319 }
320 else
321 {
322 double xx[ 3 ] = { ( i1 + 1 - xs ) * 0.25, 0.25, ( xe - i2 + 1 ) * 0.25 };
323 for( size_type j = j1 ; j < j2 ; j++ )
324 {
325 if( j == j1 )
326 {
327 yy = j1 + 1 - ys;
328 }
329 else if( j == j2 - 1 )
330 {
331 yy = ye - j;
332 }
333 else
334 {
335 yy = 1.0;
336 }
337
338 double aa[ 3 ] = { xx[ 0 ] * yy, xx[ 1 ] * yy, xx[ 2 ] * yy };
339
340 const_pointer p0 = &in( 0, j );
341 const_pointer p1 = &in( 0, j + 1 );
342 double ppp;
343 size_type i = i1;
344
345 {
346 pix += ( p0[ i ] + p1[ i ] + p0[ i + 1 ] + p1[ i + 1 ] ) * aa[ 0 ];
347 ppp = p0[ i + 1 ] + p1[ i + 1 ];
348 }
349
350 for( i = i1 + 2 ; i < i2 - 1 ; i++ )
351 {
352 ppp += ( p0[ i ] + p1[ i ] ) * 2.0;
353 }
354
355 {
356 ppp += p0[ i2 - 1 ] + p1[ i2 - 1 ];
357 pix += ppp * aa[ 1 ];
358 pix += ( p0[ i2 - 1 ] + p1[ i2 - 1 ] + p0[ i2 ] + p1[ i2 ] ) * aa[ 2 ];
359 }
360
361 sum += ( aa[ 0 ] + aa[ 2 ] ) * 4.0 + aa[ 1 ] * ( i2 - i1 - 2 ) * 4.0;
362 }
363 }
364
365 double min = type_limits< value_type >::minimum( );
366 double max = type_limits< value_type >::maximum( );
367 pix /= sum;
368 pix = pix > min ? pix : min;
369 pix = pix < max ? pix : max;
370 return( pix );
371 }
372
373 template < class T, class Allocator >
374 static double mean___( const array3< T, Allocator > &in,
375 typename array3< T, Allocator >::size_type i1,
376 typename array3< T, Allocator >::size_type i2,
377 typename array3< T, Allocator >::size_type j1,
378 typename array3< T, Allocator >::size_type j2,
379 typename array3< T, Allocator >::size_type k1,
380 typename array3< T, Allocator >::size_type k2,
381 double xs,
382 double xe,
383 double ys,
384 double ye,
385 double zs,
386 double ze )
387 {
388 typedef typename array3< T, Allocator >::value_type value_type;
389 typedef typename array3< T, Allocator >::size_type size_type;
390 typedef typename array3< T, Allocator >::const_pointer const_pointer;
391 double pix = 0.0;
392 double xx[ 3 ] = { ( i1 + 1 - xs ) * 0.125, i2 - i1 < 3 ? 0.0 : 1.0 * 0.125, ( xe - i2 + 1 ) * 0.125 };
393 double yy, zz;
394
395 for( size_type k = k1 ; k < k2 ; k++ )
396 {
397 if( k == k1 )
398 {
399 zz = k1 + 1 - zs;
400 }
401 else if( k == k2 - 1 )
402 {
403 zz = ze - k;
404 }
405 else
406 {
407 zz = 1.0;
408 }
409
410 for( size_type j = j1 ; j < j2 ; j++ )
411 {
412 if( j == j1 )
413 {
414 yy = j1 + 1 - ys;
415 }
416 else if( j == j2 - 1 )
417 {
418 yy = ye - j;
419 }
420 else
421 {
422 yy = 1.0;
423 }
424
425 double aa[ 3 ] = { xx[ 0 ] * yy * zz, xx[ 1 ] * yy * zz, xx[ 2 ] * yy * zz };
426
427 const_pointer p0 = &in( 0, j , k );
428 const_pointer p1 = &in( 0, j + 1, k );
429 const_pointer p2 = &in( 0, j , k + 1 );
430 const_pointer p3 = &in( 0, j + 1, k + 1 );
431 double ppp = 0.0;
432
433 {
434 const value_type &c1 = p0[ i1 ];
435 const value_type &c2 = p1[ i1 ];
436 const value_type &c3 = p2[ i1 ];
437 const value_type &c4 = p3[ i1 ];
438 const value_type &c5 = p0[ i1 + 1 ];
439 const value_type &c6 = p1[ i1 + 1 ];
440 const value_type &c7 = p2[ i1 + 1 ];
441 const value_type &c8 = p3[ i1 + 1 ];
442
443 pix += ( c1 + c2 + c3 + c4 + c5 + c6 + c7 + c8 ) * aa[ 0 ];
444 ppp += c5 + c6 + c7 + c8;
445 }
446
447 for( size_type i = i1 + 2 ; i < i2 - 1 ; i++ )
448 {
449 const value_type &c1 = p0[ i ];
450 const value_type &c2 = p1[ i ];
451 const value_type &c3 = p2[ i ];
452 const value_type &c4 = p3[ i ];
453
454 ppp += ( c1 + c2 + c3 + c4 ) * 2.0;
455 }
456
457 {
458 const value_type &c1 = p0[ i2 - 1 ];
459 const value_type &c2 = p1[ i2 - 1 ];
460 const value_type &c3 = p2[ i2 - 1 ];
461 const value_type &c4 = p3[ i2 - 1 ];
462 const value_type &c5 = p0[ i2 ];
463 const value_type &c6 = p1[ i2 ];
464 const value_type &c7 = p2[ i2 ];
465 const value_type &c8 = p3[ i2 ];
466
467 ppp += c1 + c2 + c3 + c4;
468 pix += ppp * aa[ 1 ];
469 pix += ( c1 + c2 + c3 + c4 + c5 + c6 + c7 + c8 ) * aa[ 2 ];
470 }
471 }
472 }
473
474 double min = type_limits< value_type >::minimum( );
475 double max = type_limits< value_type >::maximum( );
476 pix /= ( xe - xs ) * ( ye - ys ) * ( ze - zs );
477 pix = pix > min ? pix : min;
478 pix = pix < max ? pix : max;
479 return( pix );
480 }
481 };
482
483 template < >
484 struct _mean_< true >
485 {
486 template < class T, class Allocator >
487 static const typename T::template rebind< double >::other mean___( const array< T, Allocator > &in,
488 typename array< T, Allocator >::size_type i1,
489 typename array< T, Allocator >::size_type i2,
490 typename array< T, Allocator >::size_type /* j1 */,
491 typename array< T, Allocator >::size_type /* j2 */,
492 typename array< T, Allocator >::size_type /* k1 */,
493 typename array< T, Allocator >::size_type /* k2 */,
494 double xs,
495 double xe,
496 double ys,
497 double ye,
498 double zs,
499 double ze )
500 {
501 typedef typename array< T, Allocator >::value_type color;
502 typedef typename array< T, Allocator >::size_type size_type;
503 typedef typename color::value_type value_type;
504 typedef typename T::template rebind< double >::other ovalue_type;
505 double r = 0.0;
506 double g = 0.0;
507 double b = 0.0;
508
509 {
510 double rr = 0.0;
511 double gg = 0.0;
512 double bb = 0.0;
513
514 const color &c1 = in[ i1 ];
515 const color &c2 = in[ i1 + 1 ];
516
517 rr += c1.r; gg += c1.g; bb += c1.b;
518 rr += c2.r; gg += c2.g; bb += c2.b;
519
520 double a = ( i1 + 1 - xs ) * 0.5;
521 r += rr * a;
522 g += gg * a;
523 b += bb * a;
524 }
525
526 for( size_type i = i1 + 1 ; i < i2 - 1 ; i++ )
527 {
528 double rr = 0.0;
529 double gg = 0.0;
530 double bb = 0.0;
531
532 const color &c1 = in[ i ];
533 const color &c2 = in[ i + 1 ];
534
535 rr += c1.r; gg += c1.g; bb += c1.b;
536 rr += c2.r; gg += c2.g; bb += c2.b;
537
538 double a = 0.5;
539 r += rr * a;
540 g += gg * a;
541 b += bb * a;
542 }
543
544 {
545 double rr = 0.0;
546 double gg = 0.0;
547 double bb = 0.0;
548
549 const color &c1 = in[ i2 - 1 ];
550 const color &c2 = in[ i2 ];
551
552 rr += c1.r; gg += c1.g; bb += c1.b;
553 rr += c2.r; gg += c2.g; bb += c2.b;
554
555 double a = ( xe - i2 + 1 ) * 0.5;
556 r += rr * a;
557 g += gg * a;
558 b += bb * a;
559 }
560
561 double numPixels = ( xe - xs ) * ( ye - ys ) * ( ze - zs );
562 double min = type_limits< value_type >::minimum( );
563 double max = type_limits< value_type >::maximum( );
564 r /= numPixels;
565 r = r > min ? r : min;
566 r = r < max ? r : max;
567 g /= numPixels;
568 g = g > min ? g : min;
569 g = g < max ? g : max;
570 b /= numPixels;
571 b = b > min ? b : min;
572 b = b < max ? b : max;
573 return( ovalue_type( r, g, b ) );
574 }
575
576 template < class T, class Allocator >
577 static const typename T::template rebind< double >::other mean___( const array2< T, Allocator > &in,
578 typename array2< T, Allocator >::size_type i1,
579 typename array2< T, Allocator >::size_type i2,
580 typename array2< T, Allocator >::size_type j1,
581 typename array2< T, Allocator >::size_type j2,
582 typename array2< T, Allocator >::size_type /* k1 */,
583 typename array2< T, Allocator >::size_type /* k2 */,
584 double xs,
585 double xe,
586 double ys,
587 double ye,
588 double /* zs*/,
589 double /* ze */ )
590 {
591 typedef typename array2< T, Allocator >::value_type color;
592 typedef typename array2< T, Allocator >::size_type size_type;
593 typedef typename array2< T, Allocator >::const_pointer const_pointer;
594 typedef typename color::value_type value_type;
595 typedef typename T::template rebind< double >::other ovalue_type;
596 double r = 0.0;
597 double g = 0.0;
598 double b = 0.0;
599 double yy, sum = 0.0;
600
601 if( i2 - i1 < 3 )
602 {
603 double xx[ 2 ] = { ( i1 + 1 - xs ) * 0.25, ( xe - i2 + 1 ) * 0.25 };
604 for( size_type j = j1 ; j < j2 ; j++ )
605 {
606 if( j == j1 )
607 {
608 yy = j1 + 1 - ys;
609 }
610 else if( j == j2 - 1 )
611 {
612 yy = ye - j;
613 }
614 else
615 {
616 yy = 1.0;
617 }
618
619 double aa[ 2 ] = { xx[ 0 ] * yy, xx[ 1 ] * yy };
620
621 const_pointer p0 = &in( 0, j );
622 const_pointer p1 = &in( 0, j + 1 );
623
624 {
625 const color &c1 = p0[ i1 ];
626 const color &c2 = p1[ i1 ];
627 const color &c3 = p0[ i1 + 1 ];
628 const color &c4 = p1[ i2 - 1 ];
629 const color &c5 = p0[ i2 ];
630 const color &c6 = p1[ i2 ];
631
632 r += ( c1.r + c2.r + c3.r + c4.r ) * aa[ 0 ];
633 g += ( c1.g + c2.g + c3.g + c4.g ) * aa[ 0 ];
634 b += ( c1.b + c2.b + c3.b + c4.b ) * aa[ 0 ];
635
636 r += ( c3.r + c4.r + c5.r + c6.r ) * aa[ 1 ];
637 g += ( c3.g + c4.g + c5.g + c6.g ) * aa[ 1 ];
638 b += ( c3.b + c4.b + c5.b + c6.b ) * aa[ 1 ];
639 }
640
641 sum += ( aa[ 0 ] + aa[ 1 ] ) * 4.0;
642 }
643 }
644 else
645 {
646 double xx[ 3 ] = { ( i1 + 1 - xs ) * 0.25, 0.25, ( xe - i2 + 1 ) * 0.25 };
647 for( size_type j = j1 ; j < j2 ; j++ )
648 {
649 if( j == j1 )
650 {
651 yy = j1 + 1 - ys;
652 }
653 else if( j == j2 - 1 )
654 {
655 yy = ye - j;
656 }
657 else
658 {
659 yy = 1.0;
660 }
661
662 double aa[ 3 ] = { xx[ 0 ] * yy, xx[ 1 ] * yy, xx[ 2 ] * yy };
663
664 const_pointer p0 = &in( 0, j );
665 const_pointer p1 = &in( 0, j + 1 );
666 double rr, gg, bb;
667
668 {
669 const color &c1 = p0[ i1 ];
670 const color &c2 = p1[ i1 ];
671 const color &c3 = p0[ i1 + 1 ];
672 const color &c4 = p1[ i1 + 1 ];
673
674 r += ( c1.r + c2.r + c3.r + c4.r ) * aa[ 0 ];
675 g += ( c1.g + c2.g + c3.g + c4.g ) * aa[ 0 ];
676 b += ( c1.b + c2.b + c3.b + c4.b ) * aa[ 0 ];
677
678 rr = c3.r + c4.r;
679 gg = c3.g + c4.g;
680 bb = c3.b + c4.b;
681 }
682
683 for( size_type i = i1 + 2 ; i < i2 - 1 ; i++ )
684 {
685 const color &c1 = p0[ i ];
686 const color &c2 = p1[ i ];
687
688 rr += ( c1.r + c2.r ) * 2.0;
689 gg += ( c1.g + c2.g ) * 2.0;
690 bb += ( c1.b + c2.b ) * 2.0;
691 }
692
693 {
694 const color &c1 = p0[ i2 - 1 ];
695 const color &c2 = p1[ i2 - 1 ];
696 const color &c3 = p0[ i2 ];
697 const color &c4 = p1[ i2 ];
698
699 rr += c1.r + c2.r;
700 gg += c1.g + c2.g;
701 bb += c1.b + c2.b;
702
703 r += rr * aa[ 1 ];
704 g += gg * aa[ 1 ];
705 b += bb * aa[ 1 ];
706
707 r += ( c1.r + c2.r + c3.r + c4.r ) * aa[ 2 ];
708 g += ( c1.g + c2.g + c3.g + c4.g ) * aa[ 2 ];
709 b += ( c1.b + c2.b + c3.b + c4.b ) * aa[ 2 ];
710 }
711
712 sum += ( aa[ 0 ] + aa[ 2 ] ) * 4.0 + aa[ 1 ] * ( i2 - i1 - 2 ) * 4.0;
713 }
714 }
715
716 double min = type_limits< value_type >::minimum( );
717 double max = type_limits< value_type >::maximum( );
718 r /= sum;
719 r = r > min ? r : min;
720 r = r < max ? r : max;
721 g /= sum;
722 g = g > min ? g : min;
723 g = g < max ? g : max;
724 b /= sum;
725 b = b > min ? b : min;
726 b = b < max ? b : max;
727 return( ovalue_type( r, g, b ) );
728 }
729
730 template < class T, class Allocator >
731 static const typename T::template rebind< double >::other mean___( const array3< T, Allocator > &in,
732 typename array3< T, Allocator >::size_type i1,
733 typename array3< T, Allocator >::size_type i2,
734 typename array3< T, Allocator >::size_type j1,
735 typename array3< T, Allocator >::size_type j2,
736 typename array3< T, Allocator >::size_type k1,
737 typename array3< T, Allocator >::size_type k2,
738 double xs,
739 double xe,
740 double ys,
741 double ye,
742 double zs,
743 double ze )
744 {
745 typedef typename array3< T, Allocator >::value_type color;
746 typedef typename array3< T, Allocator >::size_type size_type;
747 typedef typename array3< T, Allocator >::const_pointer const_pointer;
748 typedef typename color::value_type value_type;
749 typedef typename T::template rebind< double >::other ovalue_type;
750 double r = 0.0;
751 double g = 0.0;
752 double b = 0.0;
753 double yy, zz;
754 double xx[ 3 ] = { ( i1 + 1 - xs ) * 0.125, i2 - i1 < 3 ? 0.0 : 1.0 * 0.125, ( xe - i2 + 1 ) * 0.125 };
755
756 for( size_type k = k1 ; k < k2 ; k++ )
757 {
758 if( k == k1 )
759 {
760 zz = k1 + 1 - zs;
761 }
762 else if( k == k2 - 1 )
763 {
764 zz = ze - k;
765 }
766 else
767 {
768 zz = 1.0;
769 }
770
771 for( size_type j = j1 ; j < j2 ; j++ )
772 {
773 if( j == j1 )
774 {
775 yy = j1 + 1 - ys;
776 }
777 else if( j == j2 - 1 )
778 {
779 yy = ye - j;
780 }
781 else
782 {
783 yy = 1.0;
784 }
785
786 double aa[ 3 ] = { xx[ 0 ] * yy * zz, xx[ 1 ] * yy * zz, xx[ 2 ] * yy * zz };
787
788 const_pointer p0 = &in( 0, j , k );
789 const_pointer p1 = &in( 0, j + 1, k );
790 const_pointer p2 = &in( 0, j , k + 1 );
791 const_pointer p3 = &in( 0, j + 1, k + 1 );
792 double rr = 0.0;
793 double gg = 0.0;
794 double bb = 0.0;
795
796 {
797 const color &c1 = p0[ i1 ];
798 const color &c2 = p1[ i1 ];
799 const color &c3 = p2[ i1 ];
800 const color &c4 = p3[ i1 ];
801 const color &c5 = p0[ i1 + 1 ];
802 const color &c6 = p1[ i1 + 1 ];
803 const color &c7 = p2[ i1 + 1 ];
804 const color &c8 = p3[ i1 + 1 ];
805
806 r += ( c1.r + c2.r + c3.r + c4.r + c5.r + c6.r + c7.r + c8.r ) * aa[ 0 ];
807 g += ( c1.g + c2.g + c3.g + c4.g + c5.g + c6.g + c7.g + c8.g ) * aa[ 0 ];
808 b += ( c1.b + c2.b + c3.b + c4.b + c5.b + c6.b + c7.b + c8.b ) * aa[ 0 ];
809
810 rr += c5.r + c6.r + c7.r + c8.r;
811 gg += c5.g + c6.g + c7.g + c8.g;
812 bb += c5.b + c6.b + c7.b + c8.b;
813 }
814
815 for( size_type i = i1 + 2 ; i < i2 - 1 ; i++ )
816 {
817 const color &c1 = p0[ i ];
818 const color &c2 = p1[ i ];
819 const color &c3 = p2[ i ];
820 const color &c4 = p3[ i ];
821
822 rr += ( c1.r + c2.r + c3.r + c4.r ) * 2.0;
823 gg += ( c1.g + c2.g + c3.g + c4.g ) * 2.0;
824 bb += ( c1.b + c2.b + c3.b + c4.b ) * 2.0;
825 }
826
827 {
828 const color &c1 = p0[ i2 - 1 ];
829 const color &c2 = p1[ i2 - 1 ];
830 const color &c3 = p2[ i2 - 1 ];
831 const color &c4 = p3[ i2 - 1 ];
832 const color &c5 = p0[ i2 ];
833 const color &c6 = p1[ i2 ];
834 const color &c7 = p2[ i2 ];
835 const color &c8 = p3[ i2 ];
836
837 rr += c1.r + c2.r + c3.r + c4.r;
838 gg += c1.g + c2.g + c3.g + c4.g;
839 bb += c1.b + c2.b + c3.b + c4.b;
840
841 r += rr * aa[ 1 ];
842 g += gg * aa[ 1 ];
843 b += bb * aa[ 1 ];
844
845 r += ( c1.r + c2.r + c3.r + c4.r + c5.r + c6.r + c7.r + c8.r ) * aa[ 2 ];
846 g += ( c1.g + c2.g + c3.g + c4.g + c5.g + c6.g + c7.g + c8.g ) * aa[ 2 ];
847 b += ( c1.b + c2.b + c3.b + c4.b + c5.b + c6.b + c7.b + c8.b ) * aa[ 2 ];
848 }
849 }
850 }
851
852 double numPixels = ( xe - xs ) * ( ye - ys ) * ( ze - zs );
853 double min = type_limits< value_type >::minimum( );
854 double max = type_limits< value_type >::maximum( );
855 r /= numPixels;
856 r = r > min ? r : min;
857 r = r < max ? r : max;
858 g /= numPixels;
859 g = g > min ? g : min;
860 g = g < max ? g : max;
861 b /= numPixels;
862 b = b > min ? b : min;
863 b = b < max ? b : max;
864 return( ovalue_type( r, g, b ) );
865 }
866 };
867
868 template < class T1, class Allocator1, class T2, class Allocator2 >
869 void interpolate( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
870 typename array< T1, Allocator1 >::size_type thread_idx, typename array< T1, Allocator1 >::size_type thread_numx,
871 typename array< T1, Allocator1 >::size_type thread_idy, typename array< T1, Allocator1 >::size_type thread_numy,
872 typename array< T1, Allocator1 >::size_type thread_idz, typename array< T1, Allocator1 >::size_type thread_numz )
873 {
874 typedef typename array< T1, Allocator1 >::size_type size_type;
875 typedef typename array< T1, Allocator1 >::value_type value_type;
876 typedef typename array< T2, Allocator2 >::value_type out_value_type;
877
878 size_type i, i1, i2;
879 size_type iw = in.width( );
880 size_type ow = out.width( );
881
882 double sx = static_cast< double >( iw - 1 ) / static_cast< double >( ow );
883
884 for( i = thread_idx ; i < ow ; i += thread_numx )
885 {
886 double xs = sx * i;
887 double xe = xs + sx;
888
889 i1 = static_cast< size_type >( xs );
890 i2 = static_cast< size_type >( xe );
891 i2 = xe > i2 ? i2 + 1 : i2;
892 i2 = i2 < iw - 1 ? i2 : iw - 1;
893
894 __interpolate_utility__::round( _mean_< is_color< value_type >::value >::mean___( in, i1, i2, 0, 1, 0, 1, xs, xe, 0, 1, 0, 1 ), out[ i ] );
895 }
896 }
897
898 template < class T1, class Allocator1, class T2, class Allocator2 >
899 void interpolate( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
900 typename array2< T1, Allocator1 >::size_type thread_idx, typename array2< T1, Allocator1 >::size_type thread_numx,
901 typename array2< T1, Allocator1 >::size_type thread_idy, typename array2< T1, Allocator1 >::size_type thread_numy,
902 typename array2< T1, Allocator1 >::size_type /* thread_idz */, typename array2< T1, Allocator1 >::size_type /* thread_numz */ )
903 {
904 typedef typename array2< T1, Allocator1 >::size_type size_type;
905 typedef typename array2< T1, Allocator1 >::value_type value_type;
906 typedef typename array2< T2, Allocator2 >::value_type out_value_type;
907
908 size_type i, j, i1, i2, j1, j2;
909 size_type iw = in.width( );
910 size_type ih = in.height( );
911 size_type ow = out.width( );
912 size_type oh = out.height( );
913
914 double sx = static_cast< double >( iw ) / static_cast< double >( ow );
915 double sy = static_cast< double >( ih ) / static_cast< double >( oh );
916
917 for( j = thread_idy ; j < oh ; j += thread_numy )
918 {
919 double ys = sy * j;
920 double ye = ys + sy;
921
922 j1 = static_cast< size_type >( ys );
923 j2 = static_cast< size_type >( ye );
924 j2 = ye > j2 ? j2 + 1 : j2;
925 j2 = j2 < ih - 1 ? j2 : ih - 1;
926
927 for( i = thread_idx ; i < ow ; i += thread_numx )
928 {
929 double xs = sx * i;
930 double xe = xs + sx;
931
932 i1 = static_cast< size_type >( xs );
933 i2 = static_cast< size_type >( xe );
934 i2 = xe > i2 ? i2 + 1 : i2;
935 i2 = i2 < iw - 1 ? i2 : iw - 1;
936
937 __interpolate_utility__::round( _mean_< is_color< value_type >::value >::mean___( in, i1, i2, j1, j2, 0, 1, xs, xe, ys, ye, 0, 1 ), out( i, j ) );
938 }
939 }
940 }
941
942 template < class T1, class Allocator1, class T2, class Allocator2 >
943 void interpolate( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
944 typename array3< T1, Allocator1 >::size_type thread_idx, typename array3< T1, Allocator1 >::size_type thread_numx,
945 typename array3< T1, Allocator1 >::size_type thread_idy, typename array3< T1, Allocator1 >::size_type thread_numy,
946 typename array3< T1, Allocator1 >::size_type thread_idz, typename array3< T1, Allocator1 >::size_type thread_numz )
947 {
948 typedef typename array3< T1, Allocator1 >::size_type size_type;
949 typedef typename array3< T1, Allocator1 >::value_type value_type;
950 typedef typename array3< T2, Allocator2 >::value_type out_value_type;
951
952 size_type i, j, k, i1, i2, j1, j2, k1, k2;
953 size_type iw = in.width( );
954 size_type ih = in.height( );
955 size_type id = in.depth( );
956 size_type ow = out.width( );
957 size_type oh = out.height( );
958 size_type od = out.depth( );
959
960 double sx = static_cast< double >( iw - 1 ) / static_cast< double >( ow );
961 double sy = static_cast< double >( ih - 1 ) / static_cast< double >( oh );
962 double sz = static_cast< double >( id - 1 ) / static_cast< double >( od );
963
964 for( k = thread_idz ; k < od ; k += thread_numz )
965 {
966 double zs = sz * k;
967 double ze = zs + sz;
968
969 k1 = static_cast< size_type >( zs );
970 k2 = static_cast< size_type >( ze );
971 k2 = ze > k2 ? k2 + 1 : k2;
972 k2 = k2 < id - 1 ? k2 : id - 1;
973
974 for( j = thread_idy ; j < oh ; j += thread_numy )
975 {
976 double ys = sy * j;
977 double ye = ys + sy;
978
979 j1 = static_cast< size_type >( ys );
980 j2 = static_cast< size_type >( ye );
981 j2 = ye > j2 ? j2 + 1 : j2;
982 j2 = j2 < ih - 1 ? j2 : ih - 1;
983
984 for( i = thread_idx ; i < ow ; i += thread_numx )
985 {
986 double xs = sx * i;
987 double xe = xs + sx;
988
989 i1 = static_cast< size_type >( xs );
990 i2 = static_cast< size_type >( xe );
991 i2 = xe > i2 ? i2 + 1 : i2;
992 i2 = i2 < iw - 1 ? i2 : iw - 1;
993
994 __interpolate_utility__::round( _mean_< is_color< value_type >::value >::mean___( in, i1, i2, j1, j2, k1, k2, xs, xe, ys, ye, zs, ze ), out( i, j, k ) );
995 }
996 }
997 }
998 }
999 }
1000
1001 // 線形補間
1002 namespace __linear__
1003 {
1004 template < bool b >
1005 struct _linear_
1006 {
1007 template < class T, class Allocator >
1008 static double interpolate( const array< T, Allocator > &in,
1009 typename array< T, Allocator >::size_type i1,
1010 typename array< T, Allocator >::size_type i2,
1011 typename array< T, Allocator >::size_type j1,
1012 typename array< T, Allocator >::size_type j2,
1013 typename array< T, Allocator >::size_type k1,
1014 typename array< T, Allocator >::size_type k2,
1015 double x, double y, double z )
1016 {
1017 typedef typename array< T, Allocator >::value_type value_type;
1018 double min = type_limits< value_type >::minimum( );
1019 double max = type_limits< value_type >::maximum( );
1020 double pix = in[ i1 ] * ( 1.0 - x ) + in[ i2 ] * x;
1021 pix = pix > min ? pix : min;
1022 pix = pix < max ? pix : max;
1023 return( pix );
1024 }
1025
1026 template < class T, class Allocator >
1027 static double interpolate( const array2< T, Allocator > &in,
1028 typename array2< T, Allocator >::size_type i1,
1029 typename array2< T, Allocator >::size_type i2,
1030 typename array2< T, Allocator >::size_type j1,
1031 typename array2< T, Allocator >::size_type j2,
1032 typename array2< T, Allocator >::size_type /* k1 */,
1033 typename array2< T, Allocator >::size_type /* k2 */,
1034 double x, double y, double /* z */ )
1035 {
1036 typedef typename array2< T, Allocator >::value_type value_type;
1037 double min = type_limits< value_type >::minimum( );
1038 double max = type_limits< value_type >::maximum( );
1039 double pix = ( in( i1, j1 ) * ( 1.0 - x ) + in( i2, j1 ) * x ) * ( 1.0 - y ) + ( in( i1, j2 ) * ( 1.0 - x ) + in( i2, j2 ) * x ) * y;
1040 pix = pix > min ? pix : min;
1041 pix = pix < max ? pix : max;
1042 return( pix );
1043 }
1044
1045 template < class T, class Allocator >
1046 static double interpolate( const array3< T, Allocator > &in,
1047 typename array3< T, Allocator >::size_type i1,
1048 typename array3< T, Allocator >::size_type i2,
1049 typename array3< T, Allocator >::size_type j1,
1050 typename array3< T, Allocator >::size_type j2,
1051 typename array3< T, Allocator >::size_type k1,
1052 typename array3< T, Allocator >::size_type k2,
1053 double x, double y, double z )
1054 {
1055 typedef typename array3< T, Allocator >::value_type value_type;
1056 double min = type_limits< value_type >::minimum( );
1057 double max = type_limits< value_type >::maximum( );
1058 double pix = ( ( in( i1, j1, k1 ) * ( 1.0 - x ) + in( i2, j1, k1 ) * x ) * ( 1.0 - y ) + ( in( i1, j2, k1 ) * ( 1.0 - x ) + in( i2, j2, k1 ) * x ) * y ) * ( 1.0 - z )
1059 + ( ( in( i1, j1, k2 ) * ( 1.0 - x ) + in( i2, j1, k2 ) * x ) * ( 1.0 - y ) + ( in( i1, j2, k2 ) * ( 1.0 - x ) + in( i2, j2, k2 ) * x ) * y ) * z;
1060 pix = pix > min ? pix : min;
1061 pix = pix < max ? pix : max;
1062 return( pix );
1063 }
1064 };
1065
1066 template < >
1067 struct _linear_< true >
1068 {
1069 template < class T, class Allocator >
1070 static const typename T::template rebind< double >::other interpolate( const array< T, Allocator > &in,
1071 typename array< T, Allocator >::size_type i1,
1072 typename array< T, Allocator >::size_type i2,
1073 typename array< T, Allocator >::size_type /* j1 */,
1074 typename array< T, Allocator >::size_type /* j2 */,
1075 typename array< T, Allocator >::size_type /* k1 */,
1076 typename array< T, Allocator >::size_type /* k2 */,
1077 double x, double /* y */, double /* z */ )
1078 {
1079 typedef typename array< T, Allocator >::value_type color;
1080 typedef typename color::value_type value_type;
1081 typedef typename T::template rebind< double >::other ovalue_type;
1082 double min = type_limits< value_type >::minimum( );
1083 double max = type_limits< value_type >::maximum( );
1084 double r = in[ i1 ].r * ( 1.0 - x ) + in[ i2 ].r * x;
1085 double g = in[ i1 ].g * ( 1.0 - x ) + in[ i2 ].g * x;
1086 double b = in[ i1 ].b * ( 1.0 - x ) + in[ i2 ].b * x;
1087 r = r > min ? r : min;
1088 r = r < max ? r : max;
1089 g = g > min ? g : min;
1090 g = g < max ? g : max;
1091 b = b > min ? b : min;
1092 b = b < max ? b : max;
1093 return( ovalue_type( r, g, b ) );
1094 }
1095
1096 template < class T, class Allocator >
1097 static const typename T::template rebind< double >::other interpolate( const array2< T, Allocator > &in,
1098 typename array2< T, Allocator >::size_type i1,
1099 typename array2< T, Allocator >::size_type i2,
1100 typename array2< T, Allocator >::size_type j1,
1101 typename array2< T, Allocator >::size_type j2,
1102 typename array2< T, Allocator >::size_type /* k1 */,
1103 typename array2< T, Allocator >::size_type /* k2 */,
1104 double x, double y, double /* z */ )
1105 {
1106 typedef typename array2< T, Allocator >::value_type color;
1107 typedef typename color::value_type value_type;
1108 typedef typename T::template rebind< double >::other ovalue_type;
1109 double min = type_limits< value_type >::minimum( );
1110 double max = type_limits< value_type >::maximum( );
1111 double r = ( in( i1, j1 ).r * ( 1.0 - x ) + in( i2, j1 ).r * x ) * ( 1.0 - y ) + ( in( i1, j2 ).r * ( 1.0 - x ) + in( i2, j2 ).r * x ) * y;
1112 double g = ( in( i1, j1 ).g * ( 1.0 - x ) + in( i2, j1 ).g * x ) * ( 1.0 - y ) + ( in( i1, j2 ).g * ( 1.0 - x ) + in( i2, j2 ).g * x ) * y;
1113 double b = ( in( i1, j1 ).b * ( 1.0 - x ) + in( i2, j1 ).b * x ) * ( 1.0 - y ) + ( in( i1, j2 ).b * ( 1.0 - x ) + in( i2, j2 ).b * x ) * y;
1114 r = r > min ? r : min;
1115 r = r < max ? r : max;
1116 g = g > min ? g : min;
1117 g = g < max ? g : max;
1118 b = b > min ? b : min;
1119 b = b < max ? b : max;
1120 return( ovalue_type( r, g, b ) );
1121 }
1122
1123 template < class T, class Allocator >
1124 static const typename T::template rebind< double >::other interpolate( const array3< T, Allocator > &in,
1125 typename array3< T, Allocator >::size_type i1,
1126 typename array3< T, Allocator >::size_type i2,
1127 typename array3< T, Allocator >::size_type j1,
1128 typename array3< T, Allocator >::size_type j2,
1129 typename array3< T, Allocator >::size_type k1,
1130 typename array3< T, Allocator >::size_type k2,
1131 double x, double y, double z )
1132 {
1133 typedef typename array3< T, Allocator >::value_type color;
1134 typedef typename color::value_type value_type;
1135 typedef typename T::template rebind< double >::other ovalue_type;
1136 double min = type_limits< value_type >::minimum( );
1137 double max = type_limits< value_type >::maximum( );
1138 double r = ( ( in( i1, j1, k1 ).r * ( 1.0 - x ) + in( i2, j1, k1 ).r * x ) * ( 1.0 - y ) + ( in( i1, j2, k1 ).r * ( 1.0 - x ) + in( i2, j2, k1 ).r * x ) * y ) * ( 1.0 - z )
1139 + ( ( in( i1, j1, k2 ).r * ( 1.0 - x ) + in( i2, j1, k2 ).r * x ) * ( 1.0 - y ) + ( in( i1, j2, k2 ).r * ( 1.0 - x ) + in( i2, j2, k2 ).r * x ) * y ) * z;
1140 double g = ( ( in( i1, j1, k1 ).g * ( 1.0 - x ) + in( i2, j1, k1 ).g * x ) * ( 1.0 - y ) + ( in( i1, j2, k1 ).g * ( 1.0 - x ) + in( i2, j2, k1 ).g * x ) * y ) * ( 1.0 - z )
1141 + ( ( in( i1, j1, k2 ).g * ( 1.0 - x ) + in( i2, j1, k2 ).g * x ) * ( 1.0 - y ) + ( in( i1, j2, k2 ).g * ( 1.0 - x ) + in( i2, j2, k2 ).g * x ) * y ) * z;
1142 double b = ( ( in( i1, j1, k1 ).b * ( 1.0 - x ) + in( i2, j1, k1 ).b * x ) * ( 1.0 - y ) + ( in( i1, j2, k1 ).b * ( 1.0 - x ) + in( i2, j2, k1 ).b * x ) * y ) * ( 1.0 - z )
1143 + ( ( in( i1, j1, k2 ).b * ( 1.0 - x ) + in( i2, j1, k2 ).b * x ) * ( 1.0 - y ) + ( in( i1, j2, k2 ).b * ( 1.0 - x ) + in( i2, j2, k2 ).b * x ) * y ) * z;
1144 r = r > min ? r : min;
1145 r = r < max ? r : max;
1146 g = g > min ? g : min;
1147 g = g < max ? g : max;
1148 b = b > min ? b : min;
1149 b = b < max ? b : max;
1150 return( ovalue_type( r, g, b ) );
1151 }
1152 };
1153
1154 template < class Array1, class Array2 >
1155 void interpolate( const Array1 &in, Array2 &out,
1156 typename Array1::size_type thread_idx, typename Array1::size_type thread_numx,
1157 typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
1158 typename Array1::size_type thread_idz, typename Array1::size_type thread_numz )
1159 {
1160 typedef typename Array1::size_type size_type;
1161 typedef typename Array1::value_type value_type;
1162 typedef typename Array2::value_type out_value_type;
1163
1164 size_type i, j, k, i1, i2, j1, j2, k1, k2;
1165 size_type iw = in.width( );
1166 size_type ih = in.height( );
1167 size_type id = in.depth( );
1168 size_type ow = out.width( );
1169 size_type oh = out.height( );
1170 size_type od = out.depth( );
1171
1172 double sx = static_cast< double >( iw ) / static_cast< double >( ow );
1173 double sy = static_cast< double >( ih ) / static_cast< double >( oh );
1174 double sz = static_cast< double >( id ) / static_cast< double >( od );
1175 double x, y, z;
1176
1177 for( k = thread_idz ; k < od ; k += thread_numz )
1178 {
1179 z = sz * k;
1180 k1 = static_cast< size_type >( z );
1181 z -= k1;
1182 k2 = k1 < id - 1 ? k1 + 1 : k1;
1183 for( j = thread_idy ; j < oh ; j += thread_numy )
1184 {
1185 y = sy * j;
1186 j1 = static_cast< size_type >( y );
1187 y -= j1;
1188 j2 = j1 < ih - 1 ? j1 + 1 : j1;
1189 for( i = thread_idx ; i < ow ; i += thread_numx )
1190 {
1191 x = sx * i;
1192 i1 = static_cast< size_type >( x );
1193 x -= i1;
1194 i2 = i1 < iw - 1 ? i1 + 1 : i1;
1195
1196 __interpolate_utility__::round( _linear_< is_color< value_type >::value >::interpolate( in, i1, i2, j1, j2, k1, k2, x, y, z ), out( i, j, k ) );
1197 }
1198 }
1199 }
1200 }
1201 }
1202
1203
1204
1205 // 3次のsinc関数による内挿補間
1206 namespace __cubic__
1207 {
1208 // 0 <= t < 1 の場合
1209 inline double sinc1( double t ){ return( 1.0 + ( t - 2.0 ) * t * t ); }
1210
1211 // 1 <= t < 2 の場合
1212 inline double sinc2( double t ){ return( 4.0 + ( -8.0 + ( 5.0 - t ) * t ) * t ); }
1213
1214
1215 template < bool b >
1216 struct _cubic_
1217 {
1218 template < class T, class Allocator >
1219 static double interpolate( const array< T, Allocator > &in,
1220 typename array< T, Allocator >::size_type i[4],
1221 typename array< T, Allocator >::size_type /* j */[4],
1222 typename array< T, Allocator >::size_type /* k */[4],
1223 double x, double /* y */, double /* z */ )
1224 {
1225 typedef typename array< T, Allocator >::value_type value_type;
1226
1227 double u0 = sinc2( 1 + x );
1228 double u1 = sinc1( x );
1229 double u2 = sinc1( 1 - x );
1230 double u3 = sinc2( 2 - x );
1231 double pix = in[ i[ 0 ] ] * u0 + in[ i[ 1 ] ] * u1 + in[ i[ 2 ] ] * u2 + in[ i[ 3 ] ] * u3;
1232
1233 double min = type_limits< value_type >::minimum( );
1234 double max = type_limits< value_type >::maximum( );
1235 pix = pix > min ? pix : min;
1236 pix = pix < max ? pix : max;
1237
1238 return( pix );
1239 }
1240
1241 template < class T, class Allocator >
1242 static double interpolate( const array2< T, Allocator > &in,
1243 typename array2< T, Allocator >::size_type i[4],
1244 typename array2< T, Allocator >::size_type j[4],
1245 typename array2< T, Allocator >::size_type /* k */[4],
1246 double x, double y, double /* z */ )
1247 {
1248 typedef typename array2< T, Allocator >::value_type value_type;
1249
1250 double u0 = sinc2( 1 + x );
1251 double u1 = sinc1( x );
1252 double u2 = sinc1( 1 - x );
1253 double u3 = sinc2( 2 - x );
1254 double v0 = sinc2( 1 + y );
1255 double v1 = sinc1( y );
1256 double v2 = sinc1( 1 - y );
1257 double v3 = sinc2( 2 - y );
1258 double pix = ( in( i[ 0 ], j[ 0 ] ) * u0 + in( i[ 1 ], j[ 0 ] ) * u1 + in( i[ 2 ], j[ 0 ] ) * u2 + in( i[ 3 ], j[ 0 ] ) * u3 ) * v0
1259 + ( in( i[ 0 ], j[ 1 ] ) * u0 + in( i[ 1 ], j[ 1 ] ) * u1 + in( i[ 2 ], j[ 1 ] ) * u2 + in( i[ 3 ], j[ 1 ] ) * u3 ) * v1
1260 + ( in( i[ 0 ], j[ 2 ] ) * u0 + in( i[ 1 ], j[ 2 ] ) * u1 + in( i[ 2 ], j[ 2 ] ) * u2 + in( i[ 3 ], j[ 2 ] ) * u3 ) * v2
1261 + ( in( i[ 0 ], j[ 3 ] ) * u0 + in( i[ 1 ], j[ 3 ] ) * u1 + in( i[ 2 ], j[ 3 ] ) * u2 + in( i[ 3 ], j[ 3 ] ) * u3 ) * v3;
1262
1263 double min = type_limits< value_type >::minimum( );
1264 double max = type_limits< value_type >::maximum( );
1265 pix = pix > min ? pix : min;
1266 pix = pix < max ? pix : max;
1267
1268 return( pix );
1269 }
1270
1271 template < class T, class Allocator >
1272 static double interpolate( const array3< T, Allocator > &in,
1273 typename array3< T, Allocator >::size_type i[4],
1274 typename array3< T, Allocator >::size_type j[4],
1275 typename array3< T, Allocator >::size_type k[4],
1276 double x, double y, double z )
1277 {
1278 typedef typename array3< T, Allocator >::value_type value_type;
1279
1280 double u0 = sinc2( 1 + x );
1281 double u1 = sinc1( x );
1282 double u2 = sinc1( 1 - x );
1283 double u3 = sinc2( 2 - x );
1284 double v0 = sinc2( 1 + y );
1285 double v1 = sinc1( y );
1286 double v2 = sinc1( 1 - y );
1287 double v3 = sinc2( 2 - y );
1288 double w0 = sinc2( 1 + z );
1289 double w1 = sinc1( z );
1290 double w2 = sinc1( 1 - z );
1291 double w3 = sinc2( 2 - z );
1292
1293 double p0 = ( ( in( i[ 0 ], j[ 0 ], k[ 0 ] ) * u0 + in( i[ 1 ], j[ 0 ], k[ 0 ] ) * u1 + in( i[ 2 ], j[ 0 ], k[ 0 ] ) * u2 + in( i[ 3 ], j[ 0 ], k[ 0 ] ) * u3 ) * v0
1294 + ( in( i[ 0 ], j[ 1 ], k[ 0 ] ) * u0 + in( i[ 1 ], j[ 1 ], k[ 0 ] ) * u1 + in( i[ 2 ], j[ 1 ], k[ 0 ] ) * u2 + in( i[ 3 ], j[ 1 ], k[ 0 ] ) * u3 ) * v1
1295 + ( in( i[ 0 ], j[ 2 ], k[ 0 ] ) * u0 + in( i[ 1 ], j[ 2 ], k[ 0 ] ) * u1 + in( i[ 2 ], j[ 2 ], k[ 0 ] ) * u2 + in( i[ 3 ], j[ 2 ], k[ 0 ] ) * u3 ) * v2
1296 + ( in( i[ 0 ], j[ 3 ], k[ 0 ] ) * u0 + in( i[ 1 ], j[ 3 ], k[ 0 ] ) * u1 + in( i[ 2 ], j[ 3 ], k[ 0 ] ) * u2 + in( i[ 3 ], j[ 3 ], k[ 0 ] ) * u3 ) * v3 );
1297 double p1 = ( ( in( i[ 0 ], j[ 0 ], k[ 1 ] ) * u0 + in( i[ 1 ], j[ 0 ], k[ 1 ] ) * u1 + in( i[ 2 ], j[ 0 ], k[ 1 ] ) * u2 + in( i[ 3 ], j[ 0 ], k[ 1 ] ) * u3 ) * v0
1298 + ( in( i[ 0 ], j[ 1 ], k[ 1 ] ) * u0 + in( i[ 1 ], j[ 1 ], k[ 1 ] ) * u1 + in( i[ 2 ], j[ 1 ], k[ 1 ] ) * u2 + in( i[ 3 ], j[ 1 ], k[ 1 ] ) * u3 ) * v1
1299 + ( in( i[ 0 ], j[ 2 ], k[ 1 ] ) * u0 + in( i[ 1 ], j[ 2 ], k[ 1 ] ) * u1 + in( i[ 2 ], j[ 2 ], k[ 1 ] ) * u2 + in( i[ 3 ], j[ 2 ], k[ 1 ] ) * u3 ) * v2
1300 + ( in( i[ 0 ], j[ 3 ], k[ 1 ] ) * u0 + in( i[ 1 ], j[ 3 ], k[ 1 ] ) * u1 + in( i[ 2 ], j[ 3 ], k[ 1 ] ) * u2 + in( i[ 3 ], j[ 3 ], k[ 1 ] ) * u3 ) * v3 );
1301 double p2 = ( ( in( i[ 0 ], j[ 0 ], k[ 2 ] ) * u0 + in( i[ 1 ], j[ 0 ], k[ 2 ] ) * u1 + in( i[ 2 ], j[ 0 ], k[ 2 ] ) * u2 + in( i[ 3 ], j[ 0 ], k[ 2 ] ) * u3 ) * v0
1302 + ( in( i[ 0 ], j[ 1 ], k[ 2 ] ) * u0 + in( i[ 1 ], j[ 1 ], k[ 2 ] ) * u1 + in( i[ 2 ], j[ 1 ], k[ 2 ] ) * u2 + in( i[ 3 ], j[ 1 ], k[ 2 ] ) * u3 ) * v1
1303 + ( in( i[ 0 ], j[ 2 ], k[ 2 ] ) * u0 + in( i[ 1 ], j[ 2 ], k[ 2 ] ) * u1 + in( i[ 2 ], j[ 2 ], k[ 2 ] ) * u2 + in( i[ 3 ], j[ 2 ], k[ 2 ] ) * u3 ) * v2
1304 + ( in( i[ 0 ], j[ 3 ], k[ 2 ] ) * u0 + in( i[ 1 ], j[ 3 ], k[ 2 ] ) * u1 + in( i[ 2 ], j[ 3 ], k[ 2 ] ) * u2 + in( i[ 3 ], j[ 3 ], k[ 2 ] ) * u3 ) * v3 );
1305 double p3 = ( ( in( i[ 0 ], j[ 0 ], k[ 3 ] ) * u0 + in( i[ 1 ], j[ 0 ], k[ 3 ] ) * u1 + in( i[ 2 ], j[ 0 ], k[ 3 ] ) * u2 + in( i[ 3 ], j[ 0 ], k[ 3 ] ) * u3 ) * v0
1306 + ( in( i[ 0 ], j[ 1 ], k[ 3 ] ) * u0 + in( i[ 1 ], j[ 1 ], k[ 3 ] ) * u1 + in( i[ 2 ], j[ 1 ], k[ 3 ] ) * u2 + in( i[ 3 ], j[ 1 ], k[ 3 ] ) * u3 ) * v1
1307 + ( in( i[ 0 ], j[ 2 ], k[ 3 ] ) * u0 + in( i[ 1 ], j[ 2 ], k[ 3 ] ) * u1 + in( i[ 2 ], j[ 2 ], k[ 3 ] ) * u2 + in( i[ 3 ], j[ 2 ], k[ 3 ] ) * u3 ) * v2
1308 + ( in( i[ 0 ], j[ 3 ], k[ 3 ] ) * u0 + in( i[ 1 ], j[ 3 ], k[ 3 ] ) * u1 + in( i[ 2 ], j[ 3 ], k[ 3 ] ) * u2 + in( i[ 3 ], j[ 3 ], k[ 3 ] ) * u3 ) * v3 );
1309 double pix = p0 * w0 + p1 * w1 + p2 * w2 + p3 * w3;
1310
1311 double min = type_limits< value_type >::minimum( );
1312 double max = type_limits< value_type >::maximum( );
1313 pix = pix > min ? pix : min;
1314 pix = pix < max ? pix : max;
1315
1316 return( pix );
1317 }
1318 };
1319
1320 template < >
1321 struct _cubic_< true >
1322 {
1323 template < class T, class Allocator >
1324 static const typename T::template rebind< double >::other interpolate( const array< T, Allocator > &in,
1325 typename array< T, Allocator >::size_type i[4],
1326 typename array< T, Allocator >::size_type /* j */[4],
1327 typename array< T, Allocator >::size_type /* k */[4],
1328 double x, double /* y */, double /* z */ )
1329 {
1330 typedef typename array< T, Allocator >::value_type color;
1331 typedef typename color::value_type value_type;
1332 typedef typename T::template rebind< double >::other ovalue_type;
1333
1334 double u0 = sinc2( 1 + x );
1335 double u1 = sinc1( x );
1336 double u2 = sinc1( 1 - x );
1337 double u3 = sinc2( 2 - x );
1338 double r = in[ i[ 0 ] ].r * u0 + in[ i[ 1 ] ].r * u1 + in[ i[ 2 ] ].r * u2 + in[ i[ 3 ] ].r * u3;
1339 double g = in[ i[ 0 ] ].g * u0 + in[ i[ 1 ] ].g * u1 + in[ i[ 2 ] ].g * u2 + in[ i[ 3 ] ].g * u3;
1340 double b = in[ i[ 0 ] ].b * u0 + in[ i[ 1 ] ].b * u1 + in[ i[ 2 ] ].b * u2 + in[ i[ 3 ] ].b * u3;
1341
1342 double min = type_limits< value_type >::minimum( );
1343 double max = type_limits< value_type >::maximum( );
1344
1345 r = r > min ? r : min;
1346 r = r < max ? r : max;
1347 g = g > min ? g : min;
1348 g = g < max ? g : max;
1349 b = b > min ? b : min;
1350 b = b < max ? b : max;
1351
1352 return( ovalue_type( r, g, b ) );
1353 }
1354
1355 template < class T, class Allocator >
1356 static const typename T::template rebind< double >::other interpolate( const array2< T, Allocator > &in,
1357 typename array2< T, Allocator >::size_type i[4],
1358 typename array2< T, Allocator >::size_type j[4],
1359 typename array2< T, Allocator >::size_type /* k */[4],
1360 double x, double y, double /* z */ )
1361 {
1362 typedef typename array2< T, Allocator >::value_type color;
1363 typedef typename color::value_type value_type;
1364 typedef typename T::template rebind< double >::other ovalue_type;
1365
1366 double u0 = sinc2( 1 + x );
1367 double u1 = sinc1( x );
1368 double u2 = sinc1( 1 - x );
1369 double u3 = sinc2( 2 - x );
1370 double v0 = sinc2( 1 + y );
1371 double v1 = sinc1( y );
1372 double v2 = sinc1( 1 - y );
1373 double v3 = sinc2( 2 - y );
1374
1375 double r = ( in( i[ 0 ], j[ 0 ] ).r * u0 + in( i[ 1 ], j[ 0 ] ).r * u1 + in( i[ 2 ], j[ 0 ] ).r * u2 + in( i[ 3 ], j[ 0 ] ).r * u3 ) * v0
1376 + ( in( i[ 0 ], j[ 1 ] ).r * u0 + in( i[ 1 ], j[ 1 ] ).r * u1 + in( i[ 2 ], j[ 1 ] ).r * u2 + in( i[ 3 ], j[ 1 ] ).r * u3 ) * v1
1377 + ( in( i[ 0 ], j[ 2 ] ).r * u0 + in( i[ 1 ], j[ 2 ] ).r * u1 + in( i[ 2 ], j[ 2 ] ).r * u2 + in( i[ 3 ], j[ 2 ] ).r * u3 ) * v2
1378 + ( in( i[ 0 ], j[ 3 ] ).r * u0 + in( i[ 1 ], j[ 3 ] ).r * u1 + in( i[ 2 ], j[ 3 ] ).r * u2 + in( i[ 3 ], j[ 3 ] ).r * u3 ) * v3;
1379 double g = ( in( i[ 0 ], j[ 0 ] ).g * u0 + in( i[ 1 ], j[ 0 ] ).g * u1 + in( i[ 2 ], j[ 0 ] ).g * u2 + in( i[ 3 ], j[ 0 ] ).g * u3 ) * v0
1380 + ( in( i[ 0 ], j[ 1 ] ).g * u0 + in( i[ 1 ], j[ 1 ] ).g * u1 + in( i[ 2 ], j[ 1 ] ).g * u2 + in( i[ 3 ], j[ 1 ] ).g * u3 ) * v1
1381 + ( in( i[ 0 ], j[ 2 ] ).g * u0 + in( i[ 1 ], j[ 2 ] ).g * u1 + in( i[ 2 ], j[ 2 ] ).g * u2 + in( i[ 3 ], j[ 2 ] ).g * u3 ) * v2
1382 + ( in( i[ 0 ], j[ 3 ] ).g * u0 + in( i[ 1 ], j[ 3 ] ).g * u1 + in( i[ 2 ], j[ 3 ] ).g * u2 + in( i[ 3 ], j[ 3 ] ).g * u3 ) * v3;
1383 double b = ( in( i[ 0 ], j[ 0 ] ).b * u0 + in( i[ 1 ], j[ 0 ] ).b * u1 + in( i[ 2 ], j[ 0 ] ).b * u2 + in( i[ 3 ], j[ 0 ] ).b * u3 ) * v0
1384 + ( in( i[ 0 ], j[ 1 ] ).b * u0 + in( i[ 1 ], j[ 1 ] ).b * u1 + in( i[ 2 ], j[ 1 ] ).b * u2 + in( i[ 3 ], j[ 1 ] ).b * u3 ) * v1
1385 + ( in( i[ 0 ], j[ 2 ] ).b * u0 + in( i[ 1 ], j[ 2 ] ).b * u1 + in( i[ 2 ], j[ 2 ] ).b * u2 + in( i[ 3 ], j[ 2 ] ).b * u3 ) * v2
1386 + ( in( i[ 0 ], j[ 3 ] ).b * u0 + in( i[ 1 ], j[ 3 ] ).b * u1 + in( i[ 2 ], j[ 3 ] ).b * u2 + in( i[ 3 ], j[ 3 ] ).b * u3 ) * v3;
1387
1388 double min = type_limits< value_type >::minimum( );
1389 double max = type_limits< value_type >::maximum( );
1390
1391 r = r > min ? r : min;
1392 r = r < max ? r : max;
1393 g = g > min ? g : min;
1394 g = g < max ? g : max;
1395 b = b > min ? b : min;
1396 b = b < max ? b : max;
1397
1398 return( ovalue_type( r, g, b ) );
1399 }
1400
1401 template < class T, class Allocator >
1402 static const typename T::template rebind< double >::other interpolate( const array3< T, Allocator > &in,
1403 typename array3< T, Allocator >::size_type i[4],
1404 typename array3< T, Allocator >::size_type j[4],
1405 typename array3< T, Allocator >::size_type k[4],
1406 double x, double y, double z )
1407 {
1408 typedef typename array3< T, Allocator >::value_type color;
1409 typedef typename color::value_type value_type;
1410 typedef typename T::template rebind< double >::other ovalue_type;
1411
1412 double u0 = sinc2( 1 + x );
1413 double u1 = sinc1( x );
1414 double u2 = sinc1( 1 - x );
1415 double u3 = sinc2( 2 - x );
1416 double v0 = sinc2( 1 + y );
1417 double v1 = sinc1( y );
1418 double v2 = sinc1( 1 - y );
1419 double v3 = sinc2( 2 - y );
1420 double w0 = sinc2( 1 + z );
1421 double w1 = sinc1( z );
1422 double w2 = sinc1( 1 - z );
1423 double w3 = sinc2( 2 - z );
1424
1425 double r0 = ( ( in( i[ 0 ], j[ 0 ], k[ 0 ] ).r * u0 + in( i[ 1 ], j[ 0 ], k[ 0 ] ).r * u1 + in( i[ 2 ], j[ 0 ], k[ 0 ] ).r * u2 + in( i[ 3 ], j[ 0 ], k[ 0 ] ).r * u3 ) * v0
1426 + ( in( i[ 0 ], j[ 1 ], k[ 0 ] ).r * u0 + in( i[ 1 ], j[ 1 ], k[ 0 ] ).r * u1 + in( i[ 2 ], j[ 1 ], k[ 0 ] ).r * u2 + in( i[ 3 ], j[ 1 ], k[ 0 ] ).r * u3 ) * v1
1427 + ( in( i[ 0 ], j[ 2 ], k[ 0 ] ).r * u0 + in( i[ 1 ], j[ 2 ], k[ 0 ] ).r * u1 + in( i[ 2 ], j[ 2 ], k[ 0 ] ).r * u2 + in( i[ 3 ], j[ 2 ], k[ 0 ] ).r * u3 ) * v2
1428 + ( in( i[ 0 ], j[ 3 ], k[ 0 ] ).r * u0 + in( i[ 1 ], j[ 3 ], k[ 0 ] ).r * u1 + in( i[ 2 ], j[ 3 ], k[ 0 ] ).r * u2 + in( i[ 3 ], j[ 3 ], k[ 0 ] ).r * u3 ) * v3 );
1429 double r1 = ( ( in( i[ 0 ], j[ 0 ], k[ 1 ] ).r * u0 + in( i[ 1 ], j[ 0 ], k[ 1 ] ).r * u1 + in( i[ 2 ], j[ 0 ], k[ 1 ] ).r * u2 + in( i[ 3 ], j[ 0 ], k[ 1 ] ).r * u3 ) * v0
1430 + ( in( i[ 0 ], j[ 1 ], k[ 1 ] ).r * u0 + in( i[ 1 ], j[ 1 ], k[ 1 ] ).r * u1 + in( i[ 2 ], j[ 1 ], k[ 1 ] ).r * u2 + in( i[ 3 ], j[ 1 ], k[ 1 ] ).r * u3 ) * v1
1431 + ( in( i[ 0 ], j[ 2 ], k[ 1 ] ).r * u0 + in( i[ 1 ], j[ 2 ], k[ 1 ] ).r * u1 + in( i[ 2 ], j[ 2 ], k[ 1 ] ).r * u2 + in( i[ 3 ], j[ 2 ], k[ 1 ] ).r * u3 ) * v2
1432 + ( in( i[ 0 ], j[ 3 ], k[ 1 ] ).r * u0 + in( i[ 1 ], j[ 3 ], k[ 1 ] ).r * u1 + in( i[ 2 ], j[ 3 ], k[ 1 ] ).r * u2 + in( i[ 3 ], j[ 3 ], k[ 1 ] ).r * u3 ) * v3 );
1433 double r2 = ( ( in( i[ 0 ], j[ 0 ], k[ 2 ] ).r * u0 + in( i[ 1 ], j[ 0 ], k[ 2 ] ).r * u1 + in( i[ 2 ], j[ 0 ], k[ 2 ] ).r * u2 + in( i[ 3 ], j[ 0 ], k[ 2 ] ).r * u3 ) * v0
1434 + ( in( i[ 0 ], j[ 1 ], k[ 2 ] ).r * u0 + in( i[ 1 ], j[ 1 ], k[ 2 ] ).r * u1 + in( i[ 2 ], j[ 1 ], k[ 2 ] ).r * u2 + in( i[ 3 ], j[ 1 ], k[ 2 ] ).r * u3 ) * v1
1435 + ( in( i[ 0 ], j[ 2 ], k[ 2 ] ).r * u0 + in( i[ 1 ], j[ 2 ], k[ 2 ] ).r * u1 + in( i[ 2 ], j[ 2 ], k[ 2 ] ).r * u2 + in( i[ 3 ], j[ 2 ], k[ 2 ] ).r * u3 ) * v2
1436 + ( in( i[ 0 ], j[ 3 ], k[ 2 ] ).r * u0 + in( i[ 1 ], j[ 3 ], k[ 2 ] ).r * u1 + in( i[ 2 ], j[ 3 ], k[ 2 ] ).r * u2 + in( i[ 3 ], j[ 3 ], k[ 2 ] ).r * u3 ) * v3 );
1437 double r3 = ( ( in( i[ 0 ], j[ 0 ], k[ 3 ] ).r * u0 + in( i[ 1 ], j[ 0 ], k[ 3 ] ).r * u1 + in( i[ 2 ], j[ 0 ], k[ 3 ] ).r * u2 + in( i[ 3 ], j[ 0 ], k[ 3 ] ).r * u3 ) * v0
1438 + ( in( i[ 0 ], j[ 1 ], k[ 3 ] ).r * u0 + in( i[ 1 ], j[ 1 ], k[ 3 ] ).r * u1 + in( i[ 2 ], j[ 1 ], k[ 3 ] ).r * u2 + in( i[ 3 ], j[ 1 ], k[ 3 ] ).r * u3 ) * v1
1439 + ( in( i[ 0 ], j[ 2 ], k[ 3 ] ).r * u0 + in( i[ 1 ], j[ 2 ], k[ 3 ] ).r * u1 + in( i[ 2 ], j[ 2 ], k[ 3 ] ).r * u2 + in( i[ 3 ], j[ 2 ], k[ 3 ] ).r * u3 ) * v2
1440 + ( in( i[ 0 ], j[ 3 ], k[ 3 ] ).r * u0 + in( i[ 1 ], j[ 3 ], k[ 3 ] ).r * u1 + in( i[ 2 ], j[ 3 ], k[ 3 ] ).r * u2 + in( i[ 3 ], j[ 3 ], k[ 3 ] ).r * u3 ) * v3 );
1441 double g0 = ( ( in( i[ 0 ], j[ 0 ], k[ 0 ] ).g * u0 + in( i[ 1 ], j[ 0 ], k[ 0 ] ).g * u1 + in( i[ 2 ], j[ 0 ], k[ 0 ] ).g * u2 + in( i[ 3 ], j[ 0 ], k[ 0 ] ).g * u3 ) * v0
1442 + ( in( i[ 0 ], j[ 1 ], k[ 0 ] ).g * u0 + in( i[ 1 ], j[ 1 ], k[ 0 ] ).g * u1 + in( i[ 2 ], j[ 1 ], k[ 0 ] ).g * u2 + in( i[ 3 ], j[ 1 ], k[ 0 ] ).g * u3 ) * v1
1443 + ( in( i[ 0 ], j[ 2 ], k[ 0 ] ).g * u0 + in( i[ 1 ], j[ 2 ], k[ 0 ] ).g * u1 + in( i[ 2 ], j[ 2 ], k[ 0 ] ).g * u2 + in( i[ 3 ], j[ 2 ], k[ 0 ] ).g * u3 ) * v2
1444 + ( in( i[ 0 ], j[ 3 ], k[ 0 ] ).g * u0 + in( i[ 1 ], j[ 3 ], k[ 0 ] ).g * u1 + in( i[ 2 ], j[ 3 ], k[ 0 ] ).g * u2 + in( i[ 3 ], j[ 3 ], k[ 0 ] ).g * u3 ) * v3 );
1445 double g1 = ( ( in( i[ 0 ], j[ 0 ], k[ 1 ] ).g * u0 + in( i[ 1 ], j[ 0 ], k[ 1 ] ).g * u1 + in( i[ 2 ], j[ 0 ], k[ 1 ] ).g * u2 + in( i[ 3 ], j[ 0 ], k[ 1 ] ).g * u3 ) * v0
1446 + ( in( i[ 0 ], j[ 1 ], k[ 1 ] ).g * u0 + in( i[ 1 ], j[ 1 ], k[ 1 ] ).g * u1 + in( i[ 2 ], j[ 1 ], k[ 1 ] ).g * u2 + in( i[ 3 ], j[ 1 ], k[ 1 ] ).g * u3 ) * v1
1447 + ( in( i[ 0 ], j[ 2 ], k[ 1 ] ).g * u0 + in( i[ 1 ], j[ 2 ], k[ 1 ] ).g * u1 + in( i[ 2 ], j[ 2 ], k[ 1 ] ).g * u2 + in( i[ 3 ], j[ 2 ], k[ 1 ] ).g * u3 ) * v2
1448 + ( in( i[ 0 ], j[ 3 ], k[ 1 ] ).g * u0 + in( i[ 1 ], j[ 3 ], k[ 1 ] ).g * u1 + in( i[ 2 ], j[ 3 ], k[ 1 ] ).g * u2 + in( i[ 3 ], j[ 3 ], k[ 1 ] ).g * u3 ) * v3 );
1449 double g2 = ( ( in( i[ 0 ], j[ 0 ], k[ 2 ] ).g * u0 + in( i[ 1 ], j[ 0 ], k[ 2 ] ).g * u1 + in( i[ 2 ], j[ 0 ], k[ 2 ] ).g * u2 + in( i[ 3 ], j[ 0 ], k[ 2 ] ).g * u3 ) * v0
1450 + ( in( i[ 0 ], j[ 1 ], k[ 2 ] ).g * u0 + in( i[ 1 ], j[ 1 ], k[ 2 ] ).g * u1 + in( i[ 2 ], j[ 1 ], k[ 2 ] ).g * u2 + in( i[ 3 ], j[ 1 ], k[ 2 ] ).g * u3 ) * v1
1451 + ( in( i[ 0 ], j[ 2 ], k[ 2 ] ).g * u0 + in( i[ 1 ], j[ 2 ], k[ 2 ] ).g * u1 + in( i[ 2 ], j[ 2 ], k[ 2 ] ).g * u2 + in( i[ 3 ], j[ 2 ], k[ 2 ] ).g * u3 ) * v2
1452 + ( in( i[ 0 ], j[ 3 ], k[ 2 ] ).g * u0 + in( i[ 1 ], j[ 3 ], k[ 2 ] ).g * u1 + in( i[ 2 ], j[ 3 ], k[ 2 ] ).g * u2 + in( i[ 3 ], j[ 3 ], k[ 2 ] ).g * u3 ) * v3 );
1453 double g3 = ( ( in( i[ 0 ], j[ 0 ], k[ 3 ] ).g * u0 + in( i[ 1 ], j[ 0 ], k[ 3 ] ).g * u1 + in( i[ 2 ], j[ 0 ], k[ 3 ] ).g * u2 + in( i[ 3 ], j[ 0 ], k[ 3 ] ).g * u3 ) * v0
1454 + ( in( i[ 0 ], j[ 1 ], k[ 3 ] ).g * u0 + in( i[ 1 ], j[ 1 ], k[ 3 ] ).g * u1 + in( i[ 2 ], j[ 1 ], k[ 3 ] ).g * u2 + in( i[ 3 ], j[ 1 ], k[ 3 ] ).g * u3 ) * v1
1455 + ( in( i[ 0 ], j[ 2 ], k[ 3 ] ).g * u0 + in( i[ 1 ], j[ 2 ], k[ 3 ] ).g * u1 + in( i[ 2 ], j[ 2 ], k[ 3 ] ).g * u2 + in( i[ 3 ], j[ 2 ], k[ 3 ] ).g * u3 ) * v2
1456 + ( in( i[ 0 ], j[ 3 ], k[ 3 ] ).g * u0 + in( i[ 1 ], j[ 3 ], k[ 3 ] ).g * u1 + in( i[ 2 ], j[ 3 ], k[ 3 ] ).g * u2 + in( i[ 3 ], j[ 3 ], k[ 3 ] ).g * u3 ) * v3 );
1457 double b0 = ( ( in( i[ 0 ], j[ 0 ], k[ 0 ] ).b * u0 + in( i[ 1 ], j[ 0 ], k[ 0 ] ).b * u1 + in( i[ 2 ], j[ 0 ], k[ 0 ] ).b * u2 + in( i[ 3 ], j[ 0 ], k[ 0 ] ).b * u3 ) * v0
1458 + ( in( i[ 0 ], j[ 1 ], k[ 0 ] ).b * u0 + in( i[ 1 ], j[ 1 ], k[ 0 ] ).b * u1 + in( i[ 2 ], j[ 1 ], k[ 0 ] ).b * u2 + in( i[ 3 ], j[ 1 ], k[ 0 ] ).b * u3 ) * v1
1459 + ( in( i[ 0 ], j[ 2 ], k[ 0 ] ).b * u0 + in( i[ 1 ], j[ 2 ], k[ 0 ] ).b * u1 + in( i[ 2 ], j[ 2 ], k[ 0 ] ).b * u2 + in( i[ 3 ], j[ 2 ], k[ 0 ] ).b * u3 ) * v2
1460 + ( in( i[ 0 ], j[ 3 ], k[ 0 ] ).b * u0 + in( i[ 1 ], j[ 3 ], k[ 0 ] ).b * u1 + in( i[ 2 ], j[ 3 ], k[ 0 ] ).b * u2 + in( i[ 3 ], j[ 3 ], k[ 0 ] ).b * u3 ) * v3 );
1461 double b1 = ( ( in( i[ 0 ], j[ 0 ], k[ 1 ] ).b * u0 + in( i[ 1 ], j[ 0 ], k[ 1 ] ).b * u1 + in( i[ 2 ], j[ 0 ], k[ 1 ] ).b * u2 + in( i[ 3 ], j[ 0 ], k[ 1 ] ).b * u3 ) * v0
1462 + ( in( i[ 0 ], j[ 1 ], k[ 1 ] ).b * u0 + in( i[ 1 ], j[ 1 ], k[ 1 ] ).b * u1 + in( i[ 2 ], j[ 1 ], k[ 1 ] ).b * u2 + in( i[ 3 ], j[ 1 ], k[ 1 ] ).b * u3 ) * v1
1463 + ( in( i[ 0 ], j[ 2 ], k[ 1 ] ).b * u0 + in( i[ 1 ], j[ 2 ], k[ 1 ] ).b * u1 + in( i[ 2 ], j[ 2 ], k[ 1 ] ).b * u2 + in( i[ 3 ], j[ 2 ], k[ 1 ] ).b * u3 ) * v2
1464 + ( in( i[ 0 ], j[ 3 ], k[ 1 ] ).b * u0 + in( i[ 1 ], j[ 3 ], k[ 1 ] ).b * u1 + in( i[ 2 ], j[ 3 ], k[ 1 ] ).b * u2 + in( i[ 3 ], j[ 3 ], k[ 1 ] ).b * u3 ) * v3 );
1465 double b2 = ( ( in( i[ 0 ], j[ 0 ], k[ 2 ] ).b * u0 + in( i[ 1 ], j[ 0 ], k[ 2 ] ).b * u1 + in( i[ 2 ], j[ 0 ], k[ 2 ] ).b * u2 + in( i[ 3 ], j[ 0 ], k[ 2 ] ).b * u3 ) * v0
1466 + ( in( i[ 0 ], j[ 1 ], k[ 2 ] ).b * u0 + in( i[ 1 ], j[ 1 ], k[ 2 ] ).b * u1 + in( i[ 2 ], j[ 1 ], k[ 2 ] ).b * u2 + in( i[ 3 ], j[ 1 ], k[ 2 ] ).b * u3 ) * v1
1467 + ( in( i[ 0 ], j[ 2 ], k[ 2 ] ).b * u0 + in( i[ 1 ], j[ 2 ], k[ 2 ] ).b * u1 + in( i[ 2 ], j[ 2 ], k[ 2 ] ).b * u2 + in( i[ 3 ], j[ 2 ], k[ 2 ] ).b * u3 ) * v2
1468 + ( in( i[ 0 ], j[ 3 ], k[ 2 ] ).b * u0 + in( i[ 1 ], j[ 3 ], k[ 2 ] ).b * u1 + in( i[ 2 ], j[ 3 ], k[ 2 ] ).b * u2 + in( i[ 3 ], j[ 3 ], k[ 2 ] ).b * u3 ) * v3 );
1469 double b3 = ( ( in( i[ 0 ], j[ 0 ], k[ 3 ] ).b * u0 + in( i[ 1 ], j[ 0 ], k[ 3 ] ).b * u1 + in( i[ 2 ], j[ 0 ], k[ 3 ] ).b * u2 + in( i[ 3 ], j[ 0 ], k[ 3 ] ).b * u3 ) * v0
1470 + ( in( i[ 0 ], j[ 1 ], k[ 3 ] ).b * u0 + in( i[ 1 ], j[ 1 ], k[ 3 ] ).b * u1 + in( i[ 2 ], j[ 1 ], k[ 3 ] ).b * u2 + in( i[ 3 ], j[ 1 ], k[ 3 ] ).b * u3 ) * v1
1471 + ( in( i[ 0 ], j[ 2 ], k[ 3 ] ).b * u0 + in( i[ 1 ], j[ 2 ], k[ 3 ] ).b * u1 + in( i[ 2 ], j[ 2 ], k[ 3 ] ).b * u2 + in( i[ 3 ], j[ 2 ], k[ 3 ] ).b * u3 ) * v2
1472 + ( in( i[ 0 ], j[ 3 ], k[ 3 ] ).b * u0 + in( i[ 1 ], j[ 3 ], k[ 3 ] ).b * u1 + in( i[ 2 ], j[ 3 ], k[ 3 ] ).b * u2 + in( i[ 3 ], j[ 3 ], k[ 3 ] ).b * u3 ) * v3 );
1473
1474 double r = r0 * w0 + r1 * w1 + r2 * w2 + r3 * w3;
1475 double g = g0 * w0 + g1 * w1 + g2 * w2 + g3 * w3;
1476 double b = b0 * w0 + b1 * w1 + b2 * w2 + b3 * w3;
1477
1478 double min = type_limits< value_type >::minimum( );
1479 double max = type_limits< value_type >::maximum( );
1480
1481 r = r > min ? r : min;
1482 r = r < max ? r : max;
1483 g = g > min ? g : min;
1484 g = g < max ? g : max;
1485 b = b > min ? b : min;
1486 b = b < max ? b : max;
1487
1488 return( ovalue_type( r, g, b ) );
1489 }
1490 };
1491
1492 template < class Array1, class Array2 >
1493 void interpolate( const Array1 &in, Array2 &out,
1494 typename Array1::size_type thread_idx, typename Array1::size_type thread_numx,
1495 typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
1496 typename Array1::size_type thread_idz, typename Array1::size_type thread_numz )
1497 {
1498 typedef typename Array1::size_type size_type;
1499 typedef typename Array1::value_type value_type;
1500 typedef typename Array2::value_type out_value_type;
1501
1502 size_type i, j, k, ii[ 4 ], jj[ 4 ], kk[ 4 ];
1503 size_type iw = in.width( );
1504 size_type ih = in.height( );
1505 size_type id = in.depth( );
1506 size_type ow = out.width( );
1507 size_type oh = out.height( );
1508 size_type od = out.depth( );
1509
1510 double sx = static_cast< double >( iw ) / static_cast< double >( ow );
1511 double sy = static_cast< double >( ih ) / static_cast< double >( oh );
1512 double sz = static_cast< double >( id ) / static_cast< double >( od );
1513 double x, y, z;
1514
1515 for( k = thread_idz ; k < od ; k += thread_numz )
1516 {
1517 z = sz * k;
1518 kk[ 1 ] = static_cast< size_type >( z );
1519 kk[ 0 ] = kk[ 1 ] > 0 ? kk[ 1 ] - 1 : kk[ 1 ];
1520 kk[ 2 ] = kk[ 1 ] < id - 1 ? kk[ 1 ] + 1 : kk[ 1 ];
1521 kk[ 3 ] = kk[ 2 ] < id - 1 ? kk[ 2 ] + 1 : kk[ 2 ];
1522 z -= kk[ 1 ];
1523
1524 for( j = thread_idy ; j < oh ; j += thread_numy )
1525 {
1526 y = sy * j;
1527 jj[ 1 ] = static_cast< size_type >( y );
1528 jj[ 0 ] = jj[ 1 ] > 0 ? jj[ 1 ] - 1 : jj[ 1 ];
1529 jj[ 2 ] = jj[ 1 ] < ih - 1 ? jj[ 1 ] + 1 : jj[ 1 ];
1530 jj[ 3 ] = jj[ 2 ] < ih - 1 ? jj[ 2 ] + 1 : jj[ 2 ];
1531 y -= jj[ 1 ];
1532
1533 for( i = thread_idx ; i < ow ; i += thread_numx )
1534 {
1535 x = sx * i;
1536 ii[ 1 ] = static_cast< size_type >( x );
1537 ii[ 0 ] = ii[ 1 ] > 0 ? ii[ 1 ] - 1 : ii[ 1 ];
1538 ii[ 2 ] = ii[ 1 ] < iw - 1 ? ii[ 1 ] + 1 : ii[ 1 ];
1539 ii[ 3 ] = ii[ 2 ] < iw - 1 ? ii[ 2 ] + 1 : ii[ 2 ];
1540 x -= ii[ 1 ];
1541
1542 __interpolate_utility__::round( _cubic_< is_color< value_type >::value >::interpolate( in, ii, jj, kk, x, y, z ), out( i, j, k ) );
1543 }
1544 }
1545 }
1546 }
1547 }
1548
1549
1550
1551
1552
1553
1554 // 3次のBスプライン関数による内挿補間
1555 namespace __bspline__
1556 {
1557 // 0 <= t < 1 の場合
1558 inline double bspline1( double t ){ return( 2.0 / 3.0 + ( 0.5 * t - 1.0 ) * t * t ); }
1559
1560 // 1 <= t < 2 の場合
1561 inline double bspline2( double t ){ return( 4.0 / 3.0 + ( -2.0 + ( 1.0 - t / 6.0 ) * t ) * t ); }
1562
1563
1564 template < bool b >
1565 struct _bspline_
1566 {
1567 template < class T, class Allocator >
1568 static double interpolate( const array< T, Allocator > &in,
1569 typename array< T, Allocator >::size_type i[4],
1570 typename array< T, Allocator >::size_type /* j */[4],
1571 typename array< T, Allocator >::size_type /* k */[4],
1572 double x, double /* y */, double /* z */ )
1573 {
1574 typedef typename array< T, Allocator >::value_type value_type;
1575
1576 double u0 = bspline2( 1 + x );
1577 double u1 = bspline1( x );
1578 double u2 = bspline1( 1 - x );
1579 double u3 = bspline2( 2 - x );
1580 double pix = in[ i[ 0 ] ] * u0 + in[ i[ 1 ] ] * u1 + in[ i[ 2 ] ] * u2 + in[ i[ 3 ] ] * u3;
1581
1582 double min = type_limits< value_type >::minimum( );
1583 double max = type_limits< value_type >::maximum( );
1584 pix = pix > min ? pix : min;
1585 pix = pix < max ? pix : max;
1586
1587 return( pix );
1588 }
1589
1590 template < class T, class Allocator >
1591 static double interpolate( const array2< T, Allocator > &in,
1592 typename array2< T, Allocator >::size_type i[4],
1593 typename array2< T, Allocator >::size_type j[4],
1594 typename array2< T, Allocator >::size_type /* k */[4],
1595 double x, double y, double /* z */ )
1596 {
1597 typedef typename array2< T, Allocator >::value_type value_type;
1598
1599 double u0 = bspline2( 1 + x );
1600 double u1 = bspline1( x );
1601 double u2 = bspline1( 1 - x );
1602 double u3 = bspline2( 2 - x );
1603 double v0 = bspline2( 1 + y );
1604 double v1 = bspline1( y );
1605 double v2 = bspline1( 1 - y );
1606 double v3 = bspline2( 2 - y );
1607 double pix = ( in( i[ 0 ], j[ 0 ] ) * u0 + in( i[ 1 ], j[ 0 ] ) * u1 + in( i[ 2 ], j[ 0 ] ) * u2 + in( i[ 3 ], j[ 0 ] ) * u3 ) * v0
1608 + ( in( i[ 0 ], j[ 1 ] ) * u0 + in( i[ 1 ], j[ 1 ] ) * u1 + in( i[ 2 ], j[ 1 ] ) * u2 + in( i[ 3 ], j[ 1 ] ) * u3 ) * v1
1609 + ( in( i[ 0 ], j[ 2 ] ) * u0 + in( i[ 1 ], j[ 2 ] ) * u1 + in( i[ 2 ], j[ 2 ] ) * u2 + in( i[ 3 ], j[ 2 ] ) * u3 ) * v2
1610 + ( in( i[ 0 ], j[ 3 ] ) * u0 + in( i[ 1 ], j[ 3 ] ) * u1 + in( i[ 2 ], j[ 3 ] ) * u2 + in( i[ 3 ], j[ 3 ] ) * u3 ) * v3;
1611
1612 double min = type_limits< value_type >::minimum( );
1613 double max = type_limits< value_type >::maximum( );
1614 pix = pix > min ? pix : min;
1615 pix = pix < max ? pix : max;
1616
1617 return( pix );
1618 }
1619
1620 template < class T, class Allocator >
1621 static double interpolate( const array3< T, Allocator > &in,
1622 typename array3< T, Allocator >::size_type i[4],
1623 typename array3< T, Allocator >::size_type j[4],
1624 typename array3< T, Allocator >::size_type k[4],
1625 double x, double y, double z )
1626 {
1627 typedef typename array3< T, Allocator >::value_type value_type;
1628
1629 double u0 = bspline2( 1 + x );
1630 double u1 = bspline1( x );
1631 double u2 = bspline1( 1 - x );
1632 double u3 = bspline2( 2 - x );
1633 double v0 = bspline2( 1 + y );
1634 double v1 = bspline1( y );
1635 double v2 = bspline1( 1 - y );
1636 double v3 = bspline2( 2 - y );
1637 double w0 = bspline2( 1 + z );
1638 double w1 = bspline1( z );
1639 double w2 = bspline1( 1 - z );
1640 double w3 = bspline2( 2 - z );
1641
1642 double p0 = ( ( in( i[ 0 ], j[ 0 ], k[ 0 ] ) * u0 + in( i[ 1 ], j[ 0 ], k[ 0 ] ) * u1 + in( i[ 2 ], j[ 0 ], k[ 0 ] ) * u2 + in( i[ 3 ], j[ 0 ], k[ 0 ] ) * u3 ) * v0
1643 + ( in( i[ 0 ], j[ 1 ], k[ 0 ] ) * u0 + in( i[ 1 ], j[ 1 ], k[ 0 ] ) * u1 + in( i[ 2 ], j[ 1 ], k[ 0 ] ) * u2 + in( i[ 3 ], j[ 1 ], k[ 0 ] ) * u3 ) * v1
1644 + ( in( i[ 0 ], j[ 2 ], k[ 0 ] ) * u0 + in( i[ 1 ], j[ 2 ], k[ 0 ] ) * u1 + in( i[ 2 ], j[ 2 ], k[ 0 ] ) * u2 + in( i[ 3 ], j[ 2 ], k[ 0 ] ) * u3 ) * v2
1645 + ( in( i[ 0 ], j[ 3 ], k[ 0 ] ) * u0 + in( i[ 1 ], j[ 3 ], k[ 0 ] ) * u1 + in( i[ 2 ], j[ 3 ], k[ 0 ] ) * u2 + in( i[ 3 ], j[ 3 ], k[ 0 ] ) * u3 ) * v3 );
1646 double p1 = ( ( in( i[ 0 ], j[ 0 ], k[ 1 ] ) * u0 + in( i[ 1 ], j[ 0 ], k[ 1 ] ) * u1 + in( i[ 2 ], j[ 0 ], k[ 1 ] ) * u2 + in( i[ 3 ], j[ 0 ], k[ 1 ] ) * u3 ) * v0
1647 + ( in( i[ 0 ], j[ 1 ], k[ 1 ] ) * u0 + in( i[ 1 ], j[ 1 ], k[ 1 ] ) * u1 + in( i[ 2 ], j[ 1 ], k[ 1 ] ) * u2 + in( i[ 3 ], j[ 1 ], k[ 1 ] ) * u3 ) * v1
1648 + ( in( i[ 0 ], j[ 2 ], k[ 1 ] ) * u0 + in( i[ 1 ], j[ 2 ], k[ 1 ] ) * u1 + in( i[ 2 ], j[ 2 ], k[ 1 ] ) * u2 + in( i[ 3 ], j[ 2 ], k[ 1 ] ) * u3 ) * v2
1649 + ( in( i[ 0 ], j[ 3 ], k[ 1 ] ) * u0 + in( i[ 1 ], j[ 3 ], k[ 1 ] ) * u1 + in( i[ 2 ], j[ 3 ], k[ 1 ] ) * u2 + in( i[ 3 ], j[ 3 ], k[ 1 ] ) * u3 ) * v3 );
1650 double p2 = ( ( in( i[ 0 ], j[ 0 ], k[ 2 ] ) * u0 + in( i[ 1 ], j[ 0 ], k[ 2 ] ) * u1 + in( i[ 2 ], j[ 0 ], k[ 2 ] ) * u2 + in( i[ 3 ], j[ 0 ], k[ 2 ] ) * u3 ) * v0
1651 + ( in( i[ 0 ], j[ 1 ], k[ 2 ] ) * u0 + in( i[ 1 ], j[ 1 ], k[ 2 ] ) * u1 + in( i[ 2 ], j[ 1 ], k[ 2 ] ) * u2 + in( i[ 3 ], j[ 1 ], k[ 2 ] ) * u3 ) * v1
1652 + ( in( i[ 0 ], j[ 2 ], k[ 2 ] ) * u0 + in( i[ 1 ], j[ 2 ], k[ 2 ] ) * u1 + in( i[ 2 ], j[ 2 ], k[ 2 ] ) * u2 + in( i[ 3 ], j[ 2 ], k[ 2 ] ) * u3 ) * v2
1653 + ( in( i[ 0 ], j[ 3 ], k[ 2 ] ) * u0 + in( i[ 1 ], j[ 3 ], k[ 2 ] ) * u1 + in( i[ 2 ], j[ 3 ], k[ 2 ] ) * u2 + in( i[ 3 ], j[ 3 ], k[ 2 ] ) * u3 ) * v3 );
1654 double p3 = ( ( in( i[ 0 ], j[ 0 ], k[ 3 ] ) * u0 + in( i[ 1 ], j[ 0 ], k[ 3 ] ) * u1 + in( i[ 2 ], j[ 0 ], k[ 3 ] ) * u2 + in( i[ 3 ], j[ 0 ], k[ 3 ] ) * u3 ) * v0
1655 + ( in( i[ 0 ], j[ 1 ], k[ 3 ] ) * u0 + in( i[ 1 ], j[ 1 ], k[ 3 ] ) * u1 + in( i[ 2 ], j[ 1 ], k[ 3 ] ) * u2 + in( i[ 3 ], j[ 1 ], k[ 3 ] ) * u3 ) * v1
1656 + ( in( i[ 0 ], j[ 2 ], k[ 3 ] ) * u0 + in( i[ 1 ], j[ 2 ], k[ 3 ] ) * u1 + in( i[ 2 ], j[ 2 ], k[ 3 ] ) * u2 + in( i[ 3 ], j[ 2 ], k[ 3 ] ) * u3 ) * v2
1657 + ( in( i[ 0 ], j[ 3 ], k[ 3 ] ) * u0 + in( i[ 1 ], j[ 3 ], k[ 3 ] ) * u1 + in( i[ 2 ], j[ 3 ], k[ 3 ] ) * u2 + in( i[ 3 ], j[ 3 ], k[ 3 ] ) * u3 ) * v3 );
1658 double pix = p0 * w0 + p1 * w1 + p2 * w2 + p3 * w3;
1659
1660 double min = type_limits< value_type >::minimum( );
1661 double max = type_limits< value_type >::maximum( );
1662 pix = pix > min ? pix : min;
1663 pix = pix < max ? pix : max;
1664
1665 return( pix );
1666 }
1667 };
1668
1669 template < >
1670 struct _bspline_< true >
1671 {
1672 template < class T, class Allocator >
1673 static const typename T::template rebind< double >::other interpolate( const array< T, Allocator > &in,
1674 typename array< T, Allocator >::size_type i[4],
1675 typename array< T, Allocator >::size_type /* j */[4],
1676 typename array< T, Allocator >::size_type /* k */[4],
1677 double x, double /* y */, double /* z */ )
1678 {
1679 typedef typename array< T, Allocator >::value_type color;
1680 typedef typename color::value_type value_type;
1681 typedef typename T::template rebind< double >::other ovalue_type;
1682
1683 double u0 = bspline2( 1 + x );
1684 double u1 = bspline1( x );
1685 double u2 = bspline1( 1 - x );
1686 double u3 = bspline2( 2 - x );
1687 double r = in[ i[ 0 ] ].r * u0 + in[ i[ 1 ] ].r * u1 + in[ i[ 2 ] ].r * u2 + in[ i[ 3 ] ].r * u3;
1688 double g = in[ i[ 0 ] ].g * u0 + in[ i[ 1 ] ].g * u1 + in[ i[ 2 ] ].g * u2 + in[ i[ 3 ] ].g * u3;
1689 double b = in[ i[ 0 ] ].b * u0 + in[ i[ 1 ] ].b * u1 + in[ i[ 2 ] ].b * u2 + in[ i[ 3 ] ].b * u3;
1690
1691 double min = type_limits< value_type >::minimum( );
1692 double max = type_limits< value_type >::maximum( );
1693
1694 r = r > min ? r : min;
1695 r = r < max ? r : max;
1696 g = g > min ? g : min;
1697 g = g < max ? g : max;
1698 b = b > min ? b : min;
1699 b = b < max ? b : max;
1700
1701 return( ovalue_type( r, g, b ) );
1702 }
1703
1704 template < class T, class Allocator >
1705 static const typename T::template rebind< double >::other interpolate( const array2< T, Allocator > &in,
1706 typename array2< T, Allocator >::size_type i[4],
1707 typename array2< T, Allocator >::size_type j[4],
1708 typename array2< T, Allocator >::size_type /* k */[4],
1709 double x, double y, double /* z */ )
1710 {
1711 typedef typename array2< T, Allocator >::value_type color;
1712 typedef typename color::value_type value_type;
1713 typedef typename T::template rebind< double >::other ovalue_type;
1714
1715 double u0 = bspline2( 1 + x );
1716 double u1 = bspline1( x );
1717 double u2 = bspline1( 1 - x );
1718 double u3 = bspline2( 2 - x );
1719 double v0 = bspline2( 1 + y );
1720 double v1 = bspline1( y );
1721 double v2 = bspline1( 1 - y );
1722 double v3 = bspline2( 2 - y );
1723
1724 double r = ( in( i[ 0 ], j[ 0 ] ).r * u0 + in( i[ 1 ], j[ 0 ] ).r * u1 + in( i[ 2 ], j[ 0 ] ).r * u2 + in( i[ 3 ], j[ 0 ] ).r * u3 ) * v0
1725 + ( in( i[ 0 ], j[ 1 ] ).r * u0 + in( i[ 1 ], j[ 1 ] ).r * u1 + in( i[ 2 ], j[ 1 ] ).r * u2 + in( i[ 3 ], j[ 1 ] ).r * u3 ) * v1
1726 + ( in( i[ 0 ], j[ 2 ] ).r * u0 + in( i[ 1 ], j[ 2 ] ).r * u1 + in( i[ 2 ], j[ 2 ] ).r * u2 + in( i[ 3 ], j[ 2 ] ).r * u3 ) * v2
1727 + ( in( i[ 0 ], j[ 3 ] ).r * u0 + in( i[ 1 ], j[ 3 ] ).r * u1 + in( i[ 2 ], j[ 3 ] ).r * u2 + in( i[ 3 ], j[ 3 ] ).r * u3 ) * v3;
1728 double g = ( in( i[ 0 ], j[ 0 ] ).g * u0 + in( i[ 1 ], j[ 0 ] ).g * u1 + in( i[ 2 ], j[ 0 ] ).g * u2 + in( i[ 3 ], j[ 0 ] ).g * u3 ) * v0
1729 + ( in( i[ 0 ], j[ 1 ] ).g * u0 + in( i[ 1 ], j[ 1 ] ).g * u1 + in( i[ 2 ], j[ 1 ] ).g * u2 + in( i[ 3 ], j[ 1 ] ).g * u3 ) * v1
1730 + ( in( i[ 0 ], j[ 2 ] ).g * u0 + in( i[ 1 ], j[ 2 ] ).g * u1 + in( i[ 2 ], j[ 2 ] ).g * u2 + in( i[ 3 ], j[ 2 ] ).g * u3 ) * v2
1731 + ( in( i[ 0 ], j[ 3 ] ).g * u0 + in( i[ 1 ], j[ 3 ] ).g * u1 + in( i[ 2 ], j[ 3 ] ).g * u2 + in( i[ 3 ], j[ 3 ] ).g * u3 ) * v3;
1732 double b = ( in( i[ 0 ], j[ 0 ] ).b * u0 + in( i[ 1 ], j[ 0 ] ).b * u1 + in( i[ 2 ], j[ 0 ] ).b * u2 + in( i[ 3 ], j[ 0 ] ).b * u3 ) * v0
1733 + ( in( i[ 0 ], j[ 1 ] ).b * u0 + in( i[ 1 ], j[ 1 ] ).b * u1 + in( i[ 2 ], j[ 1 ] ).b * u2 + in( i[ 3 ], j[ 1 ] ).b * u3 ) * v1
1734 + ( in( i[ 0 ], j[ 2 ] ).b * u0 + in( i[ 1 ], j[ 2 ] ).b * u1 + in( i[ 2 ], j[ 2 ] ).b * u2 + in( i[ 3 ], j[ 2 ] ).b * u3 ) * v2
1735 + ( in( i[ 0 ], j[ 3 ] ).b * u0 + in( i[ 1 ], j[ 3 ] ).b * u1 + in( i[ 2 ], j[ 3 ] ).b * u2 + in( i[ 3 ], j[ 3 ] ).b * u3 ) * v3;
1736
1737 double min = type_limits< value_type >::minimum( );
1738 double max = type_limits< value_type >::maximum( );
1739
1740 r = r > min ? r : min;
1741 r = r < max ? r : max;
1742 g = g > min ? g : min;
1743 g = g < max ? g : max;
1744 b = b > min ? b : min;
1745 b = b < max ? b : max;
1746
1747 return( ovalue_type( r, g, b ) );
1748 }
1749
1750 template < class T, class Allocator >
1751 static const typename T::template rebind< double >::other interpolate( const array3< T, Allocator > &in,
1752 typename array3< T, Allocator >::size_type i[4],
1753 typename array3< T, Allocator >::size_type j[4],
1754 typename array3< T, Allocator >::size_type k[4],
1755 double x, double y, double z )
1756 {
1757 typedef typename array3< T, Allocator >::value_type color;
1758 typedef typename color::value_type value_type;
1759 typedef typename T::template rebind< double >::other ovalue_type;
1760
1761 double u0 = bspline2( 1 + x );
1762 double u1 = bspline1( x );
1763 double u2 = bspline1( 1 - x );
1764 double u3 = bspline2( 2 - x );
1765 double v0 = bspline2( 1 + y );
1766 double v1 = bspline1( y );
1767 double v2 = bspline1( 1 - y );
1768 double v3 = bspline2( 2 - y );
1769 double w0 = bspline2( 1 + z );
1770 double w1 = bspline1( z );
1771 double w2 = bspline1( 1 - z );
1772 double w3 = bspline2( 2 - z );
1773
1774 double r0 = ( ( in( i[ 0 ], j[ 0 ], k[ 0 ] ).r * u0 + in( i[ 1 ], j[ 0 ], k[ 0 ] ).r * u1 + in( i[ 2 ], j[ 0 ], k[ 0 ] ).r * u2 + in( i[ 3 ], j[ 0 ], k[ 0 ] ).r * u3 ) * v0
1775 + ( in( i[ 0 ], j[ 1 ], k[ 0 ] ).r * u0 + in( i[ 1 ], j[ 1 ], k[ 0 ] ).r * u1 + in( i[ 2 ], j[ 1 ], k[ 0 ] ).r * u2 + in( i[ 3 ], j[ 1 ], k[ 0 ] ).r * u3 ) * v1
1776 + ( in( i[ 0 ], j[ 2 ], k[ 0 ] ).r * u0 + in( i[ 1 ], j[ 2 ], k[ 0 ] ).r * u1 + in( i[ 2 ], j[ 2 ], k[ 0 ] ).r * u2 + in( i[ 3 ], j[ 2 ], k[ 0 ] ).r * u3 ) * v2
1777 + ( in( i[ 0 ], j[ 3 ], k[ 0 ] ).r * u0 + in( i[ 1 ], j[ 3 ], k[ 0 ] ).r * u1 + in( i[ 2 ], j[ 3 ], k[ 0 ] ).r * u2 + in( i[ 3 ], j[ 3 ], k[ 0 ] ).r * u3 ) * v3 );
1778 double r1 = ( ( in( i[ 0 ], j[ 0 ], k[ 1 ] ).r * u0 + in( i[ 1 ], j[ 0 ], k[ 1 ] ).r * u1 + in( i[ 2 ], j[ 0 ], k[ 1 ] ).r * u2 + in( i[ 3 ], j[ 0 ], k[ 1 ] ).r * u3 ) * v0
1779 + ( in( i[ 0 ], j[ 1 ], k[ 1 ] ).r * u0 + in( i[ 1 ], j[ 1 ], k[ 1 ] ).r * u1 + in( i[ 2 ], j[ 1 ], k[ 1 ] ).r * u2 + in( i[ 3 ], j[ 1 ], k[ 1 ] ).r * u3 ) * v1
1780 + ( in( i[ 0 ], j[ 2 ], k[ 1 ] ).r * u0 + in( i[ 1 ], j[ 2 ], k[ 1 ] ).r * u1 + in( i[ 2 ], j[ 2 ], k[ 1 ] ).r * u2 + in( i[ 3 ], j[ 2 ], k[ 1 ] ).r * u3 ) * v2
1781 + ( in( i[ 0 ], j[ 3 ], k[ 1 ] ).r * u0 + in( i[ 1 ], j[ 3 ], k[ 1 ] ).r * u1 + in( i[ 2 ], j[ 3 ], k[ 1 ] ).r * u2 + in( i[ 3 ], j[ 3 ], k[ 1 ] ).r * u3 ) * v3 );
1782 double r2 = ( ( in( i[ 0 ], j[ 0 ], k[ 2 ] ).r * u0 + in( i[ 1 ], j[ 0 ], k[ 2 ] ).r * u1 + in( i[ 2 ], j[ 0 ], k[ 2 ] ).r * u2 + in( i[ 3 ], j[ 0 ], k[ 2 ] ).r * u3 ) * v0
1783 + ( in( i[ 0 ], j[ 1 ], k[ 2 ] ).r * u0 + in( i[ 1 ], j[ 1 ], k[ 2 ] ).r * u1 + in( i[ 2 ], j[ 1 ], k[ 2 ] ).r * u2 + in( i[ 3 ], j[ 1 ], k[ 2 ] ).r * u3 ) * v1
1784 + ( in( i[ 0 ], j[ 2 ], k[ 2 ] ).r * u0 + in( i[ 1 ], j[ 2 ], k[ 2 ] ).r * u1 + in( i[ 2 ], j[ 2 ], k[ 2 ] ).r * u2 + in( i[ 3 ], j[ 2 ], k[ 2 ] ).r * u3 ) * v2
1785 + ( in( i[ 0 ], j[ 3 ], k[ 2 ] ).r * u0 + in( i[ 1 ], j[ 3 ], k[ 2 ] ).r * u1 + in( i[ 2 ], j[ 3 ], k[ 2 ] ).r * u2 + in( i[ 3 ], j[ 3 ], k[ 2 ] ).r * u3 ) * v3 );
1786 double r3 = ( ( in( i[ 0 ], j[ 0 ], k[ 3 ] ).r * u0 + in( i[ 1 ], j[ 0 ], k[ 3 ] ).r * u1 + in( i[ 2 ], j[ 0 ], k[ 3 ] ).r * u2 + in( i[ 3 ], j[ 0 ], k[ 3 ] ).r * u3 ) * v0
1787 + ( in( i[ 0 ], j[ 1 ], k[ 3 ] ).r * u0 + in( i[ 1 ], j[ 1 ], k[ 3 ] ).r * u1 + in( i[ 2 ], j[ 1 ], k[ 3 ] ).r * u2 + in( i[ 3 ], j[ 1 ], k[ 3 ] ).r * u3 ) * v1
1788 + ( in( i[ 0 ], j[ 2 ], k[ 3 ] ).r * u0 + in( i[ 1 ], j[ 2 ], k[ 3 ] ).r * u1 + in( i[ 2 ], j[ 2 ], k[ 3 ] ).r * u2 + in( i[ 3 ], j[ 2 ], k[ 3 ] ).r * u3 ) * v2
1789 + ( in( i[ 0 ], j[ 3 ], k[ 3 ] ).r * u0 + in( i[ 1 ], j[ 3 ], k[ 3 ] ).r * u1 + in( i[ 2 ], j[ 3 ], k[ 3 ] ).r * u2 + in( i[ 3 ], j[ 3 ], k[ 3 ] ).r * u3 ) * v3 );
1790 double g0 = ( ( in( i[ 0 ], j[ 0 ], k[ 0 ] ).g * u0 + in( i[ 1 ], j[ 0 ], k[ 0 ] ).g * u1 + in( i[ 2 ], j[ 0 ], k[ 0 ] ).g * u2 + in( i[ 3 ], j[ 0 ], k[ 0 ] ).g * u3 ) * v0
1791 + ( in( i[ 0 ], j[ 1 ], k[ 0 ] ).g * u0 + in( i[ 1 ], j[ 1 ], k[ 0 ] ).g * u1 + in( i[ 2 ], j[ 1 ], k[ 0 ] ).g * u2 + in( i[ 3 ], j[ 1 ], k[ 0 ] ).g * u3 ) * v1
1792 + ( in( i[ 0 ], j[ 2 ], k[ 0 ] ).g * u0 + in( i[ 1 ], j[ 2 ], k[ 0 ] ).g * u1 + in( i[ 2 ], j[ 2 ], k[ 0 ] ).g * u2 + in( i[ 3 ], j[ 2 ], k[ 0 ] ).g * u3 ) * v2
1793 + ( in( i[ 0 ], j[ 3 ], k[ 0 ] ).g * u0 + in( i[ 1 ], j[ 3 ], k[ 0 ] ).g * u1 + in( i[ 2 ], j[ 3 ], k[ 0 ] ).g * u2 + in( i[ 3 ], j[ 3 ], k[ 0 ] ).g * u3 ) * v3 );
1794 double g1 = ( ( in( i[ 0 ], j[ 0 ], k[ 1 ] ).g * u0 + in( i[ 1 ], j[ 0 ], k[ 1 ] ).g * u1 + in( i[ 2 ], j[ 0 ], k[ 1 ] ).g * u2 + in( i[ 3 ], j[ 0 ], k[ 1 ] ).g * u3 ) * v0
1795 + ( in( i[ 0 ], j[ 1 ], k[ 1 ] ).g * u0 + in( i[ 1 ], j[ 1 ], k[ 1 ] ).g * u1 + in( i[ 2 ], j[ 1 ], k[ 1 ] ).g * u2 + in( i[ 3 ], j[ 1 ], k[ 1 ] ).g * u3 ) * v1
1796 + ( in( i[ 0 ], j[ 2 ], k[ 1 ] ).g * u0 + in( i[ 1 ], j[ 2 ], k[ 1 ] ).g * u1 + in( i[ 2 ], j[ 2 ], k[ 1 ] ).g * u2 + in( i[ 3 ], j[ 2 ], k[ 1 ] ).g * u3 ) * v2
1797 + ( in( i[ 0 ], j[ 3 ], k[ 1 ] ).g * u0 + in( i[ 1 ], j[ 3 ], k[ 1 ] ).g * u1 + in( i[ 2 ], j[ 3 ], k[ 1 ] ).g * u2 + in( i[ 3 ], j[ 3 ], k[ 1 ] ).g * u3 ) * v3 );
1798 double g2 = ( ( in( i[ 0 ], j[ 0 ], k[ 2 ] ).g * u0 + in( i[ 1 ], j[ 0 ], k[ 2 ] ).g * u1 + in( i[ 2 ], j[ 0 ], k[ 2 ] ).g * u2 + in( i[ 3 ], j[ 0 ], k[ 2 ] ).g * u3 ) * v0
1799 + ( in( i[ 0 ], j[ 1 ], k[ 2 ] ).g * u0 + in( i[ 1 ], j[ 1 ], k[ 2 ] ).g * u1 + in( i[ 2 ], j[ 1 ], k[ 2 ] ).g * u2 + in( i[ 3 ], j[ 1 ], k[ 2 ] ).g * u3 ) * v1
1800 + ( in( i[ 0 ], j[ 2 ], k[ 2 ] ).g * u0 + in( i[ 1 ], j[ 2 ], k[ 2 ] ).g * u1 + in( i[ 2 ], j[ 2 ], k[ 2 ] ).g * u2 + in( i[ 3 ], j[ 2 ], k[ 2 ] ).g * u3 ) * v2
1801 + ( in( i[ 0 ], j[ 3 ], k[ 2 ] ).g * u0 + in( i[ 1 ], j[ 3 ], k[ 2 ] ).g * u1 + in( i[ 2 ], j[ 3 ], k[ 2 ] ).g * u2 + in( i[ 3 ], j[ 3 ], k[ 2 ] ).g * u3 ) * v3 );
1802 double g3 = ( ( in( i[ 0 ], j[ 0 ], k[ 3 ] ).g * u0 + in( i[ 1 ], j[ 0 ], k[ 3 ] ).g * u1 + in( i[ 2 ], j[ 0 ], k[ 3 ] ).g * u2 + in( i[ 3 ], j[ 0 ], k[ 3 ] ).g * u3 ) * v0
1803 + ( in( i[ 0 ], j[ 1 ], k[ 3 ] ).g * u0 + in( i[ 1 ], j[ 1 ], k[ 3 ] ).g * u1 + in( i[ 2 ], j[ 1 ], k[ 3 ] ).g * u2 + in( i[ 3 ], j[ 1 ], k[ 3 ] ).g * u3 ) * v1
1804 + ( in( i[ 0 ], j[ 2 ], k[ 3 ] ).g * u0 + in( i[ 1 ], j[ 2 ], k[ 3 ] ).g * u1 + in( i[ 2 ], j[ 2 ], k[ 3 ] ).g * u2 + in( i[ 3 ], j[ 2 ], k[ 3 ] ).g * u3 ) * v2
1805 + ( in( i[ 0 ], j[ 3 ], k[ 3 ] ).g * u0 + in( i[ 1 ], j[ 3 ], k[ 3 ] ).g * u1 + in( i[ 2 ], j[ 3 ], k[ 3 ] ).g * u2 + in( i[ 3 ], j[ 3 ], k[ 3 ] ).g * u3 ) * v3 );
1806 double b0 = ( ( in( i[ 0 ], j[ 0 ], k[ 0 ] ).b * u0 + in( i[ 1 ], j[ 0 ], k[ 0 ] ).b * u1 + in( i[ 2 ], j[ 0 ], k[ 0 ] ).b * u2 + in( i[ 3 ], j[ 0 ], k[ 0 ] ).b * u3 ) * v0
1807 + ( in( i[ 0 ], j[ 1 ], k[ 0 ] ).b * u0 + in( i[ 1 ], j[ 1 ], k[ 0 ] ).b * u1 + in( i[ 2 ], j[ 1 ], k[ 0 ] ).b * u2 + in( i[ 3 ], j[ 1 ], k[ 0 ] ).b * u3 ) * v1
1808 + ( in( i[ 0 ], j[ 2 ], k[ 0 ] ).b * u0 + in( i[ 1 ], j[ 2 ], k[ 0 ] ).b * u1 + in( i[ 2 ], j[ 2 ], k[ 0 ] ).b * u2 + in( i[ 3 ], j[ 2 ], k[ 0 ] ).b * u3 ) * v2
1809 + ( in( i[ 0 ], j[ 3 ], k[ 0 ] ).b * u0 + in( i[ 1 ], j[ 3 ], k[ 0 ] ).b * u1 + in( i[ 2 ], j[ 3 ], k[ 0 ] ).b * u2 + in( i[ 3 ], j[ 3 ], k[ 0 ] ).b * u3 ) * v3 );
1810 double b1 = ( ( in( i[ 0 ], j[ 0 ], k[ 1 ] ).b * u0 + in( i[ 1 ], j[ 0 ], k[ 1 ] ).b * u1 + in( i[ 2 ], j[ 0 ], k[ 1 ] ).b * u2 + in( i[ 3 ], j[ 0 ], k[ 1 ] ).b * u3 ) * v0
1811 + ( in( i[ 0 ], j[ 1 ], k[ 1 ] ).b * u0 + in( i[ 1 ], j[ 1 ], k[ 1 ] ).b * u1 + in( i[ 2 ], j[ 1 ], k[ 1 ] ).b * u2 + in( i[ 3 ], j[ 1 ], k[ 1 ] ).b * u3 ) * v1
1812 + ( in( i[ 0 ], j[ 2 ], k[ 1 ] ).b * u0 + in( i[ 1 ], j[ 2 ], k[ 1 ] ).b * u1 + in( i[ 2 ], j[ 2 ], k[ 1 ] ).b * u2 + in( i[ 3 ], j[ 2 ], k[ 1 ] ).b * u3 ) * v2
1813 + ( in( i[ 0 ], j[ 3 ], k[ 1 ] ).b * u0 + in( i[ 1 ], j[ 3 ], k[ 1 ] ).b * u1 + in( i[ 2 ], j[ 3 ], k[ 1 ] ).b * u2 + in( i[ 3 ], j[ 3 ], k[ 1 ] ).b * u3 ) * v3 );
1814 double b2 = ( ( in( i[ 0 ], j[ 0 ], k[ 2 ] ).b * u0 + in( i[ 1 ], j[ 0 ], k[ 2 ] ).b * u1 + in( i[ 2 ], j[ 0 ], k[ 2 ] ).b * u2 + in( i[ 3 ], j[ 0 ], k[ 2 ] ).b * u3 ) * v0
1815 + ( in( i[ 0 ], j[ 1 ], k[ 2 ] ).b * u0 + in( i[ 1 ], j[ 1 ], k[ 2 ] ).b * u1 + in( i[ 2 ], j[ 1 ], k[ 2 ] ).b * u2 + in( i[ 3 ], j[ 1 ], k[ 2 ] ).b * u3 ) * v1
1816 + ( in( i[ 0 ], j[ 2 ], k[ 2 ] ).b * u0 + in( i[ 1 ], j[ 2 ], k[ 2 ] ).b * u1 + in( i[ 2 ], j[ 2 ], k[ 2 ] ).b * u2 + in( i[ 3 ], j[ 2 ], k[ 2 ] ).b * u3 ) * v2
1817 + ( in( i[ 0 ], j[ 3 ], k[ 2 ] ).b * u0 + in( i[ 1 ], j[ 3 ], k[ 2 ] ).b * u1 + in( i[ 2 ], j[ 3 ], k[ 2 ] ).b * u2 + in( i[ 3 ], j[ 3 ], k[ 2 ] ).b * u3 ) * v3 );
1818 double b3 = ( ( in( i[ 0 ], j[ 0 ], k[ 3 ] ).b * u0 + in( i[ 1 ], j[ 0 ], k[ 3 ] ).b * u1 + in( i[ 2 ], j[ 0 ], k[ 3 ] ).b * u2 + in( i[ 3 ], j[ 0 ], k[ 3 ] ).b * u3 ) * v0
1819 + ( in( i[ 0 ], j[ 1 ], k[ 3 ] ).b * u0 + in( i[ 1 ], j[ 1 ], k[ 3 ] ).b * u1 + in( i[ 2 ], j[ 1 ], k[ 3 ] ).b * u2 + in( i[ 3 ], j[ 1 ], k[ 3 ] ).b * u3 ) * v1
1820 + ( in( i[ 0 ], j[ 2 ], k[ 3 ] ).b * u0 + in( i[ 1 ], j[ 2 ], k[ 3 ] ).b * u1 + in( i[ 2 ], j[ 2 ], k[ 3 ] ).b * u2 + in( i[ 3 ], j[ 2 ], k[ 3 ] ).b * u3 ) * v2
1821 + ( in( i[ 0 ], j[ 3 ], k[ 3 ] ).b * u0 + in( i[ 1 ], j[ 3 ], k[ 3 ] ).b * u1 + in( i[ 2 ], j[ 3 ], k[ 3 ] ).b * u2 + in( i[ 3 ], j[ 3 ], k[ 3 ] ).b * u3 ) * v3 );
1822
1823 double r = r0 * w0 + r1 * w1 + r2 * w2 + r3 * w3;
1824 double g = g0 * w0 + g1 * w1 + g2 * w2 + g3 * w3;
1825 double b = b0 * w0 + b1 * w1 + b2 * w2 + b3 * w3;
1826
1827 double min = type_limits< value_type >::minimum( );
1828 double max = type_limits< value_type >::maximum( );
1829
1830 r = r > min ? r : min;
1831 r = r < max ? r : max;
1832 g = g > min ? g : min;
1833 g = g < max ? g : max;
1834 b = b > min ? b : min;
1835 b = b < max ? b : max;
1836
1837 return( ovalue_type( r, g, b ) );
1838 }
1839 };
1840
1841 template < class Array1, class Array2 >
1842 void interpolate( const Array1 &in, Array2 &out,
1843 typename Array1::size_type thread_idx, typename Array1::size_type thread_numx,
1844 typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
1845 typename Array1::size_type thread_idz, typename Array1::size_type thread_numz )
1846 {
1847 typedef typename Array1::size_type size_type;
1848 typedef typename Array1::value_type value_type;
1849 typedef typename Array2::value_type out_value_type;
1850
1851 size_type i, j, k, ii[ 4 ], jj[ 4 ], kk[ 4 ];
1852 size_type iw = in.width( );
1853 size_type ih = in.height( );
1854 size_type id = in.depth( );
1855 size_type ow = out.width( );
1856 size_type oh = out.height( );
1857 size_type od = out.depth( );
1858
1859 double sx = static_cast< double >( iw ) / static_cast< double >( ow );
1860 double sy = static_cast< double >( ih ) / static_cast< double >( oh );
1861 double sz = static_cast< double >( id ) / static_cast< double >( od );
1862 double x, y, z;
1863
1864 for( k = thread_idz ; k < od ; k += thread_numz )
1865 {
1866 z = sz * k;
1867 kk[ 1 ] = static_cast< size_type >( z );
1868 kk[ 0 ] = kk[ 1 ] > 0 ? kk[ 1 ] - 1 : kk[ 1 ];
1869 kk[ 2 ] = kk[ 1 ] < id - 1 ? kk[ 1 ] + 1 : kk[ 1 ];
1870 kk[ 3 ] = kk[ 2 ] < id - 1 ? kk[ 2 ] + 1 : kk[ 2 ];
1871 z -= kk[ 1 ];
1872
1873 for( j = thread_idy ; j < oh ; j += thread_numy )
1874 {
1875 y = sy * j;
1876 jj[ 1 ] = static_cast< size_type >( y );
1877 jj[ 0 ] = jj[ 1 ] > 0 ? jj[ 1 ] - 1 : jj[ 1 ];
1878 jj[ 2 ] = jj[ 1 ] < ih - 1 ? jj[ 1 ] + 1 : jj[ 1 ];
1879 jj[ 3 ] = jj[ 2 ] < ih - 1 ? jj[ 2 ] + 1 : jj[ 2 ];
1880 y -= jj[ 1 ];
1881
1882 for( i = thread_idx ; i < ow ; i += thread_numx )
1883 {
1884 x = sx * i;
1885 ii[ 1 ] = static_cast< size_type >( x );
1886 ii[ 0 ] = ii[ 1 ] > 0 ? ii[ 1 ] - 1 : ii[ 1 ];
1887 ii[ 2 ] = ii[ 1 ] < iw - 1 ? ii[ 1 ] + 1 : ii[ 1 ];
1888 ii[ 3 ] = ii[ 2 ] < iw - 1 ? ii[ 2 ] + 1 : ii[ 2 ];
1889 x -= ii[ 1 ];
1890
1891 __interpolate_utility__::round( _bspline_< is_color< value_type >::value >::interpolate( in, ii, jj, kk, x, y, z ), out( i, j, k ) );
1892 }
1893 }
1894 }
1895 }
1896 }
1897
1898
1899
1900 // フルスケールでのsinc関数を用いた内挿補間
1901 namespace __sinc__
1902 {
1903 // 厳密な定義での sinc 関数
1904 inline double sinc( double t )
1905 {
1906 return( t == 0.0 ? 1.0 : std::sin( 3.1415926535897932384626433832795 * t ) / ( 3.1415926535897932384626433832795 * t ) );
1907 }
1908
1909 template < bool b >
1910 struct _sinc_
1911 {
1912 template < class T, class Allocator >
1913 static double interpolate( const array< T, Allocator > &in,
1914 typename array< T, Allocator >::difference_type ix,
1915 typename array< T, Allocator >::difference_type iy,
1916 typename array< T, Allocator >::difference_type iz,
1917 double x, double y, double z )
1918 {
1919 typedef typename array< T, Allocator >::difference_type difference_type;
1920 typedef typename array< T, Allocator >::value_type value_type;
1921
1922 double pix = 0.0;
1923 difference_type i, num = in.size( );
1924 for( i = 0 ; i <= ix ; i++ )
1925 {
1926 pix += in[ i ] * sinc( ix - i - x );
1927 }
1928 for( ; i < num ; i++ )
1929 {
1930 pix += in[ i ] * sinc( i - ix + x );
1931 }
1932
1933 double min = type_limits< value_type >::minimum( );
1934 double max = type_limits< value_type >::maximum( );
1935 pix = pix > min ? pix : min;
1936 pix = pix < max ? pix : max;
1937
1938 return( pix );
1939 }
1940
1941 template < class T, class Allocator >
1942 static double interpolate( const array2< T, Allocator > &in,
1943 typename array2< T, Allocator >::difference_type ix,
1944 typename array2< T, Allocator >::difference_type iy,
1945 typename array2< T, Allocator >::difference_type /* iz */,
1946 double x, double y, double /* z */ )
1947 {
1948 typedef typename array2< T, Allocator >::difference_type difference_type;
1949 typedef typename array2< T, Allocator >::value_type value_type;
1950 typedef typename array2< T, Allocator >::const_pointer const_pointer;
1951
1952 double pix = 0.0;
1953 difference_type width = in.width( );
1954 difference_type height = in.height( );
1955
1956 difference_type j;
1957 for( j = 0 ; j < iy ; j++ )
1958 {
1959 difference_type i;
1960 const_pointer p = &in( 0, j );
1961 double tmpX = 0.0;
1962 for( i = 0 ; i <= ix ; i++ )
1963 {
1964 tmpX += p[ i ] * sinc( ix - i - x );
1965 }
1966 for( ; i < width ; i++ )
1967 {
1968 tmpX += p[ i ] * sinc( i - ix + x );
1969 }
1970
1971 pix += tmpX * sinc( iy - j - y );
1972 }
1973 for( ; j < height ; j++ )
1974 {
1975 difference_type i;
1976 const_pointer p = &in( 0, j );
1977 double tmpX = 0.0;
1978 for( i = 0 ; i <= ix ; i++ )
1979 {
1980 tmpX += p[ i ] * sinc( ix - i - x );
1981 }
1982 for( ; i < width ; i++ )
1983 {
1984 tmpX += p[ i ] * sinc( i - ix + x );
1985 }
1986
1987 pix += tmpX * sinc( j - iy + y );
1988 }
1989
1990 double min = type_limits< value_type >::minimum( );
1991 double max = type_limits< value_type >::maximum( );
1992 pix = pix > min ? pix : min;
1993 pix = pix < max ? pix : max;
1994
1995 return( pix );
1996 }
1997
1998 template < class T, class Allocator >
1999 static double interpolate( const array3< T, Allocator > &in,
2000 typename array3< T, Allocator >::difference_type ix,
2001 typename array3< T, Allocator >::difference_type iy,
2002 typename array3< T, Allocator >::difference_type iz,
2003 double x, double y, double z )
2004 {
2005 typedef typename array2< T, Allocator >::difference_type difference_type;
2006 typedef typename array2< T, Allocator >::value_type value_type;
2007 typedef typename array2< T, Allocator >::const_pointer const_pointer;
2008
2009 double pix = 0.0;
2010 difference_type width = in.width( );
2011 difference_type height = in.height( );
2012 difference_type depth = in.depth( );
2013
2014 difference_type k;
2015 for( k = 0 ; k < iz ; k++ )
2016 {
2017 difference_type j;
2018 double tmpY = 0.0;
2019 for( j = 0 ; j < iy ; j++ )
2020 {
2021 difference_type i;
2022 const_pointer p = &in( 0, j, k );
2023 double tmpX = 0.0;
2024 for( i = 0 ; i <= ix ; i++ )
2025 {
2026 tmpX += p[ i ] * sinc( ix - i - x );
2027 }
2028 for( ; i < width ; i++ )
2029 {
2030 tmpX += p[ i ] * sinc( i - ix + x );
2031 }
2032
2033 tmpY += tmpX * sinc( iy - j - y );
2034 }
2035 for( ; j < height ; j++ )
2036 {
2037 difference_type i;
2038 const_pointer p = &in( 0, j, k );
2039 double tmpX = 0.0;
2040 for( i = 0 ; i <= ix ; i++ )
2041 {
2042 tmpX += p[ i ] * sinc( ix - i - x );
2043 }
2044 for( ; i < width ; i++ )
2045 {
2046 tmpX += p[ i ] * sinc( i - ix + x );
2047 }
2048
2049 tmpY += tmpX * sinc( j - iy + y );
2050 }
2051
2052 pix += tmpY * sinc( iz - k - z );
2053 }
2054 for( ; k < depth ; k++ )
2055 {
2056 difference_type j;
2057 double tmpY = 0.0;
2058 for( j = 0 ; j < iy ; j++ )
2059 {
2060 difference_type i;
2061 const_pointer p = &in( 0, j, k );
2062 double tmpX = 0.0;
2063 for( i = 0 ; i <= ix ; i++ )
2064 {
2065 tmpX += p[ i ] * sinc( ix - i - x );
2066 }
2067 for( ; i < width ; i++ )
2068 {
2069 tmpX += p[ i ] * sinc( i - ix + x );
2070 }
2071
2072 tmpY += tmpX * sinc( iy - j - y );
2073 }
2074 for( ; j < height ; j++ )
2075 {
2076 difference_type i;
2077 const_pointer p = &in( 0, j, k );
2078 double tmpX = 0.0;
2079 for( i = 0 ; i <= ix ; i++ )
2080 {
2081 tmpX += p[ i ] * sinc( ix - i - x );
2082 }
2083 for( ; i < width ; i++ )
2084 {
2085 tmpX += p[ i ] * sinc( i - ix + x );
2086 }
2087
2088 tmpY += tmpX * sinc( j - iy + y );
2089 }
2090
2091 pix += tmpY * sinc( k - iz + z );
2092 }
2093
2094 double min = type_limits< value_type >::minimum( );
2095 double max = type_limits< value_type >::maximum( );
2096 pix = pix > min ? pix : min;
2097 pix = pix < max ? pix : max;
2098
2099 return( pix );
2100 }
2101 };
2102
2103
2104 template < >
2105 struct _sinc_< true >
2106 {
2107 template < class T, class Allocator >
2108 static typename T::template rebind< double >::other interpolate( const array< T, Allocator > &in,
2109 typename array< T, Allocator >::difference_type ix,
2110 typename array< T, Allocator >::difference_type iy,
2111 typename array< T, Allocator >::difference_type iz,
2112 double /* x */, double /* y */, double /* z */ )
2113 {
2114 return( in( ix, iy, iz ) );
2115 }
2116
2117 template < class T, class Allocator >
2118 static typename T::template rebind< double >::other interpolate( const array2< T, Allocator > &in,
2119 typename array2< T, Allocator >::difference_type ix,
2120 typename array2< T, Allocator >::difference_type iy,
2121 typename array2< T, Allocator >::difference_type iz,
2122 double /* x */, double /* y */, double /* z */ )
2123 {
2124 return( in( ix, iy, iz ) );
2125 }
2126
2127 template < class T, class Allocator >
2128 static typename T::template rebind< double >::other interpolate( const array3< T, Allocator > &in,
2129 typename array3< T, Allocator >::difference_type ix,
2130 typename array3< T, Allocator >::difference_type iy,
2131 typename array3< T, Allocator >::difference_type iz,
2132 double /* x */, double /* y */, double /* z */ )
2133 {
2134 return( in( ix, iy, iz ) );
2135 }
2136 };
2137
2138
2139 template < class Array1, class Array2 >
2140 void interpolate( const Array1 &in, Array2 &out,
2141 typename Array1::size_type thread_idx, typename Array1::size_type thread_numx,
2142 typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
2143 typename Array1::size_type thread_idz, typename Array1::size_type thread_numz )
2144 {
2145 typedef typename Array1::size_type size_type;
2146 typedef typename Array1::difference_type difference_type;
2147 typedef typename Array1::value_type value_type;
2148 typedef typename Array2::value_type out_value_type;
2149
2150 size_type i, j, k;
2151 difference_type ii, jj, kk;
2152 size_type iw = in.width( );
2153 size_type ih = in.height( );
2154 size_type id = in.depth( );
2155 size_type ow = out.width( );
2156 size_type oh = out.height( );
2157 size_type od = out.depth( );
2158
2159 double sx = static_cast< double >( iw ) / static_cast< double >( ow );
2160 double sy = static_cast< double >( ih ) / static_cast< double >( oh );
2161 double sz = static_cast< double >( id ) / static_cast< double >( od );
2162 double x, y, z;
2163
2164 for( k = thread_idz ; k < od ; k += thread_numz )
2165 {
2166 z = sz * k;
2167 kk = static_cast< difference_type >( z );
2168 z -= kk;
2169
2170 for( j = thread_idy ; j < oh ; j += thread_numy )
2171 {
2172 y = sy * j;
2173 jj = static_cast< difference_type >( y );
2174 y -= jj;
2175
2176 for( i = thread_idx ; i < ow ; i += thread_numx )
2177 {
2178 x = sx * i;
2179 ii = static_cast< difference_type >( x );
2180 x -= ii;
2181
2182 __interpolate_utility__::round( _sinc_< is_color< value_type >::value >::interpolate( in, ii, jj, kk, x, y, z ), out( i, j, k ) );
2183 }
2184 }
2185 }
2186 }
2187 }
2188
2189
2190
2191
2192 // 画像補間のスレッド実装
2193 namespace __interpolate_controller__
2194 {
2195 // 最近傍型補間
2196 template < class T1, class Allocator1, class T2, class Allocator2 >
2197 void nearest__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
2198 typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num )
2199 {
2200 __nearest__::interpolate( in, out, thread_id, thread_num, 0, 1, 0, 1 );
2201 }
2202
2203 template < class T1, class Allocator1, class T2, class Allocator2 >
2204 void nearest__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2205 typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num )
2206 {
2207 __nearest__::interpolate( in, out, 0, 1, thread_id, thread_num, 0, 1 );
2208 }
2209
2210 template < class T1, class Allocator1, class T2, class Allocator2 >
2211 void nearest__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2212 typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num )
2213 {
2214 __nearest__::interpolate( in, out, 0, 1, 0, 1, thread_id, thread_num );
2215 }
2216
2217
2218 // 平均値型補間
2219 template < class T1, class Allocator1, class T2, class Allocator2 >
2220 void mean__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
2221 typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num )
2222 {
2223 __mean__::interpolate( in, out, thread_id, thread_num, 0, 1, 0, 1 );
2224 }
2225
2226 template < class T1, class Allocator1, class T2, class Allocator2 >
2227 void mean__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2228 typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num )
2229 {
2230 __mean__::interpolate( in, out, 0, 1, thread_id, thread_num, 0, 1 );
2231 }
2232
2233 template < class T1, class Allocator1, class T2, class Allocator2 >
2234 void mean__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2235 typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num )
2236 {
2237 __mean__::interpolate( in, out, 0, 1, 0, 1, thread_id, thread_num );
2238 }
2239
2240
2241 // 線形補間
2242 template < class T1, class Allocator1, class T2, class Allocator2 >
2243 void linear__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
2244 typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num )
2245 {
2246 __linear__::interpolate( in, out, thread_id, thread_num, 0, 1, 0, 1 );
2247 }
2248
2249 template < class T1, class Allocator1, class T2, class Allocator2 >
2250 void linear__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2251 typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num )
2252 {
2253 __linear__::interpolate( in, out, 0, 1, thread_id, thread_num, 0, 1 );
2254 }
2255
2256 template < class T1, class Allocator1, class T2, class Allocator2 >
2257 void linear__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2258 typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num )
2259 {
2260 __linear__::interpolate( in, out, 0, 1, 0, 1, thread_id, thread_num );
2261 }
2262
2263
2264 // 3次のsinc関数を用いた補間
2265 template < class T1, class Allocator1, class T2, class Allocator2 >
2266 void cubic__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
2267 typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num )
2268 {
2269 __cubic__::interpolate( in, out, thread_id, thread_num, 0, 1, 0, 1 );
2270 }
2271
2272 template < class T1, class Allocator1, class T2, class Allocator2 >
2273 void cubic__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2274 typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num )
2275 {
2276 __cubic__::interpolate( in, out, 0, 1, thread_id, thread_num, 0, 1 );
2277 }
2278
2279 template < class T1, class Allocator1, class T2, class Allocator2 >
2280 void cubic__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2281 typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num )
2282 {
2283 __cubic__::interpolate( in, out, 0, 1, 0, 1, thread_id, thread_num );
2284 }
2285
2286
2287 // 3次のbpline関数を用いた補間
2288 template < class T1, class Allocator1, class T2, class Allocator2 >
2289 void bspline__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
2290 typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num )
2291 {
2292 __bspline__::interpolate( in, out, thread_id, thread_num, 0, 1, 0, 1 );
2293 }
2294
2295 template < class T1, class Allocator1, class T2, class Allocator2 >
2296 void bspline__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2297 typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num )
2298 {
2299 __bspline__::interpolate( in, out, 0, 1, thread_id, thread_num, 0, 1 );
2300 }
2301
2302 template < class T1, class Allocator1, class T2, class Allocator2 >
2303 void bspline__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2304 typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num )
2305 {
2306 __bspline__::interpolate( in, out, 0, 1, 0, 1, thread_id, thread_num );
2307 }
2308
2309
2310 // フルスケールの sinc 関数を用いた補間
2311 template < class T1, class Allocator1, class T2, class Allocator2 >
2312 void sinc__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
2313 typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num )
2314 {
2315 __sinc__::interpolate( in, out, thread_id, thread_num, 0, 1, 0, 1 );
2316 }
2317
2318 template < class T1, class Allocator1, class T2, class Allocator2 >
2319 void sinc__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2320 typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num )
2321 {
2322 __sinc__::interpolate( in, out, 0, 1, thread_id, thread_num, 0, 1 );
2323 }
2324
2325 template < class T1, class Allocator1, class T2, class Allocator2 >
2326 void sinc__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2327 typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num )
2328 {
2329 __sinc__::interpolate( in, out, 0, 1, 0, 1, thread_id, thread_num );
2330 }
2331
2332
2333 template < class T1, class T2 >
2334 class interpolate_thread : public mist::thread< interpolate_thread< T1, T2 > >
2335 {
2336 public:
2337 typedef mist::thread< interpolate_thread< T1, T2 > > base;
2338 typedef typename base::thread_exit_type thread_exit_type;
2339 typedef typename T1::size_type size_type;
2340 typedef typename T1::value_type value_type;
2341
2342 enum Mode
2343 {
2344 Nearest,
2345 Mean,
2346 Linear,
2347 Cubic,
2348 Bspline,
2349 Sinc,
2350 };
2351
2352 private:
2353 size_t thread_id_;
2354 size_t thread_num_;
2355
2356 // 入出力用の画像へのポインタ
2357 const T1 *in_;
2358 T2 *out_;
2359 Mode interpolate_;
2360
2361 public:
2362 void setup_parameters( const T1 &in, T2 &out, Mode interpolate, size_type thread_id, size_type thread_num )
2363 {
2364 in_ = &in;
2365 out_ = &out;
2366 interpolate_ = interpolate;
2367 thread_id_ = thread_id;
2368 thread_num_ = thread_num;
2369 }
2370
2371 const interpolate_thread& operator =( const interpolate_thread &p )
2372 {
2373 if( &p != this )
2374 {
2375 base::operator =( p );
2376 thread_id_ = p.thread_id_;
2377 thread_num_ = p.thread_num_;
2378 in_ = p.in_;
2379 out_ = p.out_;
2380 interpolate_ = p.interpolate_;
2381 }
2382 return( *this );
2383 }
2384
2385 interpolate_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ),
2386 in_( NULL ), out_( NULL ), interpolate_( Linear )
2387 {
2388 }
2389 interpolate_thread( const interpolate_thread &p ) : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ),
2390 in_( p.in_ ), out_( p.out_ ), interpolate_( Linear )
2391 {
2392 }
2393
2394 protected:
2395 // 継承した先で必ず実装されるスレッド関数
2396 virtual thread_exit_type thread_function( )
2397 {
2398 switch( interpolate_ )
2399 {
2400 case Nearest:
2401 nearest__( *in_, *out_, thread_id_, thread_num_ );
2402 break;
2403
2404 case Mean:
2405 mean__( *in_, *out_, thread_id_, thread_num_ );
2406 break;
2407
2408 case Cubic:
2409 cubic__( *in_, *out_, thread_id_, thread_num_ );
2410 break;
2411
2412 case Bspline:
2413 bspline__( *in_, *out_, thread_id_, thread_num_ );
2414 break;
2415
2416 case Sinc:
2417 sinc__( *in_, *out_, thread_id_, thread_num_ );
2418 break;
2419
2420 case Linear:
2421 default:
2422 linear__( *in_, *out_, thread_id_, thread_num_ );
2423 break;
2424 }
2425
2426 return( true );
2427 }
2428 };
2429 }
2430
2431
2432
2440
2441
2443 namespace nearest
2444 {
2462 template < class T1, class Allocator1, class T2, class Allocator2 >
2463 bool interpolate( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
2464 typename array< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
2465 {
2466 if( is_same_object( in, out ) || width == 0 )
2467 {
2468 return( false );
2469 }
2470
2471 typedef typename array< T1, Allocator1 >::size_type size_type;
2472 typedef __interpolate_controller__::interpolate_thread< array< T1, Allocator1 >, array< T2, Allocator2 > > interpolate_thread;
2473
2474 if( thread_num == 0 )
2475 {
2476 thread_num = static_cast< size_type >( get_cpu_num( ) );
2477 }
2478
2479 size_type i;
2480 out.resize( width );
2481
2482 if( in.width( ) == width )
2483 {
2484 for( i = 0 ; i < in.size( ) ; i++ )
2485 {
2486 out[ i ] = static_cast< typename array< T2, Allocator2 >::value_type >( in[ i ] );
2487 }
2488 return( true );
2489 }
2490
2491 interpolate_thread *thread = new interpolate_thread[ thread_num ];
2492
2493 for( i = 0 ; i < thread_num ; i++ )
2494 {
2495 thread[ i ].setup_parameters( in, out, interpolate_thread::Nearest, i, thread_num );
2496 }
2497
2498 do_threads( thread, thread_num );
2499
2500 delete [] thread;
2501
2502 return( true );
2503 }
2504
2505
2523 template < class T1, class Allocator1, class T2, class Allocator2 >
2524 bool interpolate( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out,
2525 typename array1< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
2526 {
2527 if( is_same_object( in, out ) || width == 0 )
2528 {
2529 return( false );
2530 }
2531
2532 typedef typename array1< T1, Allocator1 >::size_type size_type;
2533 typedef __interpolate_controller__::interpolate_thread< array1< T1, Allocator1 >, array1< T2, Allocator2 > > interpolate_thread;
2534
2535 if( thread_num == 0 )
2536 {
2537 thread_num = static_cast< size_type >( get_cpu_num( ) );
2538 }
2539
2540 size_type i;
2541 out.resize( width );
2542 out.reso1( in.reso1( ) * static_cast< double >( in.size( ) ) / static_cast< double >( width ) );
2543
2544 if( in.width( ) == width )
2545 {
2546 for( i = 0 ; i < in.size( ) ; i++ )
2547 {
2548 out[ i ] = static_cast< typename array1< T2, Allocator2 >::value_type >( in[ i ] );
2549 }
2550 return( true );
2551 }
2552
2553 interpolate_thread *thread = new interpolate_thread[ thread_num ];
2554
2555 for( i = 0 ; i < thread_num ; i++ )
2556 {
2557 thread[ i ].setup_parameters( in, out, interpolate_thread::Nearest, i, thread_num );
2558 }
2559
2560 do_threads( thread, thread_num );
2561
2562 delete [] thread;
2563
2564 return( true );
2565 }
2566
2567
2586 template < class T1, class Allocator1, class T2, class Allocator2 >
2587 bool interpolate( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2588 typename array2< T1, Allocator1 >::size_type width, typename array2< T1, Allocator1 >::size_type height,
2589 typename array2< T1, Allocator1 >::size_type thread_num = 0 )
2590 {
2591 if( is_same_object( in, out ) || width == 0 || height == 0 )
2592 {
2593 return( false );
2594 }
2595
2596 typedef typename array2< T1, Allocator1 >::size_type size_type;
2597 typedef __interpolate_controller__::interpolate_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 > > interpolate_thread;
2598
2599 if( thread_num == 0 )
2600 {
2601 thread_num = static_cast< size_type >( get_cpu_num( ) );
2602 }
2603
2604 size_type i;
2605 out.resize( width, height );
2606 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
2607 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
2608
2609 if( in.width( ) == width && in.height( ) == height )
2610 {
2611 for( i = 0 ; i < in.size( ) ; i++ )
2612 {
2613 out[ i ] = static_cast< typename array2< T2, Allocator2 >::value_type >( in[ i ] );
2614 }
2615 return( true );
2616 }
2617
2618 interpolate_thread *thread = new interpolate_thread[ thread_num ];
2619
2620 for( i = 0 ; i < thread_num ; i++ )
2621 {
2622 thread[ i ].setup_parameters( in, out, interpolate_thread::Nearest, i, thread_num );
2623 }
2624
2625 do_threads( thread, thread_num );
2626
2627 delete [] thread;
2628
2629 return( true );
2630 }
2631
2632
2652 template < class T1, class Allocator1, class T2, class Allocator2 >
2653 bool interpolate( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2654 typename array3< T1, Allocator1 >::size_type width,
2655 typename array3< T1, Allocator1 >::size_type height,
2656 typename array3< T1, Allocator1 >::size_type depth,
2657 typename array3< T1, Allocator1 >::size_type thread_num = 0 )
2658 {
2659 if( is_same_object( in, out ) || width == 0 || height == 0 || depth == 0 )
2660 {
2661 return( false );
2662 }
2663
2664 typedef typename array3< T1, Allocator1 >::size_type size_type;
2665 typedef __interpolate_controller__::interpolate_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 > > interpolate_thread;
2666
2667 if( thread_num == 0 )
2668 {
2669 thread_num = static_cast< size_type >( get_cpu_num( ) );
2670 }
2671
2672 size_type i;
2673 out.resize( width, height, depth );
2674 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
2675 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
2676 out.reso3( in.reso3( ) * static_cast< double >( in.depth( ) ) / static_cast< double >( depth ) );
2677
2678 if( in.width( ) == width && in.height( ) == height && in.depth( ) == depth )
2679 {
2680 for( i = 0 ; i < in.size( ) ; i++ )
2681 {
2682 out[ i ] = static_cast< typename array3< T2, Allocator2 >::value_type >( in[ i ] );
2683 }
2684 return( true );
2685 }
2686
2687 interpolate_thread *thread = new interpolate_thread[ thread_num ];
2688
2689 for( i = 0 ; i < thread_num ; i++ )
2690 {
2691 thread[ i ].setup_parameters( in, out, interpolate_thread::Nearest, i, thread_num );
2692 }
2693
2694 do_threads( thread, thread_num );
2695
2696 delete [] thread;
2697
2698 return( true );
2699 }
2700 }
2701
2703 namespace mean
2704 {
2725 template < class T1, class Allocator1, class T2, class Allocator2 >
2726 bool interpolate( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
2727 typename array< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
2728 {
2729 if( is_same_object( in, out ) || width == 0 || in.size( ) < 2 )
2730 {
2731 return( false );
2732 }
2733 else if( in.width( ) < width )
2734 {
2735 // 平均値補間は画像拡大には使えない
2736 return( false );
2737 }
2738
2739 typedef typename array< T1, Allocator1 >::size_type size_type;
2740 typedef __interpolate_controller__::interpolate_thread< array< T1, Allocator1 >, array< T2, Allocator2 > > interpolate_thread;
2741
2742 if( thread_num == 0 )
2743 {
2744 thread_num = static_cast< size_type >( get_cpu_num( ) );
2745 }
2746
2747 size_type i;
2748 out.resize( width );
2749
2750 if( in.width( ) == width )
2751 {
2752 for( i = 0 ; i < in.size( ) ; i++ )
2753 {
2754 out[ i ] = static_cast< typename array< T2, Allocator2 >::value_type >( in[ i ] );
2755 }
2756 return( true );
2757 }
2758
2759 interpolate_thread *thread = new interpolate_thread[ thread_num ];
2760
2761 for( i = 0 ; i < thread_num ; i++ )
2762 {
2763 thread[ i ].setup_parameters( in, out, interpolate_thread::Mean, i, thread_num );
2764 }
2765
2766 do_threads( thread, thread_num );
2767
2768 delete [] thread;
2769
2770 return( true );
2771 }
2772
2773
2794 template < class T1, class Allocator1, class T2, class Allocator2 >
2795 bool interpolate( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out,
2796 typename array1< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
2797 {
2798 if( is_same_object( in, out ) || width == 0 || in.width( ) < 2 )
2799 {
2800 return( false );
2801 }
2802 else if( in.width( ) < width )
2803 {
2804 // 平均値補間は画像拡大には使えない
2805 return( false );
2806 }
2807
2808 typedef typename array1< T1, Allocator1 >::size_type size_type;
2809 typedef __interpolate_controller__::interpolate_thread< array1< T1, Allocator1 >, array1< T2, Allocator2 > > interpolate_thread;
2810
2811 if( thread_num == 0 )
2812 {
2813 thread_num = static_cast< size_type >( get_cpu_num( ) );
2814 }
2815
2816 size_type i;
2817 out.resize( width );
2818 out.reso1( in.reso1( ) * static_cast< double >( in.size( ) ) / static_cast< double >( width ) );
2819
2820 if( in.width( ) == width )
2821 {
2822 for( i = 0 ; i < in.size( ) ; i++ )
2823 {
2824 out[ i ] = static_cast< typename array1< T2, Allocator2 >::value_type >( in[ i ] );
2825 }
2826 return( true );
2827 }
2828
2829 interpolate_thread *thread = new interpolate_thread[ thread_num ];
2830
2831 for( i = 0 ; i < thread_num ; i++ )
2832 {
2833 thread[ i ].setup_parameters( in, out, interpolate_thread::Mean, i, thread_num );
2834 }
2835
2836 do_threads( thread, thread_num );
2837
2838 delete [] thread;
2839
2840 return( true );
2841 }
2842
2843
2865 template < class T1, class Allocator1, class T2, class Allocator2 >
2866 bool interpolate( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2867 typename array2< T1, Allocator1 >::size_type width, typename array2< T1, Allocator1 >::size_type height,
2868 typename array2< T1, Allocator1 >::size_type thread_num = 0 )
2869 {
2870 if( is_same_object( in, out ) || width == 0 || height == 0 || in.width( ) < 2 || in.height( ) < 2 )
2871 {
2872 return( false );
2873 }
2874 else if( in.width( ) < width || in.height( ) < height )
2875 {
2876 // 平均値補間は画像拡大には使えない
2877 return( false );
2878 }
2879
2880 typedef typename array2< T1, Allocator1 >::size_type size_type;
2881 typedef __interpolate_controller__::interpolate_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 > > interpolate_thread;
2882
2883 if( thread_num == 0 )
2884 {
2885 thread_num = static_cast< size_type >( get_cpu_num( ) );
2886 }
2887
2888 size_type i;
2889 out.resize( width, height );
2890 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
2891 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
2892
2893 if( in.width( ) == width && in.height( ) == height )
2894 {
2895 for( i = 0 ; i < in.size( ) ; i++ )
2896 {
2897 out[ i ] = static_cast< typename array2< T2, Allocator2 >::value_type >( in[ i ] );
2898 }
2899 return( true );
2900 }
2901
2902 interpolate_thread *thread = new interpolate_thread[ thread_num ];
2903
2904 for( i = 0 ; i < thread_num ; i++ )
2905 {
2906 thread[ i ].setup_parameters( in, out, interpolate_thread::Mean, i, thread_num );
2907 }
2908
2909 do_threads( thread, thread_num );
2910
2911 delete [] thread;
2912
2913 return( true );
2914 }
2915
2916
2939 template < class T1, class Allocator1, class T2, class Allocator2 >
2940 bool interpolate( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2941 typename array3< T1, Allocator1 >::size_type width,
2942 typename array3< T1, Allocator1 >::size_type height,
2943 typename array3< T1, Allocator1 >::size_type depth,
2944 typename array3< T1, Allocator1 >::size_type thread_num = 0 )
2945 {
2946 if( is_same_object( in, out ) || width == 0 || height == 0 || depth == 0 || in.width( ) < 2 || in.height( ) < 2 || in.depth( ) < 2 )
2947 {
2948 return( false );
2949 }
2950 else if( in.width( ) < width || in.height( ) < height || in.depth( ) < depth )
2951 {
2952 // 平均値補間は画像拡大には使えない
2953 return( false );
2954 }
2955
2956 typedef typename array3< T1, Allocator1 >::size_type size_type;
2957 typedef __interpolate_controller__::interpolate_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 > > interpolate_thread;
2958
2959 if( thread_num == 0 )
2960 {
2961 thread_num = static_cast< size_type >( get_cpu_num( ) );
2962 }
2963
2964 size_type i;
2965 out.resize( width, height, depth );
2966 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
2967 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
2968 out.reso3( in.reso3( ) * static_cast< double >( in.depth( ) ) / static_cast< double >( depth ) );
2969
2970 if( in.width( ) == width && in.height( ) == height && in.depth( ) == depth )
2971 {
2972 for( i = 0 ; i < in.size( ) ; i++ )
2973 {
2974 out[ i ] = static_cast< typename array3< T2, Allocator2 >::value_type >( in[ i ] );
2975 }
2976 return( true );
2977 }
2978
2979 interpolate_thread *thread = new interpolate_thread[ thread_num ];
2980
2981 for( i = 0 ; i < thread_num ; i++ )
2982 {
2983 thread[ i ].setup_parameters( in, out, interpolate_thread::Mean, i, thread_num );
2984 }
2985
2986 do_threads( thread, thread_num );
2987
2988 delete [] thread;
2989
2990 return( true );
2991 }
2992 }
2993
2994
2995
2997 namespace linear
2998 {
3016 template < class T1, class Allocator1, class T2, class Allocator2 >
3017 bool interpolate( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
3018 typename array< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
3019 {
3020 if( is_same_object( in, out ) || width == 0 )
3021 {
3022 return( false );
3023 }
3024
3025 typedef typename array< T1, Allocator1 >::size_type size_type;
3026 typedef __interpolate_controller__::interpolate_thread< array< T1, Allocator1 >, array< T2, Allocator2 > > interpolate_thread;
3027
3028 if( thread_num == 0 )
3029 {
3030 thread_num = static_cast< size_type >( get_cpu_num( ) );
3031 }
3032
3033 size_type i;
3034 out.resize( width );
3035
3036 if( in.width( ) == width )
3037 {
3038 for( i = 0 ; i < in.size( ) ; i++ )
3039 {
3040 out[ i ] = static_cast< typename array< T2, Allocator2 >::value_type >( in[ i ] );
3041 }
3042 return( true );
3043 }
3044
3045 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3046
3047 for( i = 0 ; i < thread_num ; i++ )
3048 {
3049 thread[ i ].setup_parameters( in, out, interpolate_thread::Linear, i, thread_num );
3050 }
3051
3052 do_threads( thread, thread_num );
3053
3054 delete [] thread;
3055
3056 return( true );
3057 }
3058
3059
3077 template < class T1, class Allocator1, class T2, class Allocator2 >
3078 bool interpolate( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out,
3079 typename array1< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
3080 {
3081 if( is_same_object( in, out ) || width == 0 )
3082 {
3083 return( false );
3084 }
3085
3086 typedef typename array1< T1, Allocator1 >::size_type size_type;
3087 typedef __interpolate_controller__::interpolate_thread< array1< T1, Allocator1 >, array1< T2, Allocator2 > > interpolate_thread;
3088
3089 if( thread_num == 0 )
3090 {
3091 thread_num = static_cast< size_type >( get_cpu_num( ) );
3092 }
3093
3094 size_type i;
3095 out.resize( width );
3096 out.reso1( in.reso1( ) * static_cast< double >( in.size( ) ) / static_cast< double >( width ) );
3097
3098 if( in.width( ) == width )
3099 {
3100 for( i = 0 ; i < in.size( ) ; i++ )
3101 {
3102 out[ i ] = static_cast< typename array1< T2, Allocator2 >::value_type >( in[ i ] );
3103 }
3104 return( true );
3105 }
3106
3107 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3108
3109 for( i = 0 ; i < thread_num ; i++ )
3110 {
3111 thread[ i ].setup_parameters( in, out, interpolate_thread::Linear, i, thread_num );
3112 }
3113
3114 do_threads( thread, thread_num );
3115
3116 delete [] thread;
3117
3118 return( true );
3119 }
3120
3121
3140 template < class T1, class Allocator1, class T2, class Allocator2 >
3141 bool interpolate( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
3142 typename array2< T1, Allocator1 >::size_type width, typename array2< T1, Allocator1 >::size_type height,
3143 typename array2< T1, Allocator1 >::size_type thread_num = 0 )
3144 {
3145 if( is_same_object( in, out ) || width == 0 || height == 0 )
3146 {
3147 return( false );
3148 }
3149
3150 typedef typename array2< T1, Allocator1 >::size_type size_type;
3151 typedef __interpolate_controller__::interpolate_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 > > interpolate_thread;
3152
3153 if( thread_num == 0 )
3154 {
3155 thread_num = static_cast< size_type >( get_cpu_num( ) );
3156 }
3157
3158 size_type i;
3159 out.resize( width, height );
3160 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
3161 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
3162
3163 if( in.width( ) == width && in.height( ) == height )
3164 {
3165 for( i = 0 ; i < in.size( ) ; i++ )
3166 {
3167 out[ i ] = static_cast< typename array2< T2, Allocator2 >::value_type >( in[ i ] );
3168 }
3169 return( true );
3170 }
3171
3172 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3173
3174 for( i = 0 ; i < thread_num ; i++ )
3175 {
3176 thread[ i ].setup_parameters( in, out, interpolate_thread::Linear, i, thread_num );
3177 }
3178
3179 do_threads( thread, thread_num );
3180
3181 delete [] thread;
3182
3183 return( true );
3184 }
3185
3186
3206 template < class T1, class Allocator1, class T2, class Allocator2 >
3207 bool interpolate( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
3208 typename array3< T1, Allocator1 >::size_type width,
3209 typename array3< T1, Allocator1 >::size_type height,
3210 typename array3< T1, Allocator1 >::size_type depth,
3211 typename array3< T1, Allocator1 >::size_type thread_num = 0 )
3212 {
3213 if( is_same_object( in, out ) || width == 0 || height == 0 || depth == 0 )
3214 {
3215 return( false );
3216 }
3217
3218 typedef typename array3< T1, Allocator1 >::size_type size_type;
3219 typedef __interpolate_controller__::interpolate_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 > > interpolate_thread;
3220
3221 if( thread_num == 0 )
3222 {
3223 thread_num = static_cast< size_type >( get_cpu_num( ) );
3224 }
3225
3226 size_type i;
3227 out.resize( width, height, depth );
3228 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
3229 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
3230 out.reso3( in.reso3( ) * static_cast< double >( in.depth( ) ) / static_cast< double >( depth ) );
3231
3232 if( in.width( ) == width && in.height( ) == height && in.depth( ) == depth )
3233 {
3234 for( i = 0 ; i < in.size( ) ; i++ )
3235 {
3236 out[ i ] = static_cast< typename array3< T2, Allocator2 >::value_type >( in[ i ] );
3237 }
3238 return( true );
3239 }
3240
3241 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3242
3243 for( i = 0 ; i < thread_num ; i++ )
3244 {
3245 thread[ i ].setup_parameters( in, out, interpolate_thread::Linear, i, thread_num );
3246 }
3247
3248 do_threads( thread, thread_num );
3249
3250 delete [] thread;
3251
3252 return( true );
3253 }
3254 }
3255
3256
3257
3259 namespace cubic
3260 {
3278 template < class T1, class Allocator1, class T2, class Allocator2 >
3279 bool interpolate( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
3280 typename array< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
3281 {
3282 if( is_same_object( in, out ) || width == 0 )
3283 {
3284 return( false );
3285 }
3286
3287 typedef typename array< T1, Allocator1 >::size_type size_type;
3288 typedef __interpolate_controller__::interpolate_thread< array< T1, Allocator1 >, array< T2, Allocator2 > > interpolate_thread;
3289
3290 if( thread_num == 0 )
3291 {
3292 thread_num = static_cast< size_type >( get_cpu_num( ) );
3293 }
3294
3295 size_type i;
3296 out.resize( width );
3297
3298 if( in.width( ) == width )
3299 {
3300 for( i = 0 ; i < in.size( ) ; i++ )
3301 {
3302 out[ i ] = static_cast< typename array< T2, Allocator2 >::value_type >( in[ i ] );
3303 }
3304 return( true );
3305 }
3306
3307 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3308
3309 for( i = 0 ; i < thread_num ; i++ )
3310 {
3311 thread[ i ].setup_parameters( in, out, interpolate_thread::Cubic, i, thread_num );
3312 }
3313
3314 do_threads( thread, thread_num );
3315
3316 delete [] thread;
3317
3318 return( true );
3319 }
3320
3338 template < class T1, class Allocator1, class T2, class Allocator2 >
3339 bool interpolate( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out,
3340 typename array1< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
3341 {
3342 if( is_same_object( in, out ) || width == 0 )
3343 {
3344 return( false );
3345 }
3346
3347 typedef typename array1< T1, Allocator1 >::size_type size_type;
3348 typedef __interpolate_controller__::interpolate_thread< array1< T1, Allocator1 >, array1< T2, Allocator2 > > interpolate_thread;
3349
3350 if( thread_num == 0 )
3351 {
3352 thread_num = static_cast< size_type >( get_cpu_num( ) );
3353 }
3354
3355 size_type i;
3356 out.resize( width );
3357 out.reso1( in.reso1( ) * static_cast< double >( in.size( ) ) / static_cast< double >( width ) );
3358
3359 if( in.width( ) == width )
3360 {
3361 for( i = 0 ; i < in.size( ) ; i++ )
3362 {
3363 out[ i ] = static_cast< typename array1< T2, Allocator2 >::value_type >( in[ i ] );
3364 }
3365 return( true );
3366 }
3367
3368 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3369
3370 for( i = 0 ; i < thread_num ; i++ )
3371 {
3372 thread[ i ].setup_parameters( in, out, interpolate_thread::Cubic, i, thread_num );
3373 }
3374
3375 do_threads( thread, thread_num );
3376
3377 delete [] thread;
3378
3379 return( true );
3380 }
3381
3382
3401 template < class T1, class Allocator1, class T2, class Allocator2 >
3402 bool interpolate( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
3403 typename array2< T1, Allocator1 >::size_type width, typename array2< T1, Allocator1 >::size_type height,
3404 typename array2< T1, Allocator1 >::size_type thread_num = 0 )
3405 {
3406 if( is_same_object( in, out ) || width == 0 || height == 0 )
3407 {
3408 return( false );
3409 }
3410
3411 typedef typename array2< T1, Allocator1 >::size_type size_type;
3412 typedef __interpolate_controller__::interpolate_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 > > interpolate_thread;
3413
3414 if( thread_num == 0 )
3415 {
3416 thread_num = static_cast< size_type >( get_cpu_num( ) );
3417 }
3418
3419 size_type i;
3420 out.resize( width, height );
3421 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
3422 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
3423
3424 if( in.width( ) == width && in.height( ) == height )
3425 {
3426 for( i = 0 ; i < in.size( ) ; i++ )
3427 {
3428 out[ i ] = static_cast< typename array2< T2, Allocator2 >::value_type >( in[ i ] );
3429 }
3430 return( true );
3431 }
3432
3433 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3434
3435 for( i = 0 ; i < thread_num ; i++ )
3436 {
3437 thread[ i ].setup_parameters( in, out, interpolate_thread::Cubic, i, thread_num );
3438 }
3439
3440 do_threads( thread, thread_num );
3441
3442 delete [] thread;
3443
3444 return( true );
3445 }
3446
3447
3467 template < class T1, class Allocator1, class T2, class Allocator2 >
3468 bool interpolate( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
3469 typename array3< T1, Allocator1 >::size_type width,
3470 typename array3< T1, Allocator1 >::size_type height,
3471 typename array3< T1, Allocator1 >::size_type depth,
3472 typename array3< T1, Allocator1 >::size_type thread_num = 0 )
3473 {
3474 if( is_same_object( in, out ) || width == 0 || height == 0 || depth == 0 )
3475 {
3476 return( false );
3477 }
3478
3479 typedef typename array3< T1, Allocator1 >::size_type size_type;
3480 typedef __interpolate_controller__::interpolate_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 > > interpolate_thread;
3481
3482 if( thread_num == 0 )
3483 {
3484 thread_num = static_cast< size_type >( get_cpu_num( ) );
3485 }
3486
3487 size_type i;
3488 out.resize( width, height, depth );
3489 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
3490 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
3491 out.reso3( in.reso3( ) * static_cast< double >( in.depth( ) ) / static_cast< double >( depth ) );
3492
3493 if( in.width( ) == width && in.height( ) == height && in.depth( ) == depth )
3494 {
3495 for( i = 0 ; i < in.size( ) ; i++ )
3496 {
3497 out[ i ] = static_cast< typename array3< T2, Allocator2 >::value_type >( in[ i ] );
3498 }
3499 return( true );
3500 }
3501
3502 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3503
3504 for( i = 0 ; i < thread_num ; i++ )
3505 {
3506 thread[ i ].setup_parameters( in, out, interpolate_thread::Cubic, i, thread_num );
3507 }
3508
3509 do_threads( thread, thread_num );
3510
3511 delete [] thread;
3512
3513 return( true );
3514 }
3515 }
3516
3517
3518
3520 namespace BSpline
3521 {
3539 template < class T1, class Allocator1, class T2, class Allocator2 >
3540 bool interpolate( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
3541 typename array< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
3542 {
3543 if( is_same_object( in, out ) || width == 0 )
3544 {
3545 return( false );
3546 }
3547
3548 typedef typename array< T1, Allocator1 >::size_type size_type;
3549 typedef __interpolate_controller__::interpolate_thread< array< T1, Allocator1 >, array< T2, Allocator2 > > interpolate_thread;
3550
3551 if( thread_num == 0 )
3552 {
3553 thread_num = static_cast< size_type >( get_cpu_num( ) );
3554 }
3555
3556 size_type i;
3557 out.resize( width );
3558
3559 if( in.width( ) == width )
3560 {
3561 for( i = 0 ; i < in.size( ) ; i++ )
3562 {
3563 out[ i ] = static_cast< typename array< T2, Allocator2 >::value_type >( in[ i ] );
3564 }
3565 return( true );
3566 }
3567
3568 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3569
3570 for( i = 0 ; i < thread_num ; i++ )
3571 {
3572 thread[ i ].setup_parameters( in, out, interpolate_thread::Bspline, i, thread_num );
3573 }
3574
3575 do_threads( thread, thread_num );
3576
3577 delete [] thread;
3578
3579 return( true );
3580 }
3581
3599 template < class T1, class Allocator1, class T2, class Allocator2 >
3600 bool interpolate( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out,
3601 typename array1< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
3602 {
3603 if( is_same_object( in, out ) || width == 0 )
3604 {
3605 return( false );
3606 }
3607
3608 typedef typename array1< T1, Allocator1 >::size_type size_type;
3609 typedef __interpolate_controller__::interpolate_thread< array1< T1, Allocator1 >, array1< T2, Allocator2 > > interpolate_thread;
3610
3611 if( thread_num == 0 )
3612 {
3613 thread_num = static_cast< size_type >( get_cpu_num( ) );
3614 }
3615
3616 size_type i;
3617 out.resize( width );
3618 out.reso1( in.reso1( ) * static_cast< double >( in.size( ) ) / static_cast< double >( width ) );
3619
3620 if( in.width( ) == width )
3621 {
3622 for( i = 0 ; i < in.size( ) ; i++ )
3623 {
3624 out[ i ] = static_cast< typename array1< T2, Allocator2 >::value_type >( in[ i ] );
3625 }
3626 return( true );
3627 }
3628
3629 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3630
3631 for( i = 0 ; i < thread_num ; i++ )
3632 {
3633 thread[ i ].setup_parameters( in, out, interpolate_thread::Bspline, i, thread_num );
3634 }
3635
3636 do_threads( thread, thread_num );
3637
3638 delete [] thread;
3639
3640 return( true );
3641 }
3642
3643
3662 template < class T1, class Allocator1, class T2, class Allocator2 >
3663 bool interpolate( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
3664 typename array2< T1, Allocator1 >::size_type width, typename array2< T1, Allocator1 >::size_type height,
3665 typename array2< T1, Allocator1 >::size_type thread_num = 0 )
3666 {
3667 if( is_same_object( in, out ) || width == 0 || height == 0 )
3668 {
3669 return( false );
3670 }
3671
3672 typedef typename array2< T1, Allocator1 >::size_type size_type;
3673 typedef __interpolate_controller__::interpolate_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 > > interpolate_thread;
3674
3675 if( thread_num == 0 )
3676 {
3677 thread_num = static_cast< size_type >( get_cpu_num( ) );
3678 }
3679
3680 size_type i;
3681 out.resize( width, height );
3682 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
3683 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
3684
3685 if( in.width( ) == width && in.height( ) == height )
3686 {
3687 for( i = 0 ; i < in.size( ) ; i++ )
3688 {
3689 out[ i ] = static_cast< typename array2< T2, Allocator2 >::value_type >( in[ i ] );
3690 }
3691 return( true );
3692 }
3693
3694 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3695
3696 for( i = 0 ; i < thread_num ; i++ )
3697 {
3698 thread[ i ].setup_parameters( in, out, interpolate_thread::Bspline, i, thread_num );
3699 }
3700
3701 do_threads( thread, thread_num );
3702
3703 delete [] thread;
3704
3705 return( true );
3706 }
3707
3708
3728 template < class T1, class Allocator1, class T2, class Allocator2 >
3729 bool interpolate( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
3730 typename array3< T1, Allocator1 >::size_type width,
3731 typename array3< T1, Allocator1 >::size_type height,
3732 typename array3< T1, Allocator1 >::size_type depth,
3733 typename array3< T1, Allocator1 >::size_type thread_num = 0 )
3734 {
3735 if( is_same_object( in, out ) || width == 0 || height == 0 || depth == 0 )
3736 {
3737 return( false );
3738 }
3739
3740 typedef typename array3< T1, Allocator1 >::size_type size_type;
3741 typedef __interpolate_controller__::interpolate_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 > > interpolate_thread;
3742
3743 if( thread_num == 0 )
3744 {
3745 thread_num = static_cast< size_type >( get_cpu_num( ) );
3746 }
3747
3748 size_type i;
3749 out.resize( width, height, depth );
3750 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
3751 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
3752 out.reso3( in.reso3( ) * static_cast< double >( in.depth( ) ) / static_cast< double >( depth ) );
3753
3754 if( in.width( ) == width && in.height( ) == height && in.depth( ) == depth )
3755 {
3756 for( i = 0 ; i < in.size( ) ; i++ )
3757 {
3758 out[ i ] = static_cast< typename array3< T2, Allocator2 >::value_type >( in[ i ] );
3759 }
3760 return( true );
3761 }
3762
3763 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3764
3765 for( i = 0 ; i < thread_num ; i++ )
3766 {
3767 thread[ i ].setup_parameters( in, out, interpolate_thread::Bspline, i, thread_num );
3768 }
3769
3770 do_threads( thread, thread_num );
3771
3772 delete [] thread;
3773
3774 return( true );
3775 }
3776 }
3777
3778
3779
3780
3782 namespace sinc
3783 {
3801 template < class T1, class Allocator1, class T2, class Allocator2 >
3802 bool interpolate( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
3803 typename array< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
3804 {
3805 if( is_same_object( in, out ) || width == 0 )
3806 {
3807 return( false );
3808 }
3809
3810 typedef typename array< T1, Allocator1 >::size_type size_type;
3811 typedef __interpolate_controller__::interpolate_thread< array< T1, Allocator1 >, array< T2, Allocator2 > > interpolate_thread;
3812
3813 if( thread_num == 0 )
3814 {
3815 thread_num = static_cast< size_type >( get_cpu_num( ) );
3816 }
3817
3818 size_type i;
3819 out.resize( width );
3820
3821 if( in.width( ) == width )
3822 {
3823 for( i = 0 ; i < in.size( ) ; i++ )
3824 {
3825 out[ i ] = static_cast< typename array< T2, Allocator2 >::value_type >( in[ i ] );
3826 }
3827 return( true );
3828 }
3829
3830 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3831
3832 for( i = 0 ; i < thread_num ; i++ )
3833 {
3834 thread[ i ].setup_parameters( in, out, interpolate_thread::Sinc, i, thread_num );
3835 }
3836
3837 do_threads( thread, thread_num );
3838
3839 delete [] thread;
3840
3841 return( true );
3842 }
3843
3861 template < class T1, class Allocator1, class T2, class Allocator2 >
3862 bool interpolate( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out,
3863 typename array1< T1, Allocator1 >::size_type width, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
3864 {
3865 if( is_same_object( in, out ) || width == 0 )
3866 {
3867 return( false );
3868 }
3869
3870 typedef typename array1< T1, Allocator1 >::size_type size_type;
3871 typedef __interpolate_controller__::interpolate_thread< array1< T1, Allocator1 >, array1< T2, Allocator2 > > interpolate_thread;
3872
3873 if( thread_num == 0 )
3874 {
3875 thread_num = static_cast< size_type >( get_cpu_num( ) );
3876 }
3877
3878 size_type i;
3879 out.resize( width );
3880 out.reso1( in.reso1( ) * static_cast< double >( in.size( ) ) / static_cast< double >( width ) );
3881
3882 if( in.width( ) == width )
3883 {
3884 for( i = 0 ; i < in.size( ) ; i++ )
3885 {
3886 out[ i ] = static_cast< typename array1< T2, Allocator2 >::value_type >( in[ i ] );
3887 }
3888 return( true );
3889 }
3890
3891 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3892
3893 for( i = 0 ; i < thread_num ; i++ )
3894 {
3895 thread[ i ].setup_parameters( in, out, interpolate_thread::Sinc, i, thread_num );
3896 }
3897
3898 do_threads( thread, thread_num );
3899
3900 delete [] thread;
3901
3902 return( true );
3903 }
3904
3905
3924 template < class T1, class Allocator1, class T2, class Allocator2 >
3925 bool interpolate( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
3926 typename array2< T1, Allocator1 >::size_type width, typename array2< T1, Allocator1 >::size_type height,
3927 typename array2< T1, Allocator1 >::size_type thread_num = 0 )
3928 {
3929 if( is_same_object( in, out ) || width == 0 || height == 0 )
3930 {
3931 return( false );
3932 }
3933
3934 typedef typename array2< T1, Allocator1 >::size_type size_type;
3935 typedef __interpolate_controller__::interpolate_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 > > interpolate_thread;
3936
3937 if( thread_num == 0 )
3938 {
3939 thread_num = static_cast< size_type >( get_cpu_num( ) );
3940 }
3941
3942 size_type i;
3943 out.resize( width, height );
3944 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
3945 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
3946
3947 if( in.width( ) == width && in.height( ) == height )
3948 {
3949 for( i = 0 ; i < in.size( ) ; i++ )
3950 {
3951 out[ i ] = static_cast< typename array2< T2, Allocator2 >::value_type >( in[ i ] );
3952 }
3953 return( true );
3954 }
3955
3956 interpolate_thread *thread = new interpolate_thread[ thread_num ];
3957
3958 for( i = 0 ; i < thread_num ; i++ )
3959 {
3960 thread[ i ].setup_parameters( in, out, interpolate_thread::Sinc, i, thread_num );
3961 }
3962
3963 do_threads( thread, thread_num );
3964
3965 delete [] thread;
3966
3967 return( true );
3968 }
3969
3970
3990 template < class T1, class Allocator1, class T2, class Allocator2 >
3991 bool interpolate( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
3992 typename array3< T1, Allocator1 >::size_type width,
3993 typename array3< T1, Allocator1 >::size_type height,
3994 typename array3< T1, Allocator1 >::size_type depth,
3995 typename array3< T1, Allocator1 >::size_type thread_num = 0 )
3996 {
3997 if( is_same_object( in, out ) || width == 0 || height == 0 || depth == 0 )
3998 {
3999 return( false );
4000 }
4001
4002 typedef typename array3< T1, Allocator1 >::size_type size_type;
4003 typedef __interpolate_controller__::interpolate_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 > > interpolate_thread;
4004
4005 if( thread_num == 0 )
4006 {
4007 thread_num = static_cast< size_type >( get_cpu_num( ) );
4008 }
4009
4010 size_type i;
4011 out.resize( width, height, depth );
4012 out.reso1( in.reso1( ) * static_cast< double >( in.width( ) ) / static_cast< double >( width ) );
4013 out.reso2( in.reso2( ) * static_cast< double >( in.height( ) ) / static_cast< double >( height ) );
4014 out.reso3( in.reso3( ) * static_cast< double >( in.depth( ) ) / static_cast< double >( depth ) );
4015
4016 if( in.width( ) == width && in.height( ) == height && in.depth( ) == depth )
4017 {
4018 for( i = 0 ; i < in.size( ) ; i++ )
4019 {
4020 out[ i ] = static_cast< typename array3< T2, Allocator2 >::value_type >( in[ i ] );
4021 }
4022 return( true );
4023 }
4024
4025 interpolate_thread *thread = new interpolate_thread[ thread_num ];
4026
4027 for( i = 0 ; i < thread_num ; i++ )
4028 {
4029 thread[ i ].setup_parameters( in, out, interpolate_thread::Sinc, i, thread_num );
4030 }
4031
4032 do_threads( thread, thread_num );
4033
4034 delete [] thread;
4035
4036 return( true );
4037 }
4038 }
4039
4041 // 画像補間グループの終わり
4042
4043
4044 // mist名前空間の終わり
4045 _MIST_END
4046
4047
4048 #endif // __INCLUDE_MIST_INTERPOLATE__
4049