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

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

median.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 #ifndef __INCLUDE_MIST_MEDIAN_FILTER__
35 #define __INCLUDE_MIST_MEDIAN_FILTER__
36
37
38 #ifndef __INCLUDE_MIST_H__
39 #include "../mist.h"
40 #endif
41
42 #ifndef __INCLUDE_MIST_TYPE_TRAIT_H__
43 #include "../config/type_trait.h"
44 #endif
45
46 // カラー画像の設定を読み込む
47 #ifndef __INCLUDE_MIST_COLOR_H__
48 #include "../config/color.h"
49 #endif
50
51
52 #ifndef __INCLUDE_MIST_THREAD__
53 #include "../thread.h"
54 #endif
55
56 #include <algorithm>
57
58
59 // mist名前空間の始まり
60 _MIST_BEGIN
61
62 // メディアンフィルタ
63 namespace __median_filter_with_histogram__
64 {
65 // in : 入力画像. 入力画像の画素値は min と max の間とする
66 // out : 出力画像. 出力画像のメモリはあらかじめ割り当てられているものとする
67 // fw, fh, fd : フィルタサイズ
68 // min, max : 濃淡範囲
69 template < class Array1, class Array2, class Functor >
70 void median_filter( const Array1 &in, Array2 &out,
71 typename Array1::size_type fw, typename Array1::size_type fh, typename Array1::size_type fd,
72 typename Array1::value_type min, typename Array1::value_type max,
73 typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
74 typename Array1::size_type thread_idz, typename Array1::size_type thread_numz, Functor f )
75 {
76 typedef typename Array1::size_type size_type;
77 typedef typename Array1::difference_type difference_type;
78 typedef typename Array1::value_type value_type;
79 typedef typename Array2::value_type out_value_type;
80 typedef difference_type hist_value;
81
82 size_type range = static_cast< size_type >( max - min + 1 );
83 size_type a = 0, i, j, k, x, y, z, ri, leftnum;
84 difference_type pth, th, lt_med;
85 difference_type med;
86
87 size_type w = in.width( );
88 size_type h = in.height( );
89 size_type d = in.depth( );
90
91 const bool bprogress1 = thread_idy == 0 && d == 1;
92 const bool bprogress2 = thread_idz == 0 && d > 1;
93
94 size_type bw = fw / 2;
95 size_type bh = fh / 2;
96 size_type bd = fd / 2;
97
98 difference_type *leftmost = new difference_type[ fh * fd ];
99 difference_type *sort = new difference_type[ fw * fh * fd + 1 ];
100 hist_value *hist = new hist_value[ range ];
101
102 for( k = thread_idz ; k < d ; k += thread_numz )
103 {
104 for( j = thread_idy ; j < h ; j += thread_numy )
105 {
106 i = 0;
107
108 memset( hist, 0, sizeof( hist_value ) * range );
109 pth = 0;
110 leftnum = 0;
111 for( z = k < bd ? 0 : k - bd ; z <= k + bd && z < d ; z++ )
112 {
113 for( y = j < bh ? 0 : j - bh ; y <= j + bh && y < h ; y++ )
114 {
115 leftmost[ leftnum++ ] = in( 0, y, z ) - min;
116 for( x = 0 ; x <= bw ; x++ )
117 {
118 sort[ pth++ ] = in( x, y, z ) - min;
119 hist[ in( x, y, z ) - min ]++;
120 }
121 }
122 }
123
124 // 昇順に並べ替える
125 std::sort( sort, sort + pth );
126 pth--;
127
128 th = pth / 2;
129 med = sort[ th ];
130
131 if( sort[ 0 ] == sort[ th ] )
132 {
133 lt_med = 0;
134 }
135 else
136 {
137 a = th;
138 while( med <= sort[ --a ] ){}
139 lt_med = static_cast< signed int >( a + 1 );
140 }
141
142 out( i, j, k ) = static_cast< out_value_type >( med + min );
143
144 for( i = 1 ; i < w ; i++ )
145 {
146 ri = i + bw;
147
148 if( bw < i )
149 {
150 for( a = 0 ; a < leftnum ; a++ )
151 {
152 hist[ leftmost[ a ] ]--;
153 pth--;
154 if( leftmost[ a ] < med )
155 {
156 lt_med--;
157 }
158 }
159 }
160
161 leftnum = 0;
162 for( z = k < bd ? 0 : k - bd ; z <= k + bd && z < d ; z++ )
163 {
164 for( y = j < bh ? 0 : j - bh ; y <= j + bh && y < h ; y++ )
165 {
166 if( ri < w )
167 {
168 hist[ in( ri, y, z ) - min ]++;
169 pth++;
170 if( in( ri, y, z ) - min < med )
171 {
172 lt_med++;
173 }
174 }
175
176 if( bw <= i )
177 {
178 leftmost[ leftnum++ ] = in( i - bw, y, z ) - min;
179 }
180 }
181 }
182
183 th = pth / 2;
184
185 if( lt_med <= th )
186 {
187 while( lt_med + hist[ med ] <= th )
188 {
189 lt_med += hist[ med ];
190 med++;
191 }
192 }
193 else
194 {
195 while( th < lt_med )
196 {
197 med--;
198 lt_med -= hist[ med ];
199 }
200 }
201
202 out( i, j, k ) = static_cast< out_value_type >( med + min );
203 }
204
205 if( bprogress1 )
206 {
207 f( static_cast< double >( j + 1 ) / static_cast< double >( h ) * 100.0 );
208 }
209 }
210
211 if( bprogress2 )
212 {
213 f( static_cast< double >( k + 1 ) / static_cast< double >( d ) * 100.0 );
214 }
215 }
216
217 delete [] hist;
218 delete [] sort;
219 delete [] leftmost;
220 }
221 }
222
223
224 namespace __median_filter_divide_conquer__
225 {
226 // T型の配列の中央値を高速に計算するアルゴリズム
227 // 配列の要素 数はnum個である。
228 // 配列内で第c番目に小さい値を検索する
229 // 入力されるDATA型の配列は中身を入れ替えない代わりに，ワーク配列を利用する
230 template < class T >
231 inline T nth_value( const T *a, size_t num, size_t c, T *work1, T *work2, T *work3 )
232 {
233 typedef T value_type;
234
235 size_t i;
236 size_t wi1, wi2, wi3;
237 value_type s;
238 value_type *w1 = work1;
239 value_type *w2 = work2;
240 value_type *src, *www;
241
242 wi1 = wi2 = wi3 = 0;
243 s = a[ c ];
244
245 for( i = 0 ; i < num ; i++ )
246 {
247 if( a[ i ] < s )
248 {
249 w1[ wi1++ ] = a[ i ];
250 }
251 else if( a[ i ] > s )
252 {
253 w2[ wi2++ ] = a[ i ];
254 }
255 else
256 {
257 wi3++;
258 }
259 }
260
261 if( wi1 > c )
262 {
263 src = w1;
264 w1 = work3;
265 num = wi1;
266 }
267 else if( wi1 + wi3 <= c )
268 {
269 src = w2;
270 w2 = work3;
271 c -= wi1 + wi3;
272 num = wi2;
273 }
274 else
275 {
276 return ( s );
277 }
278
279 while( true )
280 {
281 wi1 = wi2 = wi3 = 0;
282 s = src[ c ];
283
284 for( i = 0 ; i < num ; i++ )
285 {
286 if( src[ i ] < s )
287 {
288 w1[ wi1++ ] = src[ i ];
289 }
290 else if( src[ i ] > s )
291 {
292 w2[ wi2++ ] = src[ i ];
293 }
294 else
295 {
296 wi3++;
297 }
298 }
299
300 if( wi1 > c )
301 {
302 www = src;
303 src = w1;
304 w1 = www;
305 num = wi1;
306 }
307 else if( wi1 + wi3 <= c )
308 {
309 www = src;
310 src = w2;
311 w2 = www;
312 c -= wi1 + wi3;
313 num = wi2;
314 }
315 else
316 {
317 break;
318 }
319 }
320
321 return ( s );
322 }
323
324
325 // in : 入力画像. 入力画像の画素値は min と max の間とする
326 // out : 出力画像. 出力画像のメモリはあらかじめ割り当てられているものとする
327 // fw, fh, fd : フィルタサイズ
328 // min, max : 濃淡範囲
329 template < class Array1, class Array2, class Functor >
330 void median_filter( const Array1 &in, Array2 &out,
331 typename Array1::size_type fw, typename Array1::size_type fh, typename Array1::size_type fd,
332 typename Array1::size_type thread_idy, typename Array1::size_type thread_numy,
333 typename Array1::size_type thread_idz, typename Array1::size_type thread_numz, Functor f )
334 {
335 typedef typename Array1::size_type size_type;
336 typedef typename Array1::value_type value_type;
337 typedef typename Array2::value_type out_value_type;
338
339 size_type i, j, k, x, y, z, ri;
340 size_type pth, c, th, windex;
341
342 size_type w = in.width( );
343 size_type h = in.height( );
344 size_type d = in.depth( );
345
346 const bool bprogress1 = thread_idy == 0 && d == 1;
347 const bool bprogress2 = thread_idz == 0 && d > 1;
348
349 size_type bw = fw / 2;
350 size_type bh = fh / 2;
351 size_type bd = fd / 2;
352
353 size_type bbh, bbd;
354
355 size_type size = fw * fh * fd;
356
357 value_type *work = new value_type[ size + 1 ];
358 value_type *work1 = new value_type[ size + 1 ];
359 value_type *work2 = new value_type[ size + 1 ];
360 value_type *work3 = new value_type[ size + 1 ];
361 value_type **sort = new value_type*[ fw ];
362
363 for( k = thread_idz ; k < d ; k += thread_numz )
364 {
365 if( k < bd )
366 {
367 bbd = k + bd + 1;
368 }
369 else if( k + bd >= d )
370 {
371 bbd = d - k + bd + 1;
372 }
373 else
374 {
375 bbd = fd;
376 }
377
378 for( j = thread_idy ; j < h ; j += thread_numy )
379 {
380 if( j < bh )
381 {
382 bbh = j + bh + 1;
383 }
384 else if( j + bh >= h )
385 {
386 bbh = h - j + bh + 1;
387 }
388 else
389 {
390 bbh = fh;
391 }
392
393 for( i = 0 ; i < fw ; i++ )
394 {
395 sort[ i ] = work + i * bbh * bbd;
396 }
397
398 i = 0;
399
400 windex = 0;
401
402 for( x = 0 ; x <= bw ; x++ )
403 {
404 windex = windex % fw;
405 c = 0;
406 for( z = k < bd ? 0 : k - bd ; z <= k + bd && z < d ; z++ )
407 {
408 for( y = j < bh ? 0 : j - bh ; y <= j + bh && y < h ; y++ )
409 {
410 sort[ windex ][ c++ ] = in( x, y, z );
411 }
412 }
413 windex++;
414 }
415
416 pth = ( bw + 1 ) * bbh * bbd;
417 th = ( pth - 1 ) / 2;
418 out( i, j, k ) = static_cast< out_value_type >( nth_value( work, pth, th, work1, work2, work3 ) );
419
420 for( i = 1 ; i < bw ; i++ )
421 {
422 ri = i + bw;
423 pth = ( i + bw + 1 ) * bbh * bbd;
424 windex = windex % fw;
425
426 c = 0;
427 for( z = k < bd ? 0 : k - bd ; z <= k + bd && z < d ; z++ )
428 {
429 for( y = j < bh ? 0 : j - bh ; y <= j + bh && y < h ; y++ )
430 {
431 sort[ windex ][ c++ ] = in( i, y, z );
432 }
433 }
434
435 th = ( pth - 1 ) / 2;
436
437 out( i, j, k ) = static_cast< out_value_type >( nth_value( work, pth, th, work1, work2, work3 ) );
438
439 windex++;
440 }
441
442 pth = fw * bbh * bbd;
443 th = ( pth - 1 ) / 2;
444 for( ; i + bw < w ; i++ )
445 {
446 ri = i + bw;
447 windex = windex % fw;
448
449 c = 0;
450 for( z = k < bd ? 0 : k - bd ; z <= k + bd && z < d ; z++ )
451 {
452 for( y = j < bh ? 0 : j - bh ; y <= j + bh && y < h ; y++ )
453 {
454 sort[ windex ][ c++ ] = in( ri, y, z );
455 }
456 }
457
458 out( i, j, k ) = static_cast< out_value_type >( nth_value( work, pth, th, work1, work2, work3 ) );
459
460 windex++;
461 }
462
463 for( ; i < w ; i++ )
464 {
465 windex = windex % fw;
466
467 c = 0;
468 for( z = k < bd ? 0 : k - bd ; z <= k + bd && z < d ; z++ )
469 {
470 for( y = j < bh ? 0 : j - bh ; y <= j + bh && y < h ; y++ )
471 {
472 sort[ windex ][ c++ ] = 0;
473 }
474 }
475
476 out( i, j, k ) = static_cast< out_value_type >( nth_value( work, pth, th, work1, work2, work3 ) );
477
478 windex++;
479 }
480
481 if( bprogress1 )
482 {
483 f( static_cast< double >( j + 1 ) / static_cast< double >( h ) * 100.0 );
484 }
485 }
486
487 if( bprogress2 )
488 {
489 f( static_cast< double >( k + 1 ) / static_cast< double >( d ) * 100.0 );
490 }
491 }
492
493 delete [] work;
494 delete [] work1;
495 delete [] work2;
496 delete [] work3;
497 delete [] sort;
498 }
499 }
500
501
502
503
504 namespace __median_filter_specialized_version__
505 {
506 template < class T >
507 inline void sort3x3( const T &v0, const T &v1, const T &v2, T v[ 3 ] )
508 {
509 if( v0 > v1 )
510 {
511 if( v1 > v2 )
512 {
513 v[ 0 ] = v0;
514 v[ 1 ] = v1;
515 v[ 2 ] = v2;
516 }
517 else
518 {
519 if( v0 > v2 )
520 {
521 v[ 0 ] = v0;
522 v[ 1 ] = v2;
523 v[ 2 ] = v1;
524 }
525 else
526 {
527 v[ 0 ] = v2;
528 v[ 1 ] = v0;
529 v[ 2 ] = v1;
530 }
531 }
532 }
533 else
534 {
535 if( v0 > v2 )
536 {
537 v[ 0 ] = v1;
538 v[ 1 ] = v0;
539 v[ 2 ] = v2;
540 }
541 else
542 {
543 if( v1 > v2 )
544 {
545 v[ 0 ] = v1;
546 v[ 1 ] = v2;
547 v[ 2 ] = v0;
548 }
549 else
550 {
551 v[ 0 ] = v2;
552 v[ 1 ] = v1;
553 v[ 2 ] = v0;
554 }
555 }
556 }
557 }
558
559 template < class T >
560 inline void sort2x2( const T &v0, const T &v1, T v[ 3 ] )
561 {
562 if( v0 > v1 )
563 {
564 v[ 0 ] = v0;
565 v[ 1 ] = v1;
566 }
567 else
568 {
569 v[ 0 ] = v1;
570 v[ 1 ] = v0;
571 }
572 }
573
574 template < class T >
575 inline void sort3x3( const T &v0, const T &v1, const T &v2, T w0[ 3 ], T w1[ 3 ], T w2[ 3 ], T *v[ 3 ] )
576 {
577 if( v0 > v1 )
578 {
579 if( v1 > v2 )
580 {
581 v[ 0 ] = w0;
582 v[ 1 ] = w1;
583 v[ 2 ] = w2;
584 }
585 else
586 {
587 if( v0 > v2 )
588 {
589 v[ 0 ] = w0;
590 v[ 1 ] = w2;
591 v[ 2 ] = w1;
592 }
593 else
594 {
595 v[ 0 ] = w2;
596 v[ 1 ] = w0;
597 v[ 2 ] = w1;
598 }
599 }
600 }
601 else
602 {
603 if( v0 > v2 )
604 {
605 v[ 0 ] = w1;
606 v[ 1 ] = w0;
607 v[ 2 ] = w2;
608 }
609 else
610 {
611 if( v1 > v2 )
612 {
613 v[ 0 ] = w1;
614 v[ 1 ] = w2;
615 v[ 2 ] = w0;
616 }
617 else
618 {
619 v[ 0 ] = w2;
620 v[ 1 ] = w1;
621 v[ 2 ] = w0;
622 }
623 }
624 }
625 }
626
627 template < class T >
628 inline void sort2x2( const T &v0, const T &v1, T w0[ 3 ], T w1[ 3 ], T *v[ 3 ] )
629 {
630 if( v0 > v1 )
631 {
632 v[ 0 ] = w0;
633 v[ 1 ] = w1;
634 }
635 else
636 {
637 v[ 0 ] = w1;
638 v[ 1 ] = w0;
639 }
640 }
641
642 template < class T >
643 inline const T &minimum( const T &v0, const T &v1 )
644 {
645 return( v0 < v1 ? v0 : v1 );
646 }
647
648 template < class T >
649 inline const T &minimum( const T &v0, const T &v1, const T &v2 )
650 {
651 return( v0 < v1 ? ( v0 < v2 ? v0 : v2 ) : ( v1 < v2 ? v1 : v2 ) );
652 }
653
654 template < class T >
655 inline const T &maximum( const T &v0, const T &v1, const T &v2 )
656 {
657 return( v0 > v1 ? ( v0 > v2 ? v0 : v2 ) : ( v1 > v2 ? v1 : v2 ) );
658 }
659
660 template < class T >
661 inline const T &maximum( const T &v0, const T &v1 )
662 {
663 return( v0 > v1 ? v0 : v1 );
664 }
665
666 template < class T >
667 inline const T &median( const T &v0, const T &v1, const T &v2 )
668 {
669 if( v0 > v1 )
670 {
671 if( v1 > v2 )
672 {
673 return( v1 );
674 }
675 else
676 {
677 if( v0 > v2 )
678 {
679 return( v2 );
680 }
681 else
682 {
683 return( v0 );
684 }
685 }
686 }
687 else
688 {
689 if( v0 > v2 )
690 {
691 return( v0 );
692 }
693 else
694 {
695 if( v1 > v2 )
696 {
697 return( v2 );
698 }
699 else
700 {
701 return( v1 );
702 }
703 }
704 }
705 }
706
707
708 /****************************************************************************************************************************************
709 **
710 ** 参考文献:
711 **
712 ** 浜村倫行, 入江文平, ``３×３メディアンフィルタの高速アルゴリズム,'' FIT（情報科学技術フォーラム）, LI-9, 2002
713 ** ``A Fast Algorithm for 3x3 Median Filtering''
714 **
715 ** ※ヒストグラムを利用するもの等と比較し，中央値を求めるのに必要な比較回数が非常に少なく，高速な動作が可能
716 **
717 ** Coded by ddeguchi.
718 **
719 ****************************************************************************************************************************************/
720 template < class T1, class T2, class Allocator1, class Allocator2, class Functor >
721 void median_filter3x3( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
722 typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num, Functor f )
723 {
724 typedef typename array2< T1, Allocator1 >::size_type size_type;
725 typedef typename array2< T1, Allocator1 >::value_type value_type;
726 typedef typename array2< T1, Allocator1 >::difference_type difference_type;
727 typedef typename array2< T2, Allocator2 >::value_type out_value_type;
728
729 size_type i, j, wi;
730 size_type w = in.width( );
731 size_type h = in.height( );
732
733 value_type work0[ 3 ];
734 value_type work1[ 3 ];
735 value_type work2[ 3 ];
736 value_type *work[ 3 ] = { work0, work1, work2 };
737 value_type *sort[ 3 ];
738
739
740 // 一番上の部分
741 for( j = thread_id ; j < 1 ; j += thread_num )
742 {
743 sort2x2( in( 0, 0 ), in( 0, 1 ), work0 );
744 sort2x2( in( 1, 0 ), in( 1, 1 ), work1 );
745
746 sort2x2( work0[ 1 ], work1[ 1 ], work0, work1, sort );
747 out( 0, 0 ) = static_cast< out_value_type >( minimum( sort[ 0 ][ 1 ], sort[ 1 ][ 0 ] ) );
748
749 for( i = 1 ; i < w - 1 ; i++ )
750 {
751 wi = ( i + 1 ) % 3;
752 sort2x2( in( i + 1, 0 ), in( i + 1, 0 + 1 ), work[ wi ] );
753 sort3x3( work0[ 1 ], work1[ 1 ], work2[ 1 ], work0, work1, work2, sort );
754
755 value_type &x = sort[ 1 ][ 1 ];
756 value_type &y = sort[ 0 ][ 1 ];
757 value_type &z = sort[ 2 ][ 0 ];
758
759 if( x < z )
760 {
761 value_type &w = sort[ 1 ][ 0 ];
762 out( i, 0 ) = static_cast< out_value_type >( minimum( y, z, w ) );
763 }
764 else
765 {
766 out( i, 0 ) = static_cast< out_value_type >( x );
767 }
768 }
769
770 sort2x2( work[ ( w - 2 ) % 3 ][ 1 ], work[ ( w - 1 ) % 3 ][ 1 ], work[ ( w - 2 ) % 3 ], work[ ( w - 1 ) % 3 ], sort );
771 out( w - 1, 0 ) = static_cast< out_value_type >( minimum( sort[ 0 ][ 1 ], sort[ 1 ][ 0 ] ) );
772
773 if( thread_id == 0 )
774 {
775 f( static_cast< double >( j + 1 ) / static_cast< double >( h ) * 100.0 );
776 }
777 }
778
779 // 真ん中部分
780 for( ; j < h - 1 ; j += thread_num )
781 {
782 sort3x3( in( 0, j - 1 ), in( 0, j ), in( 0, j + 1 ), work0 );
783 sort3x3( in( 1, j - 1 ), in( 1, j ), in( 1, j + 1 ), work1 );
784
785 sort2x2( work0[ 1 ], work1[ 1 ], work0, work1, sort );
786 out( 0, j ) = static_cast< out_value_type >( median( sort[ 0 ][ 2 ], sort[ 1 ][ 0 ], sort[ 1 ][ 1 ] ) );
787
788 for( i = 1 ; i < w - 1 ; i++ )
789 {
790 wi = ( i + 1 ) % 3;
791 sort3x3( in( i + 1, j - 1 ), in( i + 1, j ), in( i + 1, j + 1 ), work[ wi ] );
792 sort3x3( work0[ 1 ], work1[ 1 ], work2[ 1 ], work0, work1, work2, sort );
793
794 value_type &x = sort[ 1 ][ 1 ];
795 value_type &y = sort[ 0 ][ 2 ];
796 value_type &z = sort[ 2 ][ 0 ];
797
798 if( x < y && x < z )
799 {
800 value_type &w = sort[ 1 ][ 0 ];
801 out( i, j ) = static_cast< out_value_type >( minimum( y, z, w ) );
802 }
803 else if( x > y && x > z )
804 {
805 value_type &w = sort[ 1 ][ 2 ];
806 out( i, j ) = static_cast< out_value_type >( maximum( y, z, w ) );
807 }
808 else
809 {
810 out( i, j ) = static_cast< out_value_type >( x );
811 }
812 }
813
814 sort2x2( work[ ( w - 2 ) % 3 ][ 1 ], work[ ( w - 1 ) % 3 ][ 1 ], work[ ( w - 2 ) % 3 ], work[ ( w - 1 ) % 3 ], sort );
815 out( w - 1, j ) = static_cast< out_value_type >( median( sort[ 0 ][ 2 ], sort[ 1 ][ 0 ], sort[ 1 ][ 1 ] ) );
816
817 if( thread_id == 0 )
818 {
819 f( static_cast< double >( j + 1 ) / static_cast< double >( h ) * 100.0 );
820 }
821 }
822
823
824 // 一番下の部分
825 for( ; j < h ; j += thread_num )
826 {
827 sort2x2( in( 0, h - 2 ), in( 0, h - 1 ), work0 );
828 sort2x2( in( 1, h - 2 ), in( 1, h - 1 ), work1 );
829
830 sort2x2( work0[ 1 ], work1[ 1 ], work0, work1, sort );
831 out( 0, h - 1 ) = static_cast< out_value_type >( minimum( sort[ 0 ][ 1 ], sort[ 1 ][ 0 ] ) );
832
833 for( i = 1 ; i < w - 1 ; i++ )
834 {
835 wi = ( i + 1 ) % 3;
836 sort2x2( in( i + 1, h - 2 ), in( i + 1, h - 1 ), work[ wi ] );
837 sort3x3( work0[ 1 ], work1[ 1 ], work2[ 1 ], work0, work1, work2, sort );
838
839 value_type &x = sort[ 1 ][ 1 ];
840 value_type &y = sort[ 0 ][ 1 ];
841 value_type &z = sort[ 2 ][ 0 ];
842
843 if( x < z )
844 {
845 value_type &w = sort[ 1 ][ 0 ];
846 out( i, h - 1 ) = static_cast< out_value_type >( minimum( y, z, w ) );
847 }
848 else
849 {
850 out( i, h - 1 ) = static_cast< out_value_type >( x );
851 }
852 }
853
854 sort2x2( work[ ( w - 2 ) % 3 ][ 1 ], work[ ( w - 1 ) % 3 ][ 1 ], work[ ( w - 2 ) % 3 ], work[ ( w - 1 ) % 3 ], sort );
855 out( w - 1, h - 1 ) = static_cast< out_value_type >( minimum( sort[ 0 ][ 1 ], sort[ 1 ][ 0 ] ) );
856
857 if( thread_id == 0 )
858 {
859 f( static_cast< double >( j + 1 ) / static_cast< double >( h ) * 100.0 );
860 }
861 }
862 }
863 }
864
865
866 // メディアンフィルタのスレッド実装
867 namespace __median_filter_controller__
868 {
869 template < class T, class Allocator >
870 void get_min_max( const array< T, Allocator > &in, typename array< T, Allocator >::value_type &min, typename array< T, Allocator >::value_type &max )
871 {
872 min = max = in[ 0 ];
873 for( typename array< T, Allocator >::size_type i = 0 ; i < in.size( ) ; i++ )
874 {
875 if( min > in[ i ] )
876 {
877 min = in[ i ];
878 }
879 else if( max < in[ i ] )
880 {
881 max = in[ i ];
882 }
883 }
884 }
885
886 template < bool b >
887 struct __median_filter__
888 {
889 template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
890 static void median_filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
891 typename array< T1, Allocator1 >::size_type fw, typename array< T1, Allocator1 >::size_type fh, typename array< T1, Allocator1 >::size_type fd,
892 typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num, Functor f )
893 {
894 typedef typename array< T1, Allocator1 >::value_type value_type;
895 value_type min = in[ 0 ];
896 value_type max = in[ 0 ];
897 get_min_max( in, min, max );
898 __median_filter_with_histogram__::median_filter( in, out, fw, fh, fd, min, max, 0, 1, thread_id, thread_num, f );
899 }
900
901 template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
902 static void median_filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
903 typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh, typename array2< T1, Allocator1 >::size_type fd,
904 typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num, Functor f )
905 {
906 typedef typename array2< T1, Allocator1 >::value_type value_type;
907
908 if( fw == 3 && fh == 3 )
909 {
910 __median_filter_specialized_version__::median_filter3x3( in, out, thread_id, thread_num, f );
911 }
912 else
913 {
914 value_type min = in[ 0 ];
915 value_type max = in[ 0 ];
916 get_min_max( in, min, max );
917 __median_filter_with_histogram__::median_filter( in, out, fw, fh, fd, min, max, thread_id, thread_num, 0, 1, f );
918 }
919 }
920
921 template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
922 static void median_filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
923 typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
924 typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num, Functor f )
925 {
926 typedef typename array3< T1, Allocator1 >::value_type value_type;
927 value_type min = in[ 0 ];
928 value_type max = in[ 0 ];
929 get_min_max( in, min, max );
930 __median_filter_with_histogram__::median_filter( in, out, fw, fh, fd, min, max, 0, 1, thread_id, thread_num, f );
931 }
932 };
933
934 // ヒストグラムを作成できない場合のメディアンフィルタ
935 // 浮動小数点など
936 template < >
937 struct __median_filter__< false >
938 {
939 template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
940 static void median_filter( const array< T1, Allocator1 > &in, array< T2, Allocator2 > &out,
941 typename array< T1, Allocator1 >::size_type fw, typename array< T1, Allocator1 >::size_type fh, typename array< T1, Allocator1 >::size_type fd,
942 typename array< T1, Allocator1 >::size_type thread_id, typename array< T1, Allocator1 >::size_type thread_num, Functor f )
943 {
944 __median_filter_divide_conquer__::median_filter( in, out, fw, fh, fd, 0, 1, thread_id, thread_num, f );
945 }
946
947 template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
948 static void median_filter( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
949 typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh, typename array2< T1, Allocator1 >::size_type fd,
950 typename array2< T1, Allocator1 >::size_type thread_id, typename array2< T1, Allocator1 >::size_type thread_num, Functor f )
951 {
952 if( fw == 3 && fh == 3 )
953 {
954 __median_filter_specialized_version__::median_filter3x3( in, out, thread_id, thread_num, f );
955 }
956 else
957 {
958 __median_filter_divide_conquer__::median_filter( in, out, fw, fh, fd, thread_id, thread_num, 0, 1, f );
959 }
960 }
961
962 template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
963 static void median_filter( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
964 typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
965 typename array3< T1, Allocator1 >::size_type thread_id, typename array3< T1, Allocator1 >::size_type thread_num, Functor f )
966 {
967 __median_filter_divide_conquer__::median_filter( in, out, fw, fh, fd, 0, 1, thread_id, thread_num, f );
968 }
969 };
970
971
972 template < class T1, class T2, class Functor >
973 class median_thread : public mist::thread< median_thread< T1, T2, Functor > >
974 {
975 public:
976 typedef mist::thread< median_thread< T1, T2, Functor > > base;
977 typedef typename base::thread_exit_type thread_exit_type;
978 typedef typename T1::size_type size_type;
979 typedef typename T1::value_type value_type;
980
981 private:
982 size_t thread_id_;
983 size_t thread_num_;
984
985 // 入出力用の画像へのポインタ
986 const T1 *in_;
987 T2 *out_;
988 size_type fw_;
989 size_type fh_;
990 size_type fd_;
991
992 Functor f_;
993
994 public:
995 void setup_parameters( const T1 &in, T2 &out, size_t fw, size_type fh, size_type fd, size_type thread_id, size_type thread_num, Functor f )
996 {
997 in_ = &in;
998 out_ = &out;
999 fw_ = fw;
1000 fh_ = fh;
1001 fd_ = fd;
1002 thread_id_ = thread_id;
1003 thread_num_ = thread_num;
1004 f_ = f;
1005 }
1006
1007 const median_thread& operator =( const median_thread &p )
1008 {
1009 if( &p != this )
1010 {
1011 base::operator =( p );
1012 thread_id_ = p.thread_id_;
1013 thread_num_ = p.thread_num_;
1014 in_ = p.in_;
1015 out_ = p.out_;
1016 fw_ = p.fw_;
1017 fh_ = p.fh_;
1018 fd_ = p.fd_;
1019 f_ = p.f_;
1020 }
1021 return( *this );
1022 }
1023
1024 median_thread( size_type id = 0, size_type num = 1 ) : thread_id_( id ), thread_num_( num ),
1025 in_( NULL ), out_( NULL ), fw_( 3 ), fh_( 3 ), fd_( 3 )
1026 {
1027 }
1028 median_thread( const median_thread &p ) : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ),
1029 in_( p.in_ ), out_( p.out_ ), fw_( p.fw_ ), fh_( p.fh_ ), fd_( p.fd_ )
1030 {
1031 }
1032
1033 protected:
1034 // 継承した先で必ず実装されるスレッド関数
1035 virtual thread_exit_type thread_function( )
1036 {
1037 __median_filter__< type_and< is_integer< value_type >::value, !is_color< value_type >::value >::value >::median_filter( *in_, *out_, fw_, fh_, fd_, thread_id_, thread_num_, f_ );
1038 return( true );
1039 }
1040 };
1041 }
1042
1043
1051
1052
1069 template < class T1, class Allocator1, class T2, class Allocator2 >
1070 bool median( 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 )
1071 {
1072 if( is_same_object( in, out ) || in.empty( ) )
1073 {
1074 return( false );
1075 }
1076
1077 out.resize( in.size( ) );
1078
1079 typedef typename array< T1, Allocator1 >::size_type size_type;
1080 fw = static_cast< size_type >( fw / 2 ) * 2 + 1;
1081
1082 __median_filter_controller__::__median_filter__< is_integer< T1 >::value >::median_filter( in, out, fw, 1, 1, 0, 1, __mist_dmy_callback__( ) );
1083
1084 return( true );
1085 }
1086
1087
1088
1089
1106 template < class T1, class Allocator1, class T2, class Allocator2 >
1107 bool median( 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 )
1108 {
1109 if( is_same_object( in, out ) || in.empty( ) )
1110 {
1111 return( false );
1112 }
1113
1114 out.resize( in.size( ) );
1115 out.reso1( in.reso1( ) );
1116
1117 typedef typename array1< T1, Allocator1 >::size_type size_type;
1118 fw = static_cast< size_type >( fw / 2 ) * 2 + 1;
1119
1120 __median_filter_controller__::__median_filter__< is_integer< T1 >::value >::median_filter( in, out, fw, 1, 1, 0, 1, __mist_dmy_callback__( ) );
1121
1122 return( true );
1123 }
1124
1125
1126
1145 template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
1146 bool median( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
1147 typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh,
1148 Functor f, typename array2< T1, Allocator1 >::size_type thread_num )
1149 {
1150 if( is_same_object( in, out ) || in.empty( ) )
1151 {
1152 return( false );
1153 }
1154
1155 typedef typename array2< T1, Allocator1 >::size_type size_type;
1156 typedef __median_filter_controller__::median_thread< array2< T1, Allocator1 >, array2< T2, Allocator2 >, Functor > median_thread;
1157
1158 if( thread_num == 0 )
1159 {
1160 thread_num = static_cast< size_type >( get_cpu_num( ) );
1161 }
1162
1163 out.resize( in.size1( ), in.size2( ) );
1164 out.reso1( in.reso1( ) );
1165 out.reso2( in.reso2( ) );
1166
1167 median_thread *thread = new median_thread[ thread_num ];
1168
1169 fw = static_cast< size_type >( fw / 2 ) * 2 + 1;
1170 fh = static_cast< size_type >( fh / 2 ) * 2 + 1;
1171
1172 size_type i;
1173 for( i = 0 ; i < thread_num ; i++ )
1174 {
1175 thread[ i ].setup_parameters( in, out, fw, fh, 1, i, thread_num, f );
1176 }
1177
1178 f( 0.0 );
1179
1180 // スレッドを実行して，終了まで待機する
1181 do_threads_( thread, thread_num );
1182
1183 f( 100.1 );
1184
1185 delete [] thread;
1186
1187 return( true );
1188 }
1189
1190
1191
1209 template < class T1, class Allocator1, class T2, class Allocator2 >
1210 bool median( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
1211 typename array2< T1, Allocator1 >::size_type fw, typename array2< T1, Allocator1 >::size_type fh,
1212 typename array2< T1, Allocator1 >::size_type thread_num )
1213 {
1214 return( median( in, out, fw, fh, __mist_dmy_callback__( ), thread_num ) );
1215 }
1216
1217
1234 template < class T1, class Allocator1, class T2, class Allocator2 >
1235 inline bool median( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out,
1236 typename array2< T1, Allocator1 >::size_type fw,
1237 typename array2< T1, Allocator1 >::size_type thread_num = 0 )
1238 {
1239 return( median( in, out, fw, fw, __mist_dmy_callback__( ), thread_num ) );
1240 }
1241
1242
1243
1263 template < class T1, class Allocator1, class T2, class Allocator2, class Functor >
1264 bool median( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
1265 typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
1266 Functor f, typename array3< T1, Allocator1 >::size_type thread_num )
1267 {
1268 if( is_same_object( in, out ) || in.empty( ) )
1269 {
1270 return( false );
1271 }
1272
1273 typedef typename array3< T1, Allocator1 >::size_type size_type;
1274 typedef __median_filter_controller__::median_thread< array3< T1, Allocator1 >, array3< T2, Allocator2 >, Functor > median_thread;
1275
1276 if( thread_num == 0 )
1277 {
1278 thread_num = static_cast< size_type >( get_cpu_num( ) );
1279 }
1280
1281 out.resize( in.size1( ), in.size2( ), in.size3( ) );
1282 out.reso1( in.reso1( ) );
1283 out.reso2( in.reso2( ) );
1284 out.reso3( in.reso3( ) );
1285
1286 median_thread *thread = new median_thread[ thread_num ];
1287
1288 fw = static_cast< size_type >( fw / 2 ) * 2 + 1;
1289 fh = static_cast< size_type >( fh / 2 ) * 2 + 1;
1290 fd = static_cast< size_type >( fd / 2 ) * 2 + 1;
1291
1292 size_type i;
1293 for( i = 0 ; i < thread_num ; i++ )
1294 {
1295 thread[ i ].setup_parameters( in, out, fw, fh, fd, i, thread_num, f );
1296 }
1297
1298 f( 0.0 );
1299
1300 // スレッドを実行して，終了まで待機する
1301 do_threads_( thread, thread_num );
1302
1303 f( 100.1 );
1304
1305 delete [] thread;
1306
1307 return( true );
1308 }
1309
1310
1311
1330 template < class T1, class Allocator1, class T2, class Allocator2 >
1331 bool median( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
1332 typename array3< T1, Allocator1 >::size_type fw, typename array3< T1, Allocator1 >::size_type fh, typename array3< T1, Allocator1 >::size_type fd,
1333 typename array3< T1, Allocator1 >::size_type thread_num )
1334 {
1335 return( median( in, out, fw, fh, fd, __mist_dmy_callback__( ), thread_num ) );
1336 }
1337
1338
1355 template < class T1, class Allocator1, class T2, class Allocator2 >
1356 inline bool median( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out,
1357 typename array3< T1, Allocator1 >::size_type fw,
1358 typename array3< T1, Allocator1 >::size_type thread_num = 0 )
1359 {
1360 return( median( in, out, fw, fw, fw, __mist_dmy_callback__( ), thread_num ) );
1361 }
1362
1363
1365 // メディアングループの終わり
1366
1367
1368 // mist名前空間の終わり
1369 _MIST_END
1370
1371
1372 #endif // __INCLUDE_MIST_MEDIAN_FILTER__
1373