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 
Generated on Wed Nov 12 2014 19:44:16 for MIST by doxygen 1.8.1.2