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