linear.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 
34 
35 
36 #ifndef __INCLUDE_FILTER_LINEAR_FILTER_H__
37 #define __INCLUDE_FILTER_LINEAR_FILTER_H__
38 
39 
40 
41 #ifndef __INCLUDE_MIST_H__
42 #include "../mist.h"
43 #endif
44 
45 #ifndef __INCLUDE_MIST_COLOR_H__
46 #include "../config/color.h"
47 #endif
48 
49 #ifndef __INCLUDE_MIST_VECTOR__
50 #include "../vector.h"
51 #endif
52 
53 #ifndef __INCLUDE_MIST_TYPE_TRAIT_H__
54 #include "../config/type_trait.h"
55 #endif
56 
57 #ifndef __INCLUDE_MIST_THREAD__
58 #include "../thread.h"
59 #endif
60 
61 
62 #include <cmath>
63 
64 
65 // mist名前空間の始まり
66 _MIST_BEGIN
67 
68 
69 // 線形フィルタ
70 namespace __linear__
71 {
72  template < class T >
73  struct __promote_pixel_converter_
74  {
75  typedef T value_type;
76  typedef double promote_type;
77 
78  static promote_type convert_to( const value_type &pixel )
79  {
80  return( pixel );
81  }
82 
83  static value_type convert_from( const promote_type &pixel )
84  {
85  return( static_cast< value_type >( pixel ) );
86  }
87  };
88 
89  template < class T >
90  struct __promote_pixel_converter_< rgb< T > >
91  {
92  typedef rgb< T > value_type;
93  typedef rgb< double > promote_type;
94 
95  static promote_type convert_to( const value_type &pixel )
96  {
97  return( promote_type( pixel.r, pixel.g, pixel.b ) );
98  }
99 
100  static value_type convert_from( const promote_type &pixel )
101  {
102  return( value_type( static_cast< T >( pixel.r ), static_cast< T >( pixel.g ), static_cast< T >( pixel.b ) ) );
103  }
104  };
105 
106  template < class T >
107  struct __promote_pixel_converter_< bgr< T > >
108  {
109  typedef bgr< T > value_type;
110  typedef bgr< double > promote_type;
111 
112  static promote_type convert_to( const value_type &pixel )
113  {
114  return( promote_type( pixel.r, pixel.g, pixel.b ) );
115  }
116 
117  static value_type convert_from( const promote_type &pixel )
118  {
119  return( value_type( static_cast< T >( pixel.r ), static_cast< T >( pixel.g ), static_cast< T >( pixel.b ) ) );
120  }
121  };
122 
123  template < class T >
124  struct __promote_pixel_converter_< rgba< T > >
125  {
126  typedef rgba< T > value_type;
127  typedef rgba< double > promote_type;
128 
129  static promote_type convert_to( const value_type &pixel )
130  {
131  return( promote_type( pixel.r, pixel.g, pixel.b ) );
132  }
133 
134  static value_type convert_from( const promote_type &pixel )
135  {
136  return( value_type( static_cast< T >( pixel.r ), static_cast< T >( pixel.g ), static_cast< T >( pixel.b ) ) );
137  }
138  };
139 
140  template < class T >
141  struct __promote_pixel_converter_< bgra< T > >
142  {
143  typedef bgra< T > value_type;
144  typedef bgra< double > promote_type;
145 
146  static promote_type convert_to( const value_type &pixel )
147  {
148  return( promote_type( pixel.r, pixel.g, pixel.b ) );
149  }
150 
151  static value_type convert_from( const promote_type &pixel )
152  {
153  return( value_type( static_cast< T >( pixel.r ), static_cast< T >( pixel.g ), static_cast< T >( pixel.b ) ) );
154  }
155  };
156 
157  template < class T >
158  struct __promote_pixel_converter_< vector2< T > >
159  {
160  typedef vector2< T > value_type;
161  typedef vector2< double > promote_type;
162 
163  static promote_type convert_to( const value_type &pixel )
164  {
165  return( promote_type( pixel.x, pixel.y ) );
166  }
167 
168  static value_type convert_from( const promote_type &pixel )
169  {
170  return( value_type( static_cast< T >( pixel.x ), static_cast< T >( pixel.y ) ) );
171  }
172  };
173 
174  template < class T >
175  struct __promote_pixel_converter_< vector3< T > >
176  {
177  typedef vector3< T > value_type;
178  typedef vector3< double > promote_type;
179 
180  static promote_type convert_to( const value_type &pixel )
181  {
182  return( promote_type( pixel.x, pixel.y, pixel.z ) );
183  }
184 
185  static value_type convert_from( const promote_type &pixel )
186  {
187  return( value_type( static_cast< T >( pixel.x ), static_cast< T >( pixel.y ), static_cast< T >( pixel.z ) ) );
188  }
189  };
190 
191  template < int DIMENSION >
192  struct __access__
193  {
194  template < class Array >
195  inline static typename Array::value_type &at( Array &in, typename Array::size_type _1, typename Array::size_type _2, typename Array::size_type _3 )
196  {
197  return( in( _1, _2, _3 ) );
198  }
199 
200  template < class Array >
201  inline static const typename Array::value_type &at( const Array &in, typename Array::size_type _1, typename Array::size_type _2, typename Array::size_type _3 )
202  {
203  return( in( _1, _2, _3 ) );
204  }
205 
206  template < class Array >
207  inline static typename Array::size_type size1( const Array &in ){ return( in.size1( ) ); }
208 
209  template < class Array >
210  inline static typename Array::size_type size2( const Array &in ){ return( in.size2( ) ); }
211 
212  template < class Array >
213  inline static typename Array::size_type size3( const Array &in ){ return( in.size3( ) ); }
214  };
215 
216  template < >
217  struct __access__< 2 >
218  {
219  template < class Array >
220  inline static typename Array::value_type &at( Array &in, typename Array::size_type _1, typename Array::size_type _2, typename Array::size_type _3 )
221  {
222  return( in( _2, _1, _3 ) );
223  }
224 
225  template < class Array >
226  inline static const typename Array::value_type &at( const Array &in, typename Array::size_type _1, typename Array::size_type _2, typename Array::size_type _3 )
227  {
228  return( in( _2, _1, _3 ) );
229  }
230 
231  template < class Array >
232  inline static typename Array::size_type size1( const Array &in ){ return( in.size2( ) ); }
233 
234  template < class Array >
235  inline static typename Array::size_type size2( const Array &in ){ return( in.size1( ) ); }
236 
237  template < class Array >
238  inline static typename Array::size_type size3( const Array &in ){ return( in.size3( ) ); }
239  };
240 
241  template < >
242  struct __access__< 3 >
243  {
244  template < class Array >
245  inline static typename Array::value_type &at( Array &in, typename Array::size_type _1, typename Array::size_type _2, typename Array::size_type _3 )
246  {
247  return( in( _2, _3, _1 ) );
248  }
249 
250  template < class Array >
251  inline static const typename Array::value_type &at( const Array &in, typename Array::size_type _1, typename Array::size_type _2, typename Array::size_type _3 )
252  {
253  return( in( _2, _3, _1 ) );
254  }
255 
256  template < class Array >
257  inline static typename Array::size_type size1( const Array &in ){ return( in.size3( ) ); }
258 
259  template < class Array >
260  inline static typename Array::size_type size2( const Array &in ){ return( in.size1( ) ); }
261 
262  template < class Array >
263  inline static typename Array::size_type size3( const Array &in ){ return( in.size2( ) ); }
264  };
265 
266  // in : 入力画像. 入力画像の画素値は min と max の間とする
267  // out : 出力画像. 出力画像のメモリはあらかじめ割り当てられているものとする
268  // fw, fh, fd : フィルタサイズ
269  template < class Array1, class Array2, class Kernel, class Functor >
270  void linear_filter( const Array1 &in, Array2 &out, const Kernel &kernel,
271  typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
272  typename Array1::size_type thread_idz, typename Array1::size_type thread_numz, Functor f )
273  {
274  typedef typename Array1::size_type size_type;
275  typedef typename Array1::difference_type difference_type;
276  typedef __promote_pixel_converter_< typename Array2::value_type > promote_pixel_converter;
277  typedef typename promote_pixel_converter::promote_type promote_type;
278  typedef typename Array2::value_type out_value_type;
279  typedef typename Array1::const_pointer ipointer_type;
280  typedef typename Array2::pointer opointer_type;
281  typedef typename Kernel::const_pointer kpointer_type;
282  typedef typename Kernel::value_type kvalue_type;
283 
284  difference_type w = in.width( );
285  difference_type h = in.height( );
286  difference_type d = in.depth( );
287 
288  difference_type fw = kernel.width( );
289  difference_type fh = kernel.height( );
290  difference_type fd = kernel.depth( );
291 
292  const bool bprogress1 = thread_idy == 0 && d == 1;
293  const bool bprogress2 = thread_idz == 0 && d > 1;
294 
295  difference_type rw = fw / 2;
296  difference_type rh = fh / 2;
297  difference_type rd = fd / 2;
298 
299  difference_type i, j, k;
300 
301  size_type fsize = fw * fh * fd;
302  difference_type *pf = new difference_type[ fsize ];
303 
304  {
305  difference_type count = 0;
306  difference_type cx = w / 2;
307  difference_type cy = h / 2;
308  difference_type cz = d / 2;
309  ipointer_type pc = &in( cx, cy, cz );
310 
311  for( difference_type z = 0 ; z < fd ; z++ )
312  {
313  for( difference_type y = 0 ; y < fh ; y++ )
314  {
315  for( difference_type x = 0 ; x < fw ; x++ )
316  {
317  pf[ count++ ] = &in( x + cx - rw, y + cy - rh, z + cz - rd ) - pc;
318  }
319  }
320  }
321  }
322 
323  // 画像の縁の処理を行う
324  for( k = thread_idz ; k < rd ; k += thread_numz )
325  {
326  for( j = thread_idy ; j < h ; j += thread_numy )
327  {
328  opointer_type op = &out( 0, j, k );
329 
330  for( i = 0 ; i < w ; i++ )
331  {
332  promote_type value = promote_type( );
333  kvalue_type sum = 0;
334  difference_type sl = i < rw ? rw - i : 0;
335  difference_type el = w - i <= rw ? rw + w - i : fw;
336  difference_type sm = j < rh ? rh - j : 0;
337  difference_type em = h - j <= rh ? rh + h - j : fh;
338  difference_type sn = k < rd ? rd - k : 0;
339 
340  for( difference_type n = sn ; n < fd ; n++ )
341  {
342  for( difference_type m = sm ; m < em ; m++ )
343  {
344  for( difference_type l = sl ; l < el ; l++ )
345  {
346  kvalue_type kv = kernel( l, m, n );
347  value += kv * in( i + l - rw, j + m - rh, k + n - rd );
348  sum += kv;
349  }
350  }
351  }
352 
353  op[ i ] = promote_pixel_converter::convert_from( sum == 0 ? value : value / sum );
354  }
355  }
356  }
357 
358  // 画像の中心部分の処理を行う
359  for( ; k + rd < d ; k += thread_numz )
360  {
361  for( j = thread_idy ; j < rh ; j += thread_numy )
362  {
363  opointer_type op = &out( 0, j, k );
364 
365  for( i = 0 ; i < w ; i++ )
366  {
367  promote_type value = promote_type( );
368  kvalue_type sum = 0;
369  difference_type sl = i < rw ? rw - i : 0;
370  difference_type el = w - i <= rw ? rw + w - i : fw;
371  difference_type sm = rh - j;
372 
373  for( difference_type n = 0 ; n < fd ; n++ )
374  {
375  for( difference_type m = sm ; m < fh ; m++ )
376  {
377  for( difference_type l = sl ; l < el ; l++ )
378  {
379  kvalue_type kv = kernel( l, m, n );
380  value += kv * in( i + l - rw, j + m - rh, k + n - rd );
381  sum += kv;
382  }
383  }
384  }
385 
386  op[ i ] = promote_pixel_converter::convert_from( sum == 0 ? value : value / sum );
387  }
388  }
389 
390  // 中心部分の処理を行う
391  for( ; j + rh < h ; j += thread_numy )
392  {
393  ipointer_type ip = &in( 0, j, k );
394  opointer_type op = &out( 0, j, k );
395 
396  for( i = 0 ; i < rw ; i++ )
397  {
398  difference_type *ppf = pf;
399  ipointer_type p = ip + i;
400  kpointer_type kp = &kernel[ 0 ];
401  promote_type value = promote_type( );
402  kvalue_type sum = 0;
403  difference_type lnum = fh * fd;
404  difference_type sn = rw - i;
405 
406  for( difference_type l = 0 ; l < lnum ; l++ )
407  {
408  for( difference_type n = sn ; n < fw ; n++ )
409  {
410  value += kp[ n ] * p[ ppf[ n ] ];
411  sum += kp[ n ];
412  }
413 
414  kp += fw;
415  ppf += fw;
416  }
417 
418  op[ i ] = promote_pixel_converter::convert_from( sum == 0 ? value : value / sum );
419  }
420 
421  for( ; i + rw < w ; i++ )
422  {
423  ipointer_type p = ip + i;
424  promote_type value = promote_type( );
425 
426  for( size_type n = 0 ; n < fsize ; n++ )
427  {
428  value += kernel[ n ] * p[ pf[ n ] ];
429  }
430 
431  op[ i ] = promote_pixel_converter::convert_from( value );
432  }
433 
434  for( ; i < w ; i++ )
435  {
436  difference_type *ppf = pf;
437  ipointer_type p = ip + i;
438  kpointer_type kp = &kernel[ 0 ];
439  promote_type value = promote_type( );
440  kvalue_type sum = 0;
441  difference_type lnum = fh * fd;
442  difference_type en = rw + w - i;
443 
444  for( difference_type l = 0 ; l < lnum ; l++ )
445  {
446  for( difference_type n = 0 ; n < en ; n++ )
447  {
448  value += kp[ n ] * p[ ppf[ n ] ];
449  sum += kp[ n ];
450  }
451 
452  kp += fw;
453  ppf += fw;
454  }
455 
456  op[ i ] = promote_pixel_converter::convert_from( sum == 0 ? value : value / sum );
457  }
458 
459  if( bprogress1 )
460  {
461  f( static_cast< double >( j + 1 ) / static_cast< double >( h ) * 100.0 );
462  }
463  }
464 
465  for( ; j < h ; j += thread_numy )
466  {
467  opointer_type op = &out( 0, j, k );
468 
469  for( i = 0 ; i < w ; i++ )
470  {
471  promote_type value = promote_type( );
472  kvalue_type sum = 0;
473  difference_type sl = i < rw ? rw - i : 0;
474  difference_type el = w - i <= rw ? rw + w - i : fw;
475  difference_type em = rh + h - j;
476 
477  for( difference_type n = 0 ; n < fd ; n++ )
478  {
479  for( difference_type m = 0 ; m < em ; m++ )
480  {
481  for( difference_type l = sl ; l < el ; l++ )
482  {
483  kvalue_type kv = kernel( l, m, n );
484  value += kv * in( i + l - rw, j + m - rh, k + n - rd );
485  sum += kv;
486  }
487  }
488  }
489 
490  op[ i ] = promote_pixel_converter::convert_from( sum == 0 ? value : value / sum );
491  }
492  }
493 
494  if( bprogress2 )
495  {
496  f( static_cast< double >( k + 1 ) / static_cast< double >( d ) * 100.0 );
497  }
498  }
499 
500  // 画像の縁の処理を行う
501  for( ; k < d ; k += thread_numz )
502  {
503  for( j = thread_idy ; j < h ; j += thread_numy )
504  {
505  opointer_type op = &out( 0, j, k );
506 
507  for( i = 0 ; i < w ; i++ )
508  {
509  promote_type value = promote_type( );
510  kvalue_type sum = 0;
511  difference_type sl = i < rw ? rw - i : 0;
512  difference_type el = w - i <= rw ? rw + w - i : fw;
513  difference_type sm = j < rh ? rh - j : 0;
514  difference_type em = h - j <= rh ? rh + h - j : fh;
515  difference_type en = d - k <= rd ? rd + d - k : fd;
516 
517  for( difference_type n = 0 ; n < en ; n++ )
518  {
519  for( difference_type m = sm ; m < em ; m++ )
520  {
521  for( difference_type l = sl ; l < el ; l++ )
522  {
523  kvalue_type kv = kernel( l, m, n );
524  value += kv * in( i + l - rw, j + m - rh, k + n - rd );
525  sum += kv;
526  }
527  }
528  }
529 
530  op[ i ] = promote_pixel_converter::convert_from( sum == 0 ? value : value / sum );
531  }
532  }
533  }
534 
535  delete [] pf;
536  }
537 
538  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
539  static void __linear_filter__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out, const Kernel &kernel,
540  typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num, Functor f )
541  {
542  linear_filter( in, out, kernel, 0, 1, thread_id, thread_num, f );
543  }
544 
545  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
546  static void __linear_filter__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, const Kernel &kernel,
547  typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num, Functor f )
548  {
549  linear_filter( in, out, kernel, thread_id, thread_num, 0, 1, f );
550  }
551 
552  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
553  static void __linear_filter__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, const Kernel &kernel,
554  typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num, Functor f )
555  {
556  linear_filter( in, out, kernel, 0, 1, thread_id, thread_num, f );
557  }
558 
559  template < class T1, class T2, class Kernel, class Functor >
560  class linear_thread : public mist::thread< linear_thread< T1, T2, Kernel, Functor > >
561  {
562  public:
563  typedef mist::thread< linear_thread< T1, T2, Kernel, Functor > > base;
564  typedef typename base::thread_exit_type thread_exit_type;
565  typedef typename T1::size_type size_type;
566  typedef typename T1::value_type value_type;
567 
568  private:
569  size_t thread_id_;
570  size_t thread_num_;
571 
572  // 入出力用の画像へのポインタ
573  const T1 *in_;
574  T2 *out_;
575  const Kernel *kernel_;
576 
577  Functor f_;
578 
579  public:
580  void setup_parameters( const T1 &in, T2 &out, const Kernel &kernel, size_type thread_id, size_type thread_num, Functor f )
581  {
582  in_ = &in;
583  out_ = &out;
584  kernel_ = &kernel;
585  thread_id_ = thread_id;
586  thread_num_ = thread_num;
587  f_ = f;
588  }
589 
590  linear_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ),
591  in_( NULL ), out_( NULL ), kernel_( NULL )
592  {
593  }
594 
595  linear_thread( const linear_thread &p ) : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ),
596  in_( p.in_ ), out_( p.out_ ), kernel_( p.kernel_ )
597  {
598  }
599 
600  protected:
601  // 継承した先で必ず実装されるスレッド関数
602  virtual thread_exit_type thread_function( )
603  {
604  __linear_filter__( *in_, *out_, *kernel_, thread_id_, thread_num_, f_ );
605  return( true );
606  }
607  };
608 
609 
610  template < int DIMENSION >
611  struct _1D_linear_filter_
612  {
613  template < class Array1, class Array2, class Kernel, class Functor >
614  static void linear_filter( const Array1 &in, Array2 &out, const Kernel &kernel,
615  typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
616  typename Array1::size_type thread_idz, typename Array1::size_type thread_numz, Functor f, double sp, double ep )
617  {
618  typedef typename Array1::size_type size_type;
619  typedef typename Array1::value_type value_type;
620  typedef typename Kernel::value_type kvalue_type;
621  typedef typename Array1::const_pointer ipointer;
622  typedef typename Array2::pointer opointer;
623  typedef typename Array1::difference_type difference_type;
624  typedef __promote_pixel_converter_< typename Array2::value_type > promote_pixel_converter;
625  typedef typename promote_pixel_converter::promote_type promote_type;
626  typedef __access__< DIMENSION > access;
627 
628  difference_type _1e = access::size1( in );
629  difference_type _2e = access::size2( in );
630  difference_type _3e = access::size3( in );
631 
632  const bool bprogress1 = thread_idy == 0 && _3e == 1;
633  const bool bprogress2 = thread_idz == 0 && _3e > 1;
634 
635  difference_type fw = kernel.size( );
636  difference_type rw = fw / 2;
637 
638  difference_type i1, i2, i3;
639  value_type *tmp = new value_type[ _1e ];
640  difference_type diff = &access::at( in, 1, 0, 0 ) - &access::at( in, 0, 0, 0 );
641 
642  for( i3 = thread_idz ; i3 < _3e ; i3 += thread_numz )
643  {
644  for( i2 = thread_idy ; i2 < _2e ; i2 += thread_numy )
645  {
646  ipointer p0 = &access::at( in, 0, i2, i3 );
647  opointer op = &access::at( out, 0, i2, i3 );
648  ipointer p = p0;
649 
650  // テンポラリ領域に画素値を一時的に記憶する
651  for( i1 = 0 ; i1 + rw + 1 < fw ; i1++ )
652  {
653  tmp[ i1 ] = *p;
654  p += diff;
655  }
656 
657  for( i1 = 0 ; i1 < rw ; i1++ )
658  {
659  tmp[ i1 + rw ] = *p;
660  p += diff;
661 
662  promote_type value = promote_type( );
663 
664  value_type *ptmp = tmp + i1 - rw;
665  kvalue_type sum = 0;
666  for( size_type n = rw - i1 ; n < kernel.size( ) ; n++ )
667  {
668  value += kernel[ n ] * ptmp[ n ];
669  sum += kernel[ n ];
670  }
671 
672  *op = promote_pixel_converter::convert_from( sum == 0 ? value : value / sum );
673  op += diff;
674  }
675 
676  for( ; i1 + rw < _1e ; i1++ )
677  {
678  tmp[ i1 + rw ] = *p;
679  p += diff;
680 
681  promote_type value = promote_type( );
682 
683  value_type *ptmp = tmp + i1 - rw;
684  for( size_type n = 0 ; n < kernel.size( ) ; n++ )
685  {
686  value += kernel[ n ] * ptmp[ n ];
687  }
688 
689  *op = promote_pixel_converter::convert_from( value );
690  op += diff;
691  }
692 
693  for( ; i1 < _1e ; i1++ )
694  {
695  promote_type value = promote_type( );
696 
697  value_type *ptmp = tmp + i1 - rw;
698  kvalue_type sum = 0;
699  size_type ne = rw + _1e - i1;
700  for( size_type n = 0 ; n < ne ; n++ )
701  {
702  value += kernel[ n ] * ptmp[ n ];
703  sum += kernel[ n ];
704  }
705 
706  *op = promote_pixel_converter::convert_from( sum == 0 ? value : value / sum );
707  op += diff;
708  }
709 
710  if( bprogress1 )
711  {
712  f( sp + static_cast< double >( i2 + 1 ) / static_cast< double >( _2e ) * ( ep - sp ) );
713  }
714  }
715 
716  if( bprogress2 )
717  {
718  f( sp + static_cast< double >( i3 + 1 ) / static_cast< double >( _3e ) * ( ep - sp ) );
719  }
720  }
721 
722  delete [] tmp;
723  }
724  };
725 
726  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
727  static void __1D_linear_filter__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out, const Kernel &k1, const Kernel & /* k2 */, const Kernel & /* k3 */,
728  typename array3< T1, Allocator1 >::size_type axis, typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num, Functor f )
729  {
730  linear_filter( in, out, k1, 0, 1, thread_id, thread_num, f );
731  }
732 
733  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
734  static void __1D_linear_filter__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, const Kernel &k1, const Kernel &k2, const Kernel & /* k3 */,
735  typename array3< T1, Allocator1 >::size_type axis, typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num, Functor f )
736  {
737  if( axis == 0 )
738  {
739  _1D_linear_filter_< 1 >::linear_filter( in, out, k1, thread_id, thread_num, 0, 1, f, 0.0, 50.0 );
740  }
741  else
742  {
743  _1D_linear_filter_< 2 >::linear_filter( out, out, k2, thread_id, thread_num, 0, 1, f, 50.0, 100.0 );
744  }
745  }
746 
747  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
748  static void __1D_linear_filter__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, const Kernel &k1, const Kernel &k2, const Kernel &k3,
749  typename array3< T1, Allocator1 >::size_type axis, typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num, Functor f )
750  {
751  if( axis == 0 )
752  {
753  _1D_linear_filter_< 1 >::linear_filter( in, out, k1, 0, 1, thread_id, thread_num, f, 0.0, 100.0 / 3.0 );
754  }
755  else if( axis == 1 )
756  {
757  _1D_linear_filter_< 2 >::linear_filter( out, out, k2, 0, 1, thread_id, thread_num, f, 100.0 / 3.0, 200.0 / 3.0 );
758  }
759  else
760  {
761  _1D_linear_filter_< 3 >::linear_filter( out, out, k3, 0, 1, thread_id, thread_num, f, 200.0 / 3.0, 100.0 );
762  }
763  }
764 
765  template < class T1, class T2, class Kernel, class Functor >
766  class _1D_linear_thread : public mist::thread< _1D_linear_thread< T1, T2, Kernel, Functor > >
767  {
768  public:
769  typedef mist::thread< _1D_linear_thread< T1, T2, Kernel, Functor > > base;
770  typedef typename base::thread_exit_type thread_exit_type;
771  typedef typename T1::size_type size_type;
772  typedef typename T1::value_type value_type;
773 
774  private:
775  size_t thread_id_;
776  size_t thread_num_;
777 
778  // 入出力用の画像へのポインタ
779  const T1 *in_;
780  T2 *out_;
781  const Kernel *k1_;
782  const Kernel *k2_;
783  const Kernel *k3_;
784  size_type axis_;
785 
786  Functor f_;
787 
788  public:
789  void setup_parameters( const T1 &in, T2 &out, const Kernel &k1, const Kernel &k2, const Kernel &k3, size_type axis, size_type thread_id, size_type thread_num, Functor f )
790  {
791  in_ = &in;
792  out_ = &out;
793  k1_ = &k1;
794  k2_ = &k2;
795  k3_ = &k3;
796  axis_ = axis;
797  thread_id_ = thread_id;
798  thread_num_ = thread_num;
799  f_ = f;
800  }
801 
802  _1D_linear_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ),
803  in_( NULL ), out_( NULL ), k1_( NULL ), k2_( NULL ), k3_( NULL )
804  {
805  }
806 
807  _1D_linear_thread( const _1D_linear_thread &p ) : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ),
808  in_( p.in_ ), out_( p.out_ ), k1_( p.k1_ ), k2_( p.k2_ ), k3_( p.k3_ )
809  {
810  }
811 
812  protected:
813  // 継承した先で必ず実装されるスレッド関数
814  virtual thread_exit_type thread_function( )
815  {
816  __1D_linear_filter__( *in_, *out_, *k1_, *k2_, *k3_, axis_, thread_id_, thread_num_, f_ );
817  return( true );
818  }
819  };
820 
821  // in : 入力画像. 入力画像の画素値は min と max の間とする
822  // out : 出力画像. 出力画像のメモリはあらかじめ割り当てられているものとする
823  // fw, fh, fd : フィルタサイズ
824  template < class Array1, class Array2, class Functor >
825  void average_filter( const Array1 &in, Array2 &out,
826  typename Array1::size_type fw, typename Array1::size_type fh, typename Array1::size_type fd,
827  typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
828  typename Array1::size_type thread_idz, typename Array1::size_type thread_numz, Functor f )
829  {
830  typedef typename Array1::size_type size_type;
831  typedef typename Array1::difference_type difference_type;
832  typedef __promote_pixel_converter_< typename Array2::value_type > promote_pixel_converter;
833  typedef typename promote_pixel_converter::promote_type promote_type;
834  typedef typename Array2::value_type out_value_type;
835  typedef typename Array1::const_pointer ipointer_type;
836  typedef typename Array2::pointer opointer_type;
837 
838  difference_type w = in.width( );
839  difference_type h = in.height( );
840  difference_type d = in.depth( );
841 
842  const bool bprogress1 = thread_idy == 0 && d == 1;
843  const bool bprogress2 = thread_idz == 0 && d > 1;
844 
845  difference_type rw = fw / 2;
846  difference_type rh = fh / 2;
847  difference_type rd = fd / 2;
848  difference_type rrw = fw - rw;
849 
850  difference_type i, j, k;
851 
852  promote_type *sum = new promote_type[ fw ];
853 
854  for( k = thread_idz ; k < d ; k += thread_numz )
855  {
856  difference_type sz = k - rd;
857  difference_type ez = sz + fd - 1;
858 
859  if( sz < 0 )
860  {
861  sz = 0;
862  }
863  if( ez >= d )
864  {
865  ez = d - 1;
866  }
867 
868  difference_type fdd = ez - sz + 1;
869 
870  for( j = thread_idy ; j < h ; j += thread_numy )
871  {
872  difference_type sy = j - rh;
873  difference_type ey = sy + fh - 1;
874 
875  if( sy < 0 )
876  {
877  sy = 0;
878  }
879  if( ey >= h )
880  {
881  ey = h - 1;
882  }
883 
884  difference_type fhh = ey - sy + 1;
885 
886  // キャッシュ用の総和を計算しておく
887  promote_type __sum__ = promote_type( );
888  for( i = 0 ; i < rrw - 1 ; i++ )
889  {
890  promote_type value = promote_type( );
891  for( difference_type z = sz ; z <= ez ; z++ )
892  {
893  for( difference_type y = sy ; y <= ey ; y++ )
894  {
895  value += in( i, y, z );
896  }
897  }
898 
899  __sum__ += value;
900  sum[ i ] = value;
901  }
902 
903  double fsize = static_cast< double >( fhh * fdd );
904  opointer_type op = &out( 0, j, k );
905  size_type indx = rrw - 1;
906 
907  // 左端の処理を済ませておく
908  for( i = 0 ; i <= rw ; i++, indx++ )
909  {
910  promote_type value = promote_type( );
911  for( difference_type z = sz ; z <= ez ; z++ )
912  {
913  for( difference_type y = sy ; y <= ey ; y++ )
914  {
915  value += in( indx, y, z );
916  }
917  }
918 
919  __sum__ += value;
920  sum[ indx ] = value;
921 
922  op[ i ] = promote_pixel_converter::convert_from( __sum__ / ( fsize * ( rrw + i ) ) );
923  }
924 
925  for( ; i + rrw <= w ; i++, indx++ )
926  {
927  if( indx >= fw )
928  {
929  indx = 0;
930  }
931 
932  __sum__ -= sum[ indx ];
933 
934  promote_type value = promote_type( );
935  for( difference_type z = sz ; z <= ez ; z++ )
936  {
937  for( difference_type y = sy ; y <= ey ; y++ )
938  {
939  value += in( i + rrw - 1, y, z );
940  }
941  }
942 
943  __sum__ += value;
944  sum[ indx ] = value;
945 
946  op[ i ] = promote_pixel_converter::convert_from( __sum__ / ( fsize * fw ) );
947  }
948 
949  for( ; i < w ; i++, indx++ )
950  {
951  if( indx >= fw )
952  {
953  indx = 0;
954  }
955 
956  __sum__ -= sum[ indx ];
957 
958  op[ i ] = promote_pixel_converter::convert_from( __sum__ / ( fsize * ( rw + w - i ) ) );
959  }
960 
961  if( bprogress1 )
962  {
963  f( static_cast< double >( j + 1 ) / static_cast< double >( h ) * 100.0 );
964  }
965  }
966 
967  if( bprogress2 )
968  {
969  f( static_cast< double >( k + 1 ) / static_cast< double >( d ) * 100.0 );
970  }
971  }
972 
973  delete [] sum;
974  }
975 
976 
977  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
978  static void __average_filter__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
979  typename array< T1, Allocator1 >::size_type fw, typename array< T1, Allocator1 >::size_type fh, typename array< T1, Allocator1 >::size_type fd,
980  typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num, Functor f )
981  {
982  average_filter( in, out, fw, fh, fd, 0, 1, thread_id, thread_num, f );
983  }
984 
985  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
986  static void __average_filter__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
987  typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh, typename array2< T1, Allocator1 >::size_type fd,
988  typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num, Functor f )
989  {
990  average_filter( in, out, fw, fh, fd, thread_id, thread_num, 0, 1, f );
991  }
992 
993  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
994  static void __average_filter__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
995  typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
996  typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num, Functor f )
997  {
998  average_filter( in, out, fw, fh, fd, 0, 1, thread_id, thread_num, f );
999  }
1000 
1001  template < class T1, class T2, class Functor >
1002  class average_thread : public mist::thread< average_thread< T1, T2, Functor > >
1003  {
1004  public:
1005  typedef mist::thread< average_thread< T1, T2, Functor > > base;
1006  typedef typename base::thread_exit_type thread_exit_type;
1007  typedef typename T1::size_type size_type;
1008  typedef typename T1::value_type value_type;
1009 
1010  private:
1011  size_t thread_id_;
1012  size_t thread_num_;
1013 
1014  // 入出力用の画像へのポインタ
1015  const T1 *in_;
1016  T2 *out_;
1017  size_type fw_;
1018  size_type fh_;
1019  size_type fd_;
1020 
1021  Functor f_;
1022 
1023  public:
1024  void setup_parameters( const T1 &in, T2 &out, size_type fw, size_type fh, size_type fd, size_type thread_id, size_type thread_num, Functor f )
1025  {
1026  in_ = &in;
1027  out_ = &out;
1028  fw_ = fw;
1029  fh_ = fh;
1030  fd_ = fd;
1031  thread_id_ = thread_id;
1032  thread_num_ = thread_num;
1033  f_ = f;
1034  }
1035 
1036  average_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ),
1037  in_( NULL ), out_( NULL ), fw_( 0 ), fh_( 0 ), fd_( 0 )
1038  {
1039  }
1040 
1041  average_thread( const average_thread &p ) : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ),
1042  in_( p.in_ ), out_( p.out_ ), fw_( p.fw_ ), fh_( p.fh_ ), fd_( p.fd_ )
1043  {
1044  }
1045 
1046  protected:
1047  // 継承した先で必ず実装されるスレッド関数
1048  virtual thread_exit_type thread_function( )
1049  {
1050  __average_filter__( *in_, *out_, fw_, fh_, fd_, thread_id_, thread_num_, f_ );
1051  return( true );
1052  }
1053  };
1054 
1055  template < int DIMENSION >
1056  struct _1D_average_filter_
1057  {
1058  template < class Array1, class Array2, class Functor >
1059  static void average_filter( const Array1 &in, Array2 &out, typename Array1::difference_type fw,
1060  typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
1061  typename Array1::size_type thread_idz, typename Array1::size_type thread_numz, Functor f, double sp, double ep )
1062  {
1063  typedef typename Array1::size_type size_type;
1064  typedef typename Array1::value_type value_type;
1065  typedef typename Array1::const_pointer ipointer;
1066  typedef typename Array2::pointer opointer;
1067  typedef typename Array1::difference_type difference_type;
1068  typedef __promote_pixel_converter_< typename Array2::value_type > promote_pixel_converter;
1069  typedef typename promote_pixel_converter::promote_type promote_type;
1070  typedef __access__< DIMENSION > access;
1071 
1072  difference_type _1e = access::size1( in );
1073  difference_type _2e = access::size2( in );
1074  difference_type _3e = access::size3( in );
1075 
1076  const bool bprogress1 = thread_idy == 0 && _3e == 1;
1077  const bool bprogress2 = thread_idz == 0 && _3e > 1;
1078 
1079  difference_type rw = fw / 2;
1080 
1081  difference_type i1, i2, i3;
1082  value_type *tmp = new value_type[ _1e ];
1083  difference_type diff = &access::at( in, 1, 0, 0 ) - &access::at( in, 0, 0, 0 );
1084 
1085  for( i3 = thread_idz ; i3 < _3e ; i3 += thread_numz )
1086  {
1087  for( i2 = thread_idy ; i2 < _2e ; i2 += thread_numy )
1088  {
1089  ipointer p0 = &access::at( in, 0, i2, i3 );
1090  opointer op = &access::at( out, 0, i2, i3 );
1091  ipointer p = p0;
1092 
1093  promote_type value = promote_type( );
1094 
1095  // テンポラリ領域に画素値を一時的に記憶する
1096  for( i1 = 0 ; i1 + rw + 1 < fw ; i1++ )
1097  {
1098  value += *p;
1099  tmp[ i1 ] = *p;
1100  p += diff;
1101  }
1102 
1103  for( i1 = 0 ; i1 <= rw ; i1++ )
1104  {
1105  value += *p;
1106  tmp[ i1 + rw ] = *p;
1107  p += diff;
1108 
1109  *op = promote_pixel_converter::convert_from( value / static_cast< double >( fw - rw + i1 ) );
1110  op += diff;
1111  }
1112 
1113  for( ; i1 + rw < _1e ; i1++ )
1114  {
1115  value += *p;
1116  value -= tmp[ i1 - rw - 1 ];
1117  tmp[ i1 + rw ] = *p;
1118  p += diff;
1119 
1120  *op = promote_pixel_converter::convert_from( value / static_cast< double >( fw ) );
1121  op += diff;
1122  }
1123 
1124  for( ; i1 < _1e ; i1++ )
1125  {
1126  value -= tmp[ i1 - rw - 1 ];
1127 
1128  *op = promote_pixel_converter::convert_from( value / static_cast< double >( rw + _1e - i1 ) );
1129  op += diff;
1130  }
1131 
1132  if( bprogress1 )
1133  {
1134  f( sp + static_cast< double >( i2 + 1 ) / static_cast< double >( _2e ) * ( ep - sp ) );
1135  }
1136  }
1137 
1138  if( bprogress2 )
1139  {
1140  f( sp + static_cast< double >( i3 + 1 ) / static_cast< double >( _3e ) * ( ep - sp ) );
1141  }
1142  }
1143 
1144  delete [] tmp;
1145  }
1146  };
1147 
1148  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
1149  static void __1D_average_filter__( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
1150  typename array< T1, Allocator1 >::size_type fw, typename array< T1, Allocator1 >::size_type fh, typename array< T1, Allocator1 >::size_type fd,
1151  typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num, Functor f )
1152  {
1153  average_filter( in, out, fw, fh, fd, 0, 1, thread_id, thread_num, f );
1154  }
1155 
1156  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
1157  static void __1D_average_filter__( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
1158  typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh, typename array2< T1, Allocator1 >::size_type fd,
1159  typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num, Functor f )
1160  {
1161  _1D_average_filter_< 1 >::average_filter( in, out, fw, thread_id, thread_num, 0, 1, f, 0.0, 50.0 );
1162  _1D_average_filter_< 2 >::average_filter( out, out, fh, thread_id, thread_num, 0, 1, f, 50.0, 100.0 );
1163  }
1164 
1165  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
1166  static void __1D_average_filter__( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
1167  typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
1168  typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num, Functor f )
1169  {
1170  _1D_average_filter_< 1 >::average_filter( in, out, fw, 0, 1, thread_id, thread_num, f, 0.0, 100.0 / 3.0 );
1171  _1D_average_filter_< 2 >::average_filter( out, out, fh, 0, 1, thread_id, thread_num, f, 100.0 / 3.0, 200.0 / 3.0 );
1172  _1D_average_filter_< 3 >::average_filter( out, out, fd, 0, 1, thread_id, thread_num, f, 200.0 / 3.0, 100.0 );
1173  }
1174 
1175  template < class T1, class T2, class Functor >
1176  class _1D_average_thread : public mist::thread< _1D_average_thread< T1, T2, Functor > >
1177  {
1178  public:
1179  typedef mist::thread< _1D_average_thread< T1, T2, Functor > > base;
1180  typedef typename base::thread_exit_type thread_exit_type;
1181  typedef typename T1::size_type size_type;
1182  typedef typename T1::value_type value_type;
1183 
1184  private:
1185  size_t thread_id_;
1186  size_t thread_num_;
1187 
1188  // 入出力用の画像へのポインタ
1189  const T1 *in_;
1190  T2 *out_;
1191  size_type fw_;
1192  size_type fh_;
1193  size_type fd_;
1194 
1195  Functor f_;
1196 
1197  public:
1198  void setup_parameters( const T1 &in, T2 &out, size_type fw, size_type fh, size_type fd, size_type thread_id, size_type thread_num, Functor f )
1199  {
1200  in_ = &in;
1201  out_ = &out;
1202  fw_ = fw;
1203  fh_ = fh;
1204  fd_ = fd;
1205  thread_id_ = thread_id;
1206  thread_num_ = thread_num;
1207  f_ = f;
1208  }
1209 
1210  _1D_average_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ),
1211  in_( NULL ), out_( NULL ), fw_( 0 ), fh_( 0 ), fd_( 0 )
1212  {
1213  }
1214 
1215  _1D_average_thread( const _1D_average_thread &p ) : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ),
1216  in_( p.in_ ), out_( p.out_ ), fw_( p.fw_ ), fh_( p.fh_ ), fd_( p.fd_ )
1217  {
1218  }
1219 
1220  protected:
1221  // 継承した先で必ず実装されるスレッド関数
1222  virtual thread_exit_type thread_function( )
1223  {
1224  __1D_average_filter__( *in_, *out_, fw_, fh_, fd_, thread_id_, thread_num_, f_ );
1225  return( true );
1226  }
1227  };
1228 
1229 
1230  template < class Array >
1231  inline void compute_normalized_kernel( Array &kernel )
1232  {
1233  typedef typename Array::value_type value_type;
1234  typedef typename Array::size_type size_type;
1235 
1236  value_type sum = value_type( );
1237  for( size_type i = 0 ; i < kernel.size( ) ; i++ )
1238  {
1239  sum += kernel[ i ];
1240  }
1241 
1242  for( size_type i = 0 ; i < kernel.size( ) ; i++ )
1243  {
1244  kernel[ i ] /= sum;
1245  }
1246  }
1247 
1248  template < class Array >
1249  inline void compute_gaussian_kernel( Array &kernel, double sigma )
1250  {
1251  typedef typename Array::size_type size_type;
1252  typedef typename Array::difference_type difference_type;
1253  typedef typename Array::value_type value_type;
1254 
1255  if( sigma == 0.0 )
1256  {
1257  sigma = 1.0;
1258  }
1259 
1260  double radius = sigma * 2.0;
1261  double ax = kernel.reso1( );
1262  double ay = kernel.reso2( );
1263  double az = kernel.reso3( );
1264  double axx = ax * ax;
1265  double ayy = ay * ay;
1266  double azz = az * az;
1267 
1268  difference_type Rx = static_cast< difference_type >( std::ceil( radius / ax ) );
1269  difference_type Ry = static_cast< difference_type >( std::ceil( radius / ay ) );
1270  difference_type Rz = static_cast< difference_type >( std::ceil( radius / az ) );
1271 
1272  kernel.resize( Rx * 2 + 1, Ry * 2 + 1, Rz * 2 + 1 );
1273 
1274  difference_type fw = kernel.width( );
1275  difference_type fh = kernel.height( );
1276  difference_type fd = kernel.depth( );
1277 
1278  difference_type rw = fw / 2;
1279  difference_type rh = fh / 2;
1280  difference_type rd = fd / 2;
1281 
1282  double _1_sig2 = 1.0 / ( sigma * sigma * 2.0 );
1283  double sum = 0.0;
1284  for( difference_type z = -rd ; z <= rd ; z++ )
1285  {
1286  double zz = static_cast< double >( z * z * azz );
1287  for( difference_type y = -rh ; y <= rh ; y++ )
1288  {
1289  double yy = static_cast< double >( y * y * ayy );
1290  for( difference_type x = -rw ; x <= rw ; x++ )
1291  {
1292  double xx = static_cast< double >( x * x * axx );
1293  double g = std::exp( -( xx + yy + zz ) * _1_sig2 );
1294  sum += g;
1295  kernel( x + rw, y + rh, z + rd ) = static_cast< value_type >( g );
1296  }
1297  }
1298  }
1299 
1300  for( size_type i = 0 ; i < kernel.size( ) ; i++ )
1301  {
1302  kernel[ i ] /= sum;
1303  }
1304  }
1305 
1306  template < class Array >
1307  inline void compute_gaussian_kernel( Array &k1, Array &k2, Array &k3, double sigma )
1308  {
1309  typedef typename Array::size_type size_type;
1310  typedef typename Array::difference_type difference_type;
1311  typedef typename Array::value_type value_type;
1312 
1313  if( sigma == 0.0 )
1314  {
1315  sigma = 1.0;
1316  }
1317 
1318  double radius = sigma * 2.0;
1319  double ax = k1.reso1( );
1320  double ay = k2.reso1( );
1321  double az = k3.reso1( );
1322  double axx = ax * ax;
1323  double ayy = ay * ay;
1324  double azz = az * az;
1325 
1326  difference_type Rx = static_cast< difference_type >( std::ceil( radius / ax ) );
1327  difference_type Ry = static_cast< difference_type >( std::ceil( radius / ay ) );
1328  difference_type Rz = static_cast< difference_type >( std::ceil( radius / az ) );
1329 
1330  k1.resize( Rx * 2 + 1 );
1331  k2.resize( Ry * 2 + 1 );
1332  k3.resize( Rz * 2 + 1 );
1333 
1334  difference_type fw = k1.size( );
1335  difference_type fh = k2.size( );
1336  difference_type fd = k3.size( );
1337 
1338  difference_type rw = fw / 2;
1339  difference_type rh = fh / 2;
1340  difference_type rd = fd / 2;
1341 
1342  double _1_sig2 = 1.0 / ( sigma * sigma * 2.0 );
1343  {
1344  double sum = 0.0;
1345  for( difference_type x = -rw ; x <= rw ; x++ )
1346  {
1347  double xx = static_cast< double >( x * x * axx );
1348  double g = std::exp( -xx * _1_sig2 );
1349  sum += g;
1350  k1[ x + rw ] = static_cast< value_type >( g );
1351  }
1352 
1353  for( size_type i = 0 ; i < k1.size( ) ; i++ )
1354  {
1355  k1[ i ] /= sum;
1356  }
1357  }
1358 
1359  {
1360  double sum = 0.0;
1361  for( difference_type y = -rh ; y <= rh ; y++ )
1362  {
1363  double yy = static_cast< double >( y * y * ayy );
1364  double g = std::exp( -yy * _1_sig2 );
1365  sum += g;
1366  k2[ y + rh ] = static_cast< value_type >( g );
1367  }
1368 
1369  for( size_type i = 0 ; i < k2.size( ) ; i++ )
1370  {
1371  k2[ i ] /= sum;
1372  }
1373  }
1374 
1375  {
1376  double sum = 0.0;
1377  for( difference_type z = -rd ; z <= rd ; z++ )
1378  {
1379  double zz = static_cast< double >( z * z * azz );
1380  double g = std::exp( -zz * _1_sig2 );
1381  sum += g;
1382  k3[ z + rd ] = static_cast< value_type >( g );
1383  }
1384 
1385  for( size_type i = 0 ; i < k3.size( ) ; i++ )
1386  {
1387  k3[ i ] /= sum;
1388  }
1389  }
1390  }
1391 }
1392 
1393 
1394 
1395 
1396 
1426 
1427 
1429 namespace linear
1430 {
1449  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
1450  bool filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out, const Kernel &kernel, Functor f, typename array< T1, Allocator1 >::size_type thread_num )
1451  {
1452  if( is_same_object( in, out ) || in.empty( ) )
1453  {
1454  return( false );
1455  }
1456 
1457  out.resize( in.size( ) );
1458 
1459  __linear__::__linear_filter__( in, out, kernel, 0, 1, __mist_dmy_callback__( ) );
1460 
1461  return( true );
1462  }
1463 
1464 
1483  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
1484  bool filter( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out, const Kernel &kernel, Functor f, typename array1< T1, Allocator1 >::size_type thread_num )
1485  {
1486  if( is_same_object( in, out ) || in.empty( ) )
1487  {
1488  return( false );
1489  }
1490 
1491  out.resize( in.size( ) );
1492  out.reso1( in.reso1( ) );
1493 
1494  __linear__::__linear_filter__( in, out, kernel, 0, 1, __mist_dmy_callback__( ) );
1495 
1496  return( true );
1497  }
1498 
1499 
1518  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
1519  bool filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, const Kernel &kernel, Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
1520  {
1521  if( is_same_object( in, out ) || in.empty( ) )
1522  {
1523  return( false );
1524  }
1525 
1526  typedef typename array2< T1, Allocator1 >::size_type size_type;
1527  typedef __linear__::linear_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 >, Kernel, Functor > linear_thread;
1528 
1529  if( thread_num == 0 )
1530  {
1531  thread_num = static_cast< size_type >( get_cpu_num( ) );
1532  }
1533 
1534  out.resize( in.size1( ), in.size2( ) );
1535  out.reso1( in.reso1( ) );
1536  out.reso2( in.reso2( ) );
1537 
1538  linear_thread *thread = new linear_thread[ thread_num ];
1539 
1540  size_type i;
1541  for( i = 0 ; i < thread_num ; i++ )
1542  {
1543  thread[ i ].setup_parameters( in, out, kernel, i, thread_num, f );
1544  }
1545 
1546  f( 0.0 );
1547 
1548  // スレッドを実行して，終了まで待機する
1549  do_threads_( thread, thread_num );
1550 
1551  f( 100.1 );
1552 
1553  delete [] thread;
1554 
1555  return( true );
1556  }
1557 
1558 
1579  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
1580  bool filter1d( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, const Kernel &kernel1, const Kernel &kernel2, Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
1581  {
1582  if( in.empty( ) )
1583  {
1584  return( false );
1585  }
1586 
1587  typedef typename array2< T1, Allocator1 >::size_type size_type;
1588  typedef __linear__::_1D_linear_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 >, Kernel, Functor > _1D_linear_thread;
1589 
1590  if( thread_num == 0 )
1591  {
1592  thread_num = static_cast< size_type >( get_cpu_num( ) );
1593  }
1594 
1595  out.resize( in.size1( ), in.size2( ) );
1596  out.reso1( in.reso1( ) );
1597  out.reso2( in.reso2( ) );
1598 
1599  _1D_linear_thread *thread = new _1D_linear_thread[ thread_num ];
1600 
1601  f( 0.0 );
1602 
1603  {
1604  for( size_type i = 0 ; i < thread_num ; i++ )
1605  {
1606  thread[ i ].setup_parameters( in, out, kernel1, kernel2, kernel2, 0, i, thread_num, f );
1607  }
1608 
1609  // スレッドを実行して，終了まで待機する
1610  do_threads_( thread, thread_num );
1611  }
1612 
1613  {
1614  for( size_type i = 0 ; i < thread_num ; i++ )
1615  {
1616  thread[ i ].setup_parameters( in, out, kernel1, kernel2, kernel2, 1, i, thread_num, f );
1617  }
1618 
1619  // スレッドを実行して，終了まで待機する
1620  do_threads_( thread, thread_num );
1621  }
1622 
1623  f( 100.1 );
1624 
1625  delete [] thread;
1626 
1627  return( true );
1628  }
1629 
1630 
1649  template < class T, class Allocator, class Kernel, class Functor >
1650  bool filter1d( array2< T, Allocator > &in, const Kernel &kernel1, const Kernel &kernel2, Functor f, typename array2< T, Allocator >::size_type thread_num )
1651  {
1652  if( in.empty( ) )
1653  {
1654  return( false );
1655  }
1656 
1657  typedef typename array2< T, Allocator >::size_type size_type;
1658  typedef __linear__::_1D_linear_thread< array2< T, Allocator >, array2< T, Allocator >, Kernel, Functor > _1D_linear_thread;
1659 
1660  if( thread_num == 0 )
1661  {
1662  thread_num = static_cast< size_type >( get_cpu_num( ) );
1663  }
1664 
1665  _1D_linear_thread *thread = new _1D_linear_thread[ thread_num ];
1666 
1667  f( 0.0 );
1668 
1669  {
1670  for( size_type i = 0 ; i < thread_num ; i++ )
1671  {
1672  thread[ i ].setup_parameters( in, in, kernel1, kernel2, kernel2, 0, i, thread_num, f );
1673  }
1674 
1675  // スレッドを実行して，終了まで待機する
1676  do_threads_( thread, thread_num );
1677  }
1678 
1679  {
1680  for( size_type i = 0 ; i < thread_num ; i++ )
1681  {
1682  thread[ i ].setup_parameters( in, in, kernel1, kernel2, kernel2, 1, i, thread_num, f );
1683  }
1684 
1685  // スレッドを実行して，終了まで待機する
1686  do_threads_( thread, thread_num );
1687  }
1688 
1689  f( 100.1 );
1690 
1691  delete [] thread;
1692 
1693  return( true );
1694  }
1695 
1696 
1715  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
1716  bool filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, const Kernel &kernel, Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
1717  {
1718  if( is_same_object( in, out ) || in.empty( ) )
1719  {
1720  return( false );
1721  }
1722 
1723  typedef typename array3< T1, Allocator1 >::size_type size_type;
1724  typedef __linear__::linear_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 >, Kernel, Functor > linear_thread;
1725 
1726  if( thread_num == 0 )
1727  {
1728  thread_num = static_cast< size_type >( get_cpu_num( ) );
1729  }
1730 
1731  out.resize( in.size1( ), in.size2( ), in.size3( ) );
1732  out.reso1( in.reso1( ) );
1733  out.reso2( in.reso2( ) );
1734  out.reso3( in.reso3( ) );
1735 
1736  linear_thread *thread = new linear_thread[ thread_num ];
1737 
1738  size_type i;
1739  for( i = 0 ; i < thread_num ; i++ )
1740  {
1741  thread[ i ].setup_parameters( in, out, kernel, i, thread_num, f );
1742  }
1743 
1744  f( 0.0 );
1745 
1746  // スレッドを実行して，終了まで待機する
1747  do_threads_( thread, thread_num );
1748 
1749  f( 100.1 );
1750 
1751  delete [] thread;
1752 
1753  return( true );
1754  }
1755 
1756 
1778  template < class T1, class Allocator1, class T2, class Allocator2, class Kernel, class Functor >
1779  bool filter1d( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, const Kernel &kernel1, const Kernel &kernel2, const Kernel &kernel3, Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
1780  {
1781  if( in.empty( ) )
1782  {
1783  return( false );
1784  }
1785 
1786  typedef typename array3< T1, Allocator1 >::size_type size_type;
1787  typedef __linear__::_1D_linear_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 >, Kernel, Functor > _1D_linear_thread;
1788 
1789  if( thread_num == 0 )
1790  {
1791  thread_num = static_cast< size_type >( get_cpu_num( ) );
1792  }
1793 
1794  out.resize( in.size1( ), in.size2( ), in.size3( ) );
1795  out.reso1( in.reso1( ) );
1796  out.reso2( in.reso2( ) );
1797  out.reso3( in.reso3( ) );
1798 
1799  _1D_linear_thread *thread = new _1D_linear_thread[ thread_num ];
1800 
1801  f( 0.0 );
1802 
1803  {
1804  for( size_type i = 0 ; i < thread_num ; i++ )
1805  {
1806  thread[ i ].setup_parameters( in, out, kernel1, kernel2, kernel3, 0, i, thread_num, f );
1807  }
1808 
1809  // スレッドを実行して，終了まで待機する
1810  do_threads_( thread, thread_num );
1811  }
1812 
1813  {
1814  for( size_type i = 0 ; i < thread_num ; i++ )
1815  {
1816  thread[ i ].setup_parameters( in, out, kernel1, kernel2, kernel3, 1, i, thread_num, f );
1817  }
1818 
1819  // スレッドを実行して，終了まで待機する
1820  do_threads_( thread, thread_num );
1821  }
1822 
1823  {
1824  for( size_type i = 0 ; i < thread_num ; i++ )
1825  {
1826  thread[ i ].setup_parameters( in, out, kernel1, kernel2, kernel3, 2, i, thread_num, f );
1827  }
1828 
1829  // スレッドを実行して，終了まで待機する
1830  do_threads_( thread, thread_num );
1831  }
1832 
1833  f( 100.1 );
1834 
1835  delete [] thread;
1836 
1837  return( true );
1838  }
1839 
1840 
1860  template < class T, class Allocator, class Kernel, class Functor >
1861  bool filter1d( array3< T, Allocator > &in, const Kernel &kernel1, const Kernel &kernel2, const Kernel &kernel3, Functor f, typename array3< T, Allocator >::size_type thread_num )
1862  {
1863  if( in.empty( ) )
1864  {
1865  return( false );
1866  }
1867 
1868  typedef typename array3< T, Allocator >::size_type size_type;
1869  typedef __linear__::_1D_linear_thread< array3< T, Allocator >, array3< T, Allocator >, Kernel, Functor > _1D_linear_thread;
1870 
1871  if( thread_num == 0 )
1872  {
1873  thread_num = static_cast< size_type >( get_cpu_num( ) );
1874  }
1875 
1876  _1D_linear_thread *thread = new _1D_linear_thread[ thread_num ];
1877 
1878  f( 0.0 );
1879 
1880  {
1881  for( size_type i = 0 ; i < thread_num ; i++ )
1882  {
1883  thread[ i ].setup_parameters( in, in, kernel1, kernel2, kernel3, 0, i, thread_num, f );
1884  }
1885 
1886  // スレッドを実行して，終了まで待機する
1887  do_threads_( thread, thread_num );
1888  }
1889 
1890  {
1891  for( size_type i = 0 ; i < thread_num ; i++ )
1892  {
1893  thread[ i ].setup_parameters( in, in, kernel1, kernel2, kernel3, 1, i, thread_num, f );
1894  }
1895 
1896  // スレッドを実行して，終了まで待機する
1897  do_threads_( thread, thread_num );
1898  }
1899 
1900  {
1901  for( size_type i = 0 ; i < thread_num ; i++ )
1902  {
1903  thread[ i ].setup_parameters( in, in, kernel1, kernel2, kernel3, 2, i, thread_num, f );
1904  }
1905 
1906  // スレッドを実行して，終了まで待機する
1907  do_threads_( thread, thread_num );
1908  }
1909 
1910  f( 100.1 );
1911 
1912  delete [] thread;
1913 
1914  return( true );
1915  }
1916 }
1917 
1918 
1936 template < class T1, class Allocator1, class T2, class Allocator2, class Kernel >
1937 bool linear_filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out, const Kernel &kernel, typename array< T1, Allocator1 >::size_type thread_num = 0 )
1938 {
1939  return( linear::filter( in, out, kernel, thread_num ) );
1940 }
1941 
1942 
1943 
1944 
1962 template < class T1, class Allocator1, class T2, class Allocator2, class Kernel >
1963 bool linear_filter( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out, const Kernel &kernel, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
1964 {
1965  return( linear::filter( in, out, kernel, thread_num ) );
1966 }
1967 
1968 
1986 template < class T1, class Allocator1, class T2, class Allocator2, class Kernel >
1987 bool linear_filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, const Kernel &kernel, typename array2< T1, Allocator1 >::size_type thread_num = 0 )
1988 {
1989  return( linear::filter( in, out, kernel, __mist_dmy_callback__( ), thread_num ) );
1990 }
1991 
1992 
2010 template < class T1, class Allocator1, class T2, class Allocator2, class Kernel >
2011 bool linear_filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, const Kernel &kernel, typename array3< T1, Allocator1 >::size_type thread_num = 0 )
2012 {
2013  return( linear::filter( in, out, kernel, __mist_dmy_callback__( ), thread_num ) );
2014 }
2015 
2016 
2017 
2031 
2033 namespace laplacian
2034 {
2050  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2051  bool filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out, Functor f, typename array< T1, Allocator1 >::size_type thread_num )
2052  {
2053  array< double > a( 3 );
2054  a[ 0 ] = 1.0; a[ 1 ] = -2.0; a[ 2 ] = 1.0;
2055 
2056  return( linear::filter( in, out, a, f, thread_num ) );
2057  }
2058 
2074  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2075  bool filter( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out, Functor f, typename array1< T1, Allocator1 >::size_type thread_num )
2076  {
2077  array< double > a( 3 );
2078  a[ 0 ] = 1.0; a[ 1 ] = -2.0; a[ 2 ] = 1.0;
2079 
2080  return( linear::filter( in, out, a, f, thread_num ) );
2081  }
2082 
2098  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2099  bool filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
2100  {
2101  array2< double > a( 3, 3 );
2102 
2103  a( 0, 0 ) = 1.0; a( 1, 0 ) = 1.0; a( 2, 0 ) = 1.0;
2104  a( 0, 1 ) = 1.0; a( 1, 1 ) = -8.0; a( 2, 1 ) = 1.0;
2105  a( 0, 2 ) = 1.0; a( 1, 2 ) = 1.0; a( 2, 2 ) = 1.0;
2106 
2107  return( linear::filter( in, out, a, f, thread_num ) );
2108  }
2109 
2125  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2126  bool filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
2127  {
2128  array3< double > a( 3, 3, 3 );
2129 
2130  a( 0, 0, 0 ) = 1.0; a( 1, 0, 0 ) = 1.0; a( 2, 0, 0 ) = 1.0;
2131  a( 0, 1, 0 ) = 1.0; a( 1, 1, 0 ) = 1.0; a( 2, 1, 0 ) = 1.0;
2132  a( 0, 2, 0 ) = 1.0; a( 1, 2, 0 ) = 1.0; a( 2, 2, 0 ) = 1.0;
2133  a( 0, 0, 1 ) = 1.0; a( 1, 0, 1 ) = 1.0; a( 2, 0, 1 ) = 1.0;
2134  a( 0, 1, 1 ) = 1.0; a( 1, 1, 1 ) = -26.0; a( 2, 1, 1 ) = 1.0;
2135  a( 0, 2, 1 ) = 1.0; a( 1, 2, 1 ) = 1.0; a( 2, 2, 1 ) = 1.0;
2136  a( 0, 0, 2 ) = 1.0; a( 1, 0, 2 ) = 1.0; a( 2, 0, 2 ) = 1.0;
2137  a( 0, 1, 2 ) = 1.0; a( 1, 1, 2 ) = 1.0; a( 2, 1, 2 ) = 1.0;
2138  a( 0, 2, 2 ) = 1.0; a( 1, 2, 2 ) = 1.0; a( 2, 2, 2 ) = 1.0;
2139 
2140  return( linear::filter( in, out, a, f, thread_num ) );
2141  }
2142 }
2143 
2144 
2159 template < class T1, class Allocator1, class T2, class Allocator2 >
2160 bool laplacian_filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out, typename array< T1, Allocator1 >::size_type thread_num = 0 )
2161 {
2162  return( laplacian::filter( in, out, __mist_dmy_callback__( ), thread_num ) );
2163 }
2164 
2179 template < class T1, class Allocator1, class T2, class Allocator2 >
2180 bool laplacian_filter( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
2181 {
2182  return( laplacian::filter( in, out, __mist_dmy_callback__( ), thread_num ) );
2183 }
2184 
2199 template < class T1, class Allocator1, class T2, class Allocator2 >
2200 bool laplacian_filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, typename array2< T1, Allocator1 >::size_type thread_num = 0 )
2201 {
2202  return( laplacian::filter( in, out, __mist_dmy_callback__( ), thread_num ) );
2203 }
2204 
2219 template < class T1, class Allocator1, class T2, class Allocator2 >
2220 bool laplacian_filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, typename array3< T1, Allocator1 >::size_type thread_num = 0 )
2221 {
2222  return( laplacian::filter( in, out, __mist_dmy_callback__( ), thread_num ) );
2223 }
2224 
2226 // ラプラシアングループの終わり
2227 
2228 
2229 
2230 
2231 
2244 
2246 namespace gaussian
2247 {
2248  template < bool b >
2249  struct gaussian_filter_helper
2250  {
2251  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2252  static bool filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, const double sigma, Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
2253  {
2254  array2< double > a( 3, 3, in.reso1( ), in.reso2( ) );
2255 
2256  __linear__::compute_gaussian_kernel( a, sigma );
2257 
2258  return( linear::filter( in, out, a, f, thread_num ) );
2259  }
2260 
2261  template < class T, class Allocator, class Functor >
2262  static bool filter( array2< T, Allocator > &in, const double sigma, Functor f, typename array2< T, Allocator >::size_type thread_num )
2263  {
2264  array2< double > a( 3, 3, in.reso1( ), in.reso2( ) );
2265 
2266  __linear__::compute_gaussian_kernel( a, sigma );
2267 
2268  array2< T, Allocator > tmp( in );
2269  return( linear::filter( tmp, in, a, f, thread_num ) );
2270  }
2271 
2272  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2273  static bool filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, const double sigma, Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
2274  {
2275  array3< double > a( 3, 3, 3, in.reso1( ), in.reso2( ), in.reso3( ) );
2276 
2277  __linear__::compute_gaussian_kernel( a, sigma );
2278 
2279  return( linear::filter( in, out, a, f, thread_num ) );
2280  }
2281 
2282  template < class T, class Allocator, class Functor >
2283  static bool filter( array3< T, Allocator > &in, const double sigma, Functor f, typename array3< T, Allocator >::size_type thread_num )
2284  {
2285  array3< double > a( 3, 3, 3, in.reso1( ), in.reso2( ), in.reso3( ) );
2286 
2287  __linear__::compute_gaussian_kernel( a, sigma );
2288 
2289  array2< T, Allocator > tmp( in );
2290  return( linear::filter( tmp, in, a, f, thread_num ) );
2291  }
2292  };
2293 
2294  template < >
2295  struct gaussian_filter_helper< true >
2296  {
2297  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2298  static bool filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, const double sigma, Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
2299  {
2300  // データの方が浮動小数点型の場合は高速なバージョンを利用する
2301  array1< double > a1( 3, in.reso1( ) ), a2( 3, in.reso2( ) ), a3( 3, in.reso3( ) );
2302 
2303  __linear__::compute_gaussian_kernel( a1, a2, a3, sigma );
2304 
2305  return( linear::filter1d( in, out, a1, a2, f, thread_num ) );
2306  }
2307 
2308  template < class T, class Allocator, class Functor >
2309  static bool filter( array2< T, Allocator > &in, const double sigma, Functor f, typename array2< T, Allocator >::size_type thread_num )
2310  {
2311  // データの方が浮動小数点型の場合は高速なバージョンを利用する
2312  array1< double > a1( 3, in.reso1( ) ), a2( 3, in.reso2( ) ), a3( 3, in.reso3( ) );
2313 
2314  __linear__::compute_gaussian_kernel( a1, a2, a3, sigma );
2315 
2316  return( linear::filter1d( in, a1, a2, f, thread_num ) );
2317  }
2318 
2319  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2320  static bool filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, const double sigma, Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
2321  {
2322  // データの方が浮動小数点型の場合は高速なバージョンを利用する
2323  array1< double > a1( 3, in.reso1( ) ), a2( 3, in.reso2( ) ), a3( 3, in.reso3( ) );
2324 
2325  __linear__::compute_gaussian_kernel( a1, a2, a3, sigma );
2326 
2327  return( linear::filter1d( in, out, a1, a2, a3, f, thread_num ) );
2328  }
2329 
2330  template < class T, class Allocator, class Functor >
2331  static bool filter( array3< T, Allocator > &in, const double sigma, Functor f, typename array3< T, Allocator >::size_type thread_num )
2332  {
2333  // データの方が浮動小数点型の場合は高速なバージョンを利用する
2334  array1< double > a1( 3, in.reso1( ) ), a2( 3, in.reso2( ) ), a3( 3, in.reso3( ) );
2335 
2336  __linear__::compute_gaussian_kernel( a1, a2, a3, sigma );
2337 
2338  return( linear::filter1d( in, in, a1, a2, a3, f, thread_num ) );
2339  }
2340  };
2341 
2358  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2359  bool filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out, const double sigma, Functor f, typename array< T1, Allocator1 >::size_type thread_num )
2360  {
2361  array< double > a( 3 );
2362 
2363  __linear__::compute_gaussian_kernel( a, sigma );
2364 
2365  return( linear::filter( in, out, a, f, thread_num ) );
2366  }
2367 
2384  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2385  bool filter( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out, const double sigma, Functor f, typename array1< T1, Allocator1 >::size_type thread_num )
2386  {
2387  array< double > a( 3, in.reso1( ) );
2388 
2389  __linear__::compute_gaussian_kernel( a, sigma );
2390 
2391  return( linear::filter( in, out, a, f, thread_num ) );
2392  }
2393 
2414  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2415  bool filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, const double sigma, Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
2416  {
2417  // データの方が浮動小数点型の場合は高速なバージョンを利用する
2418  return( gaussian_filter_helper< is_float< T2 >::value >::filter( in, out, sigma, f, thread_num ) );
2419  }
2420 
2439  template < class T, class Allocator, class Functor >
2440  bool filter( array2< T, Allocator > &in, const double sigma, Functor f, typename array2< T, Allocator >::size_type thread_num )
2441  {
2442  // データの方が浮動小数点型の場合は高速なバージョンを利用する
2443  return( gaussian_filter_helper< is_float< T >::value >::filter( in, sigma, f, thread_num ) );
2444  }
2445 
2466  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2467  bool filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, const double sigma, Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
2468  {
2469  // データの方が浮動小数点型の場合は高速なバージョンを利用する
2470  return( gaussian_filter_helper< is_float< T2 >::value >::filter( in, out, sigma, f, thread_num ) );
2471  }
2472 
2491  template < class T, class Allocator, class Functor >
2492  bool filter( array3< T, Allocator > &in, const double sigma, Functor f, typename array3< T, Allocator >::size_type thread_num )
2493  {
2494  // データの方が浮動小数点型の場合は高速なバージョンを利用する
2495  return( gaussian_filter_helper< is_float< T >::value >::filter( in, sigma, f, thread_num ) );
2496  }
2497 }
2498 
2499 
2515 template < class T1, class Allocator1, class T2, class Allocator2 >
2516 bool gaussian_filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out, const double sigma = 1.0, typename array< T1, Allocator1 >::size_type thread_num = 0 )
2517 {
2518  return( gaussian::filter( in, out, sigma, __mist_dmy_callback__( ), thread_num ) );
2519 }
2520 
2536 template < class T1, class Allocator1, class T2, class Allocator2 >
2537 bool gaussian_filter( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out, const double sigma = 1.0, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
2538 {
2539  return( gaussian::filter( in, out, sigma, __mist_dmy_callback__( ), thread_num ) );
2540 }
2541 
2561 template < class T1, class Allocator1, class T2, class Allocator2 >
2562 bool gaussian_filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out, const double sigma = 1.0, typename array2< T1, Allocator1 >::size_type thread_num = 0 )
2563 {
2564  return( gaussian::filter( in, out, sigma, __mist_dmy_callback__( ), thread_num ) );
2565 }
2566 
2585 template < class T, class Allocator >
2586 bool gaussian_filter( array2< T, Allocator > &in, const double sigma = 1.0, typename array2< T, Allocator >::size_type thread_num = 0 )
2587 {
2588  return( gaussian::filter( in, sigma, __mist_dmy_callback__( ), thread_num ) );
2589 }
2590 
2610 template < class T1, class Allocator1, class T2, class Allocator2 >
2611 bool gaussian_filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out, const double sigma = 1.0, typename array3< T1, Allocator1 >::size_type thread_num = 0 )
2612 {
2613  return( gaussian::filter( in, out, sigma, __mist_dmy_callback__( ), thread_num ) );
2614 }
2615 
2634 template < class T, class Allocator >
2635 bool gaussian_filter( array3< T, Allocator > &in, const double sigma = 1.0, typename array3< T, Allocator >::size_type thread_num = 0 )
2636 {
2637  return( gaussian::filter( in, sigma, __mist_dmy_callback__( ), thread_num ) );
2638 }
2639 
2641 // ガウシアングループの終わり
2642 
2643 
2644 
2657 
2659 namespace average
2660 {
2661  template < bool b >
2662  struct average_thread_helper
2663  {
2664  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2665  static bool filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2666  typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh,
2667  Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
2668  {
2669  if( is_same_object( in, out ) || in.empty( ) )
2670  {
2671  return( false );
2672  }
2673 
2674  typedef typename array2< T1, Allocator1 >::size_type size_type;
2675  typedef __linear__::average_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 >, Functor > average_thread;
2676 
2677  if( thread_num == 0 )
2678  {
2679  thread_num = static_cast< size_type >( get_cpu_num( ) );
2680  }
2681 
2682  out.resize( in.size1( ), in.size2( ) );
2683  out.reso1( in.reso1( ) );
2684  out.reso2( in.reso2( ) );
2685 
2686  average_thread *thread = new average_thread[ thread_num ];
2687 
2688  size_type i;
2689  for( i = 0 ; i < thread_num ; i++ )
2690  {
2691  thread[ i ].setup_parameters( in, out, fw, fh, 1, i, thread_num, f );
2692  }
2693 
2694  f( 0.0 );
2695 
2696  // スレッドを実行して，終了まで待機する
2697  do_threads_( thread, thread_num );
2698 
2699  f( 100.1 );
2700 
2701  delete [] thread;
2702 
2703  return( true );
2704  }
2705 
2706  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2707  static bool filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2708  typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
2709  Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
2710  {
2711  if( is_same_object( in, out ) || in.empty( ) )
2712  {
2713  return( false );
2714  }
2715 
2716  typedef typename array3< T1, Allocator1 >::size_type size_type;
2717  typedef __linear__::average_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 >, Functor > average_thread;
2718 
2719  if( thread_num == 0 )
2720  {
2721  thread_num = static_cast< size_type >( get_cpu_num( ) );
2722  }
2723 
2724  out.resize( in.size1( ), in.size2( ), in.size3( ) );
2725  out.reso1( in.reso1( ) );
2726  out.reso2( in.reso2( ) );
2727  out.reso3( in.reso3( ) );
2728 
2729  average_thread *thread = new average_thread[ thread_num ];
2730 
2731  size_type i;
2732  for( i = 0 ; i < thread_num ; i++ )
2733  {
2734  thread[ i ].setup_parameters( in, out, fw, fh, fd, i, thread_num, f );
2735  }
2736 
2737  f( 0.0 );
2738 
2739  // スレッドを実行して，終了まで待機する
2740  do_threads_( thread, thread_num );
2741 
2742  f( 100.1 );
2743 
2744  delete [] thread;
2745 
2746  return( true );
2747  }
2748  };
2749 
2750  template < >
2751  struct average_thread_helper< true >
2752  {
2753  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2754  static bool filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2755  typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh,
2756  Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
2757  {
2758  if( in.empty( ) )
2759  {
2760  return( false );
2761  }
2762 
2763  typedef typename array2< T1, Allocator1 >::size_type size_type;
2764  typedef __linear__::_1D_average_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 >, Functor > average_thread;
2765 
2766  if( thread_num == 0 )
2767  {
2768  thread_num = static_cast< size_type >( get_cpu_num( ) );
2769  }
2770 
2771  out.resize( in.size1( ), in.size2( ) );
2772  out.reso1( in.reso1( ) );
2773  out.reso2( in.reso2( ) );
2774 
2775  average_thread *thread = new average_thread[ thread_num ];
2776 
2777  size_type i;
2778  for( i = 0 ; i < thread_num ; i++ )
2779  {
2780  thread[ i ].setup_parameters( in, out, fw, fh, 1, i, thread_num, f );
2781  }
2782 
2783  f( 0.0 );
2784 
2785  // スレッドを実行して，終了まで待機する
2786  do_threads_( thread, thread_num );
2787 
2788  f( 100.1 );
2789 
2790  delete [] thread;
2791 
2792  return( true );
2793  }
2794 
2795  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2796  static bool filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2797  typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
2798  Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
2799  {
2800  if( in.empty( ) )
2801  {
2802  return( false );
2803  }
2804 
2805  typedef typename array3< T1, Allocator1 >::size_type size_type;
2806  typedef __linear__::_1D_average_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 >, Functor > average_thread;
2807 
2808  if( thread_num == 0 )
2809  {
2810  thread_num = static_cast< size_type >( get_cpu_num( ) );
2811  }
2812 
2813  out.resize( in.size1( ), in.size2( ), in.size3( ) );
2814  out.reso1( in.reso1( ) );
2815  out.reso2( in.reso2( ) );
2816  out.reso3( in.reso3( ) );
2817 
2818  average_thread *thread = new average_thread[ thread_num ];
2819 
2820  size_type i;
2821  for( i = 0 ; i < thread_num ; i++ )
2822  {
2823  thread[ i ].setup_parameters( in, out, fw, fh, fd, i, thread_num, f );
2824  }
2825 
2826  f( 0.0 );
2827 
2828  // スレッドを実行して，終了まで待機する
2829  do_threads_( thread, thread_num );
2830 
2831  f( 100.1 );
2832 
2833  delete [] thread;
2834 
2835  return( true );
2836  }
2837  };
2838 
2855  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2856  bool filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
2857  typename array< T1, Allocator1 >::size_type fw, Functor f, typename array< T1, Allocator1 >::size_type thread_num )
2858  {
2859  if( is_same_object( in, out ) || in.empty( ) )
2860  {
2861  return( false );
2862  }
2863 
2864  out.resize( in.size( ) );
2865 
2866  __linear__::__average_filter__( in, out, fw, 1, 1, 0, 1, __mist_dmy_callback__( ) );
2867 
2868  return( true );
2869  }
2870 
2887  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2888  bool filter( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out,
2889  typename array1< T1, Allocator1 >::size_type fw, Functor f, typename array1< T1, Allocator1 >::size_type thread_num )
2890  {
2891  if( is_same_object( in, out ) || in.empty( ) )
2892  {
2893  return( false );
2894  }
2895 
2896  out.resize( in.size( ) );
2897  out.reso1( in.reso1( ) );
2898 
2899  __linear__::__average_filter__( in, out, fw, 1, 1, 0, 1, __mist_dmy_callback__( ) );
2900 
2901  return( true );
2902  }
2903 
2921  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2922  bool filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
2923  typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh,
2924  Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
2925  {
2926  return( average_thread_helper< is_float< T2 >::value >::filter( in, out, fw, fh, f, thread_num ) );
2927  }
2928 
2947  template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
2948  bool filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
2949  typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
2950  Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
2951  {
2952  return( average_thread_helper< is_float< T2 >::value >::filter( in, out, fw, fh, fd, f, thread_num ) );
2953  }
2954 }
2955 
2956 
2972 template < class T1, class Allocator1, class T2, class Allocator2 >
2973 bool average_filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out, typename array< T1, Allocator1 >::size_type fw, typename array< T1, Allocator1 >::size_type thread_num = 0 )
2974 {
2975  return( average::filter( in, out, fw, __mist_dmy_callback__( ), thread_num ) );
2976 }
2977 
2993 template < class T1, class Allocator1, class T2, class Allocator2 >
2994 bool average_filter( const array1< T1, Allocator1 > &in, array1< T2, Allocator2 > &out, typename array1< T1, Allocator1 >::size_type fw, typename array1< T1, Allocator1 >::size_type thread_num = 0 )
2995 {
2996  return( average::filter( in, out, fw, __mist_dmy_callback__( ), thread_num ) );
2997 }
2998 
3015 template < class T1, class Allocator1, class T2, class Allocator2 >
3016 bool average_filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
3017  typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh, typename array2< T1, Allocator1 >::size_type thread_num = 0 )
3018 {
3019  return( average::filter( in, out, fw, fh, __mist_dmy_callback__( ), thread_num ) );
3020 }
3021 
3039 template < class T1, class Allocator1, class T2, class Allocator2 >
3040 bool average_filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
3041  typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
3042  typename array3< T1, Allocator1 >::size_type thread_num = 0 )
3043 {
3044  return( average::filter( in, out, fw, fh, fd, __mist_dmy_callback__( ), thread_num ) );
3045 }
3046 
3048 // 一様重みグループの終わり
3049 
3050 
3051 
3053 // 線形グループの終わり
3054 
3055 
3056 // mist名前空間の終わり
3057 _MIST_END
3058 
3059 
3060 
3061 #endif // __INCLUDE_FILTER_LINEAR_FILTER__
3062 
Generated on Wed Nov 12 2014 19:44:17 for MIST by doxygen 1.8.1.2