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

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

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