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

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

distance.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
42
43 #ifndef __INCLUDE_MIST_DISTANCE_TRANSFORM__
44 #define __INCLUDE_MIST_DISTANCE_TRANSFORM__
45
46 #include <cmath>
47
48 #ifndef __INCLUDE_MIST_H__
49 #include "../mist.h"
50 #endif
51
52 #ifndef __INCLUDE_MIST_THREAD__
53 #include "../thread.h"
54 #endif
55
56 #ifndef __INCLUDE_MIST_LIMITS__
57 #include "../limits.h"
58 #endif
59
60
61
62 // mist名前空間の始まり
63 _MIST_BEGIN
64
65
66 namespace _distance_utility_
67 {
68 template < int DIMENSION >
69 struct __access__
70 {
71 template < class Array >
72 inline static typename Array::value_type &at( Array &in, typename Array::size_type _1, typename Array::size_type _2, typename Array::size_type _3 )
73 {
74 return( in( _1, _2, _3 ) );
75 }
76
77 template < class Array >
78 inline static typename Array::size_type size1( const Array &in ){ return( in.size1( ) ); }
79
80 template < class Array >
81 inline static typename Array::size_type size2( const Array &in ){ return( in.size2( ) ); }
82
83 template < class Array >
84 inline static typename Array::size_type size3( const Array &in ){ return( in.size3( ) ); }
85
86 template < class Array >
87 inline static double aspect( const Array &in ){ return( 1.0 ); }
88 };
89
90 template < >
91 struct __access__< 2 >
92 {
93 template < class Array >
94 inline static typename Array::value_type &at( Array &in, typename Array::size_type _1, typename Array::size_type _2, typename Array::size_type _3 )
95 {
96 return( in( _2, _1, _3 ) );
97 }
98
99 template < class Array >
100 inline static typename Array::size_type size1( const Array &in ){ return( in.size2( ) ); }
101
102 template < class Array >
103 inline static typename Array::size_type size2( const Array &in ){ return( in.size1( ) ); }
104
105 template < class Array >
106 inline static typename Array::size_type size3( const Array &in ){ return( in.size3( ) ); }
107
108 template < class Array >
109 inline static double aspect( const Array &in ){ return( in.reso2( ) / in.reso1( ) ); }
110 };
111
112 template < >
113 struct __access__< 3 >
114 {
115 template < class Array >
116 inline static typename Array::value_type &at( Array &in, typename Array::size_type _1, typename Array::size_type _2, typename Array::size_type _3 )
117 {
118 return( in( _2, _3, _1 ) );
119 }
120
121 template < class Array >
122 inline static typename Array::size_type size1( const Array &in ){ return( in.size3( ) ); }
123
124 template < class Array >
125 inline static typename Array::size_type size2( const Array &in ){ return( in.size1( ) ); }
126
127 template < class Array >
128 inline static typename Array::size_type size3( const Array &in ){ return( in.size2( ) ); }
129
130 template < class Array >
131 inline static double aspect( const Array &in ){ return( in.reso3( ) / in.reso1( ) ); }
132 };
133
134 template < int DIMENSION >
135 struct __range__
136 {
137 typedef size_t size_type;
138
139 size_type sx;
140 size_type ex;
141 size_type sy;
142 size_type ey;
143 size_type sz;
144 size_type ez;
145
146 __range__( ) : sx( 0 ), ex( 0 ), sy( 0 ), ey( 0 ), sz( 0 ), ez( 0 )
147 {
148 }
149
150 __range__( size_type ssx, size_type eex, size_type ssy, size_type eey, size_type ssz, size_type eez )
151 : sx( ssx ), ex( eex ), sy( ssy ), ey( eey ), sz( ssz ), ez( eez )
152 {
153 }
154
155 template < int DIM >
156 __range__( const __range__< DIM > &r )
157 : sx( r.sx ), ex( r.ex ), sy( r.sy ), ey( r.ey ), sz( r.sz ), ez( r.ez )
158 {
159 }
160
161 size_type begin1( ) const { return( sx ); }
162 size_type end1( ) const { return( ex ); }
163
164 size_type begin2( ) const { return( sy ); }
165 size_type end2( ) const { return( ey ); }
166
167 size_type begin3( ) const { return( sz ); }
168 size_type end3( ) const { return( ez ); }
169 };
170
171 template < >
172 struct __range__< 2 >
173 {
174 typedef size_t size_type;
175
176 size_type sx;
177 size_type ex;
178 size_type sy;
179 size_type ey;
180 size_type sz;
181 size_type ez;
182
183 __range__( ) : sx( 0 ), ex( 0 ), sy( 0 ), ey( 0 ), sz( 0 ), ez( 0 )
184 {
185 }
186
187 __range__( size_type ssx, size_type eex, size_type ssy, size_type eey, size_type ssz, size_type eez )
188 : sx( ssx ), ex( eex ), sy( ssy ), ey( eey ), sz( ssz ), ez( eez )
189 {
190 }
191
192 template < int DIM >
193 __range__( const __range__< DIM > &r )
194 : sx( r.sx ), ex( r.ex ), sy( r.sy ), ey( r.ey ), sz( r.sz ), ez( r.ez )
195 {
196 }
197
198 size_type begin1( ) const { return( sy ); }
199 size_type end1( ) const { return( ey ); }
200
201 size_type begin2( ) const { return( sx ); }
202 size_type end2( ) const { return( ex ); }
203
204 size_type begin3( ) const { return( sz ); }
205 size_type end3( ) const { return( ez ); }
206 };
207
208 template < >
209 struct __range__< 3 >
210 {
211 typedef size_t size_type;
212
213 size_type sx;
214 size_type ex;
215 size_type sy;
216 size_type ey;
217 size_type sz;
218 size_type ez;
219
220 __range__( ) : sx( 0 ), ex( 0 ), sy( 0 ), ey( 0 ), sz( 0 ), ez( 0 )
221 {
222 }
223
224 __range__( size_type ssx, size_type eex, size_type ssy, size_type eey, size_type ssz, size_type eez )
225 : sx( ssx ), ex( eex ), sy( ssy ), ey( eey ), sz( ssz ), ez( eez )
226 {
227 }
228
229 template < int DIM >
230 __range__( const __range__< DIM > &r )
231 : sx( r.sx ), ex( r.ex ), sy( r.sy ), ey( r.ey ), sz( r.sz ), ez( r.ez )
232 {
233 }
234
235 size_type begin1( ) const { return( sz ); }
236 size_type end1( ) const { return( ez ); }
237
238 size_type begin2( ) const { return( sx ); }
239 size_type end2( ) const { return( ex ); }
240
241 size_type begin3( ) const { return( sy ); }
242 size_type end3( ) const { return( ey ); }
243 };
244 }
245
246
247 // 斉藤先生によるユークリッド型距離変換の実装部分
248 namespace __saito__
249 {
250 template < class T >
251 void euclidean_distance_transform_x( T &in, double max_length = -1.0, typename T::size_type thread_id = 0, typename T::size_type thread_num = 1 )
252 {
253 typedef typename T::difference_type difference_type;
254 typedef typename T::value_type value_type;
255
256 difference_type i, j, k;
257 value_type nd;
258
259 const difference_type w = in.width( );
260 const difference_type h = in.height( );
261 const difference_type d = in.depth( );
262
263 value_type max = static_cast< value_type >( max_length <= 0 ? std::sqrt( static_cast< double >( type_limits< value_type >::maximum( ) ) ) : max_length );
264
265 for( k = 0 ; k < d ; k++ )
266 {
267 for( j = thread_id ; j < h ; j += thread_num )
268 {
269 if( in( 0, j, k ) != 0 )
270 {
271 nd = static_cast< value_type >( w ) < max ? static_cast< value_type >( w ) : max;
272 in( 0, j, k ) = nd * nd;
273 }
274 else
275 {
276 nd = 0;
277 }
278
279 for( i = 1 ; i < w ; i++ )
280 {
281 if( in( i, j, k ) != 0 )
282 {
283 nd = nd + 1 < max ? nd + 1 : nd;
284 in( i, j, k ) = nd * nd;
285 }
286 else
287 {
288 nd = 0;
289 }
290 }
291
292 if( in( w - 1, j, k ) != 0 )
293 {
294 nd = static_cast< value_type >( w ) < max ? static_cast< value_type >( w ) : max;
295 in( w - 1, j, k ) = in( w - 1, j, k ) < nd * nd ? in( w - 1, j, k ) : nd * nd;
296 }
297 else
298 {
299 nd = 0;
300 }
301
302 for( i = w - 2 ; i >= 0 ; i-- )
303 {
304 if( in( i, j, k ) != 0 )
305 {
306 nd = nd + 1 < max ? nd + 1 : nd;
307 in( i, j, k ) = in( i, j, k ) < nd * nd ? in( i, j, k ) : nd * nd;
308 }
309 else
310 {
311 nd = 0;
312 }
313 }
314 }
315 }
316 }
317
318 template < class T >
319 void euclidean_distance_transform_y( T &in, double max_length = -1.0, typename T::size_type thread_id = 0, typename T::size_type thread_num = 1 )
320 {
321 typedef typename T::size_type size_type;
322 typedef typename T::value_type value_type;
323 typedef typename T::difference_type difference_type;
324
325 size_type i, j, k;
326 difference_type n;
327
328 double wmin, wtmp;
329 size_type irange;
330 difference_type irange1, irange2;
331 double *work;
332 double vy, vyvy;
333
334 size_type w = in.width( );
335 size_type h = in.height( );
336 size_type d = in.depth( );
337
338 value_type max = static_cast< value_type >( max_length <= 0 ? std::sqrt( static_cast< double >( type_limits< value_type >::maximum( ) ) ) : max_length );
339
340 vy = in.reso2( ) / in.reso1( );
341 vyvy = vy * vy;
342 work = new double[ h ];
343
344 for( k = 0 ; k < d ; k++ )
345 {
346 for( i = thread_id ; i < w ; i += thread_num )
347 {
348 for( j = 0 ; j < h ; j++ )
349 {
350 work[ j ] = static_cast< double >( in( i, j, k ) );
351 }
352
353 for( j = 0 ; j < h ; j++ )
354 {
355 wmin = work[ j ];
356
357 if( wmin != 0.0 )
358 {
359 irange = static_cast< size_type >( std::sqrt( wmin ) / vy ) + 1;
360 irange = static_cast< size_type >( static_cast< value_type >( irange ) < max ? irange : max );
361
362 irange1 = j < irange ? j : irange;
363 irange2 = j + irange >= h ? h - 1 - j : irange;
364 for( n = -irange1 ; n <= irange2 ; n++ )
365 {
366 wtmp = work[ j + n ] + static_cast< double >( n * n * vyvy );
367 wmin = wmin > wtmp ? wtmp : wmin;
368 }
369
370 in( i, j, k ) = static_cast< value_type >( wmin );
371 }
372 }
373 }
374 }
375
376 delete [] work;
377 }
378
379
380 template < class T >
381 void euclidean_distance_transform_z( T &in, double max_length = -1.0, typename T::size_type thread_id = 0, typename T::size_type thread_num = 1 )
382 {
383 typedef typename T::size_type size_type;
384 typedef typename T::value_type value_type;
385 typedef typename T::difference_type difference_type;
386
387 size_type i, j, k;
388 difference_type n;
389 double wmin, wtmp;
390 size_type irange;
391 difference_type irange1, irange2;
392 double *work;
393 double vz, vzvz;
394
395 size_type w = in.width( );
396 size_type h = in.height( );
397 size_type d = in.depth( );
398
399 vz = in.reso3( ) / in.reso1( );
400 vzvz = vz * vz;
401 work = new double[ d ];
402
403 value_type max = static_cast< value_type >( max_length <= 0 ? std::sqrt( static_cast< double >( type_limits< value_type >::maximum( ) ) ) : max_length );
404
405 for( j = 0 ; j < h ; j++ )
406 {
407 for( i = thread_id ; i < w ; i += thread_num )
408 {
409 for( k = 0 ; k < d ; k++ )
410 {
411 work[ k ] = static_cast< double >( in( i, j, k ) );
412 }
413
414 for( k = 0 ; k < d ; k++ )
415 {
416 wmin = work[ k ];
417
418 if( wmin != 0.0 )
419 {
420 irange = static_cast< size_type >( std::sqrt( wmin ) / vz ) + 1;
421 irange = static_cast< size_type >( static_cast< value_type >( irange ) < max ? irange : max );
422
423 irange1 = k < irange ? k : irange;
424 irange2 = k + irange >= d ? d - 1 - k : irange;
425 for( n = -irange1 ; n <= irange2 ; n++ )
426 {
427 wtmp = work[ k + n ] + static_cast< double >( n * n * vzvz );
428 wmin = wmin > wtmp ? wtmp : wmin;
429 }
430
431 in( i, j, k ) = static_cast< value_type >( wmin );
432 }
433 }
434 }
435 }
436
437 delete [] work;
438 }
439
440
442 template < int DIMENSION >
443 struct __inverse_distance_transform__
444 {
445 template < class Array >
446 static void distance_transform( Array &in, const _distance_utility_::__range__< DIMENSION > &range, typename Array::size_type thread_id = 0, typename Array::size_type thread_num = 1 )
447 {
448 typedef typename Array::size_type size_type;
449 typedef typename Array::value_type value_type;
450 typedef typename Array::pointer pointer;
451 typedef typename Array::difference_type difference_type;
452 typedef _distance_utility_::__access__< DIMENSION > access;
453
454 const value_type max = type_limits< value_type >::maximum( );
455 double len = 0.0;
456
457 size_type _1s = range.begin1( );
458 size_type _2s = range.begin2( );
459 size_type _3s = range.begin3( );
460 size_type _1e = range.end1( );
461 size_type _2e = range.end2( );
462 size_type _3e = range.end3( );
463 difference_type inum = _1e - _1s + 1;
464
465 value_type *g = new value_type[ inum + 1 ];
466
467 double as = access::aspect( in );
468 double as2 = as * as;
469 double _2as = 1.0 / ( 2.0 * as2 );
470 difference_type diff = &access::at( in, 1, 0, 0 ) - &access::at( in, 0, 0, 0 );
471
472 for( size_type i3 = _3s + thread_id ; i3 <= _3e ; i3 += thread_num )
473 {
474 for( size_type i2 = _2s ; i2 <= _2e ; i2++ )
475 {
476 pointer sp = &access::at( in, _1s, i2, i3 );
477 pointer ep = &access::at( in, _1e, i2, i3 );
478
479 memset( g, 0, sizeof( value_type ) * ( _1e - _1s + 1 + 1 ) );
480
481 {
482 pointer op = sp;
483 pointer p = sp + diff;
484 g[ 0 ] = *sp;
485 for( difference_type i1 = 1 ; i1 < inum ; i1++, p += diff )
486 {
487 if( *p < *op )
488 {
489 difference_type num = static_cast< difference_type >( ( *op - *p - 1 ) * _2as );
490 pointer pp = p;
491 for( difference_type n = 0 ; n <= num && pp <= ep ; n++, pp += diff )
492 {
493 double nn = *op - ( n + 1 ) * ( n + 1 ) * as2;
494 if( *pp >= nn )
495 {
496 break;
497 }
498 else if( g[ i1 + n ] < nn )
499 {
500 g[ i1 + n ] = static_cast< value_type >( nn );
501 }
502 }
503 }
504 else if( g[ i1 ] < *p )
505 {
506 g[ i1 ] = *p;
507 }
508
509 op = p;
510 }
511 }
512
513 {
514 pointer p = ep;
515 for( difference_type i1 = inum - 1 ; i1 >= 0 ; i1--, p -= diff )
516 {
517 if( g[ i1 ] < g[ i1 + 1 ] )
518 {
519 difference_type num = static_cast< difference_type >( ( g[ i1 + 1 ] - g[ i1 ] - 1 ) * _2as );
520 pointer pp = p;
521 for( difference_type n = 0 ; n <= num && pp >= sp ; n++, pp -= diff )
522 {
523 double nn = g[ i1 + 1 ] - ( n + 1 ) * ( n + 1 ) * as2;
524 if( g[ i1 - n ] >= nn )
525 {
526 break;
527 }
528 else if( *pp < nn )
529 {
530 *pp = static_cast< value_type >( nn );
531 }
532 }
533 }
534 else if( *p < g[ i1 ] )
535 {
536 *p = g[ i1 ];
537 }
538 }
539 }
540 }
541 }
542
543 delete [] g;
544 }
545 };
546
547
548 // 斉藤先生によるユークリッド型距離変換のスレッド部分
549 template < class T >
550 class saito_distance_transform_thread : public mist::thread< saito_distance_transform_thread< T > >
551 {
552 public:
553 typedef mist::thread< saito_distance_transform_thread< T > > base;
554 typedef typename base::thread_exit_type thread_exit_type;
555 typedef typename T::size_type size_type;
556 typedef typename T::value_type value_type;
557
558 private:
559 size_t thread_id_;
560 size_t thread_num_;
561
562 // 入出力用の画像へのポインタ
563 T *in_;
564 double max_length_;
565 size_type axis_;
566
567 public:
568 void setup_parameters( T &in, double max_length, size_type axis, size_type thread_id, size_type thread_num )
569 {
570 in_ = &in;
571 max_length_ = max_length;
572 axis_ = axis;
573 thread_id_ = thread_id;
574 thread_num_ = thread_num;
575 }
576
577 void setup_axis( size_type axis )
578 {
579 axis_ = axis;
580 }
581
582 const saito_distance_transform_thread& operator =( const saito_distance_transform_thread &p )
583 {
584 if( &p != this )
585 {
586 base::operator =( p );
587 thread_id_ = p.thread_id_;
588 thread_num_ = p.thread_num_;
589 in_ = p.in_;
590 max_length_ = p.max_length;
591 axis_ = p.axis_;
592 }
593 return( *this );
594 }
595
596 saito_distance_transform_thread( size_type id = 0, size_type num = 1 )
597 : thread_id_( id ), thread_num_( num ), in_( NULL ), max_length_( -1.0 ), axis_( 0 )
598 {
599 }
600 saito_distance_transform_thread( const saito_distance_transform_thread &p )
601 : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ), in_( NULL ), max_length_( -1.0 ), axis_( p.axis_ )
602 {
603 }
604
605 protected:
606 // 継承した先で必ず実装されるスレッド関数
607 virtual thread_exit_type thread_function( )
608 {
609 switch( axis_ )
610 {
611 case 1:
612 __saito__::euclidean_distance_transform_y( *in_, max_length_, thread_id_, thread_num_ );
613 break;
614
615 case 2:
616 __saito__::euclidean_distance_transform_z( *in_, max_length_, thread_id_, thread_num_ );
617 break;
618
619 case 0:
620 default:
621 __saito__::euclidean_distance_transform_x( *in_, max_length_, thread_id_, thread_num_ );
622 break;
623 }
624 return( true );
625 }
626 };
627
628 template < class T >
629 class saito_inverse_distance_transform_thread : public mist::thread< saito_inverse_distance_transform_thread< T > >
630 {
631 public:
632 typedef mist::thread< saito_inverse_distance_transform_thread< T > > base;
633 typedef typename base::thread_exit_type thread_exit_type;
634 typedef typename T::size_type size_type;
635 typedef typename T::value_type value_type;
636
637 protected:
638 size_t thread_id_;
639 size_t thread_num_;
640
641 // 入出力用の画像へのポインタ
642 T *in_;
643 size_type axis_;
644 _distance_utility_::__range__< 1 > range_;
645
646 public:
647 void setup_parameters( T &in, const _distance_utility_::__range__< 1 > &range, size_type axis, size_type thread_id, size_type thread_num )
648 {
649 in_ = &in;
650 range_ = range;
651 axis_ = axis;
652 thread_id_ = thread_id;
653 thread_num_ = thread_num;
654 }
655
656 void setup_axis( size_type axis )
657 {
658 axis_ = axis;
659 }
660
661 const saito_inverse_distance_transform_thread& operator =( const saito_inverse_distance_transform_thread &p )
662 {
663 if( &p != this )
664 {
665 base::operator =( p );
666 thread_id_ = p.thread_id_;
667 thread_num_ = p.thread_num_;
668 in_ = p.in_;
669 range_ = p.range_;
670 axis_ = p.axis_;
671 }
672 return( *this );
673 }
674
675 saito_inverse_distance_transform_thread( size_type id = 0, size_type num = 1 )
676 : thread_id_( id ), thread_num_( num ), in_( NULL ), axis_( 0 )
677 {
678 }
679 saito_inverse_distance_transform_thread( const saito_inverse_distance_transform_thread &p )
680 : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ), in_( NULL ), axis_( p.axis_ )
681 {
682 }
683
684 protected:
685 // 継承した先で必ず実装されるスレッド関数
686 virtual thread_exit_type thread_function( )
687 {
688 switch( axis_ )
689 {
690 case 0:
691 __inverse_distance_transform__< 1 >::distance_transform( *in_, _distance_utility_::__range__< 1 >( range_ ), thread_id_, thread_num_ );
692 break;
693
694 case 1:
695 __inverse_distance_transform__< 2 >::distance_transform( *in_, _distance_utility_::__range__< 2 >( range_ ), thread_id_, thread_num_ );
696 break;
697
698 case 2:
699 default:
700 __inverse_distance_transform__< 3 >::distance_transform( *in_, _distance_utility_::__range__< 3 >( range_ ), thread_id_, thread_num_ );
701 break;
702 }
703 return( true );
704 }
705 };
706 }
707
708
709
710 // Calvinによるユークリッド型距離変換の実装部分
711 namespace __calvin__
712 {
713 inline bool remove_edt( const double uR, const double vR, const double wR, const double ud, const double vd, const double wd )
714 {
715 const double a = vd - ud;
716 const double b = wd - vd;
717 const double c = a + b;
718 return( c * vR - b * uR - a * wR - a * b * c > 0 );
719 }
720
721
723 template < int DIMENSION >
724 struct __distance_transform__
725 {
726 template < class Array >
727 static void distance_transform( Array &in, const _distance_utility_::__range__< DIMENSION > &range, typename Array::size_type thread_id = 0, typename Array::size_type thread_num = 1 )
728 {
729 typedef typename Array::size_type size_type;
730 typedef typename Array::value_type value_type;
731 typedef typename Array::pointer pointer;
732 typedef typename Array::difference_type difference_type;
733 typedef _distance_utility_::__access__< DIMENSION > access;
734
735 size_type _1s = range.begin1( );
736 size_type _2s = range.begin2( );
737 size_type _3s = range.begin3( );
738 size_type _1e = range.end1( );
739 size_type _2e = range.end2( );
740 size_type _3e = range.end3( );
741
742 double as = access::aspect( in );
743 difference_type diff = &access::at( in, 1, 0, 0 ) - &access::at( in, 0, 0, 0 );
744
745 double *g = new double[ _1e - _1s + 2 ];
746 double *h = new double[ _1e - _1s + 2 ];
747
748 for( size_type i3 = _3s + thread_id ; i3 <= _3e ; i3 += thread_num )
749 {
750 for( size_type i2 = _2s ; i2 <= _2e ; i2++ )
751 {
752 pointer p = &access::at( in, 0, i2, i3 );
753 for( size_type _i1s = _1s ; _i1s <= _1e ; )
754 {
755 // 最初に画素が始まる位置を検索する
756 for( ; _i1s <= _1e ; _i1s++ )
757 {
758 if( p[ _i1s * diff ] != 0 )
759 {
760 break;
761 }
762 }
763
764 if( _i1s > _1e )
765 {
766 // 処理対象が存在しないのでスキップ
767 break;
768 }
769
770 // 連続する画素の終端を検索する
771 size_type _i1e;
772 for( _i1e = _i1s ; _i1e <= _1e ; _i1e++ )
773 {
774 if( p[ _i1e * diff ] == 0 )
775 {
776 break;
777 }
778 }
779
780 if( _i1s > 0 )
781 {
782 _i1s--;
783 }
784 if( _i1e > _1e )
785 {
786 _i1e = _1e;
787 }
788
789 difference_type l = 0;
790
791 for( size_type i1 = _i1s ; i1 <= _i1e ; i1++ )
792 {
793 double nd = static_cast< double >( i1 - _i1s + 1 ) * as;
794 double fn = static_cast< double >( p[ i1 * diff ] );
795
796 if( l < 2 )
797 {
798 l++;
799 g[ l ] = fn;
800 h[ l ] = nd;
801 }
802 else
803 {
804 while( l >= 2 && remove_edt( g[ l - 1 ], g[ l ], fn, h[ l - 1 ], h[ l ], nd ) )
805 {
806 l--;
807 }
808 l++;
809 g[ l ] = fn;
810 h[ l ] = nd;
811 }
812 }
813
814 if( l == 0 )
815 {
816 break;
817 }
818
819 difference_type ns = l;
820 l = 1;
821
822 for( size_type i1 = _i1s ; i1 <= _i1e ; i1++ )
823 {
824 double nd = ( i1 - _i1s + 1 ) * as;
825 double len = h[ l ] - nd;
826
827 len = g[ l ] + len * len;
828 for( ; l < ns ; l++ )
829 {
830 double len_ = h[ l + 1 ] - nd;
831 len_ = g[ l + 1 ] + len_ * len_;
832 if( len > len_ )
833 {
834 len = len_;
835 }
836 else
837 {
838 break;
839 }
840 }
841 p[ i1 * diff ] = static_cast< value_type >( len );
842 }
843
844 // 次の位置に設定する
845 _i1s = _i1e + 1;
846 }
847 }
848 }
849
850 delete [] g;
851 delete [] h;
852 }
853 };
854
855
857 template < >
858 struct __distance_transform__< 1 >
859 {
860 template < class Array >
861 static void distance_transform( Array &in, const _distance_utility_::__range__< 1 > &range, typename Array::size_type thread_id = 0, typename Array::size_type thread_num = 1 )
862 {
863 typedef typename Array::size_type size_type;
864 typedef typename Array::value_type value_type;
865 typedef typename Array::pointer pointer;
866 typedef typename Array::difference_type difference_type;
867
868 const difference_type w = in.width( );
869 value_type max = type_limits< value_type >::maximum( );
870 value_type len;
871
872 max = static_cast< value_type >( w ) < max ? static_cast< value_type >( w ) : max;
873
874 size_type sx = range.begin1( );
875 size_type sy = range.begin2( );
876 size_type sz = range.begin3( );
877 size_type ex = range.end1( );
878 size_type ey = range.end2( );
879 size_type ez = range.end3( );
880
881 for( size_type k = sz ; k <= ez ; k++ )
882 {
883 for( size_type j = sy + thread_id ; j <= ey ; j += thread_num )
884 {
885 pointer sp = &in( sx, j, k );
886 pointer ep = &in( ex, j, k );
887
888 if( sp[ 0 ] != 0 )
889 {
890 len = max;
891 sp[ 0 ] = len * len;
892 }
893 else
894 {
895 len = 0;
896 }
897
898 for( pointer p = sp + 1 ; p <= ep ; p++ )
899 {
900 if( *p != 0 )
901 {
902 len = len + 1 < max ? len + 1 : len;
903 *p = len * len;
904 }
905 else
906 {
907 len = 0;
908 }
909 }
910
911 if( ep[ 0 ] != 0 )
912 {
913 len = max;
914 ep[ 0 ] = ep[ 0 ] < len * len ? ep[ 0 ] : len * len;
915 }
916 else
917 {
918 len = 0;
919 }
920
921 for( pointer p = ep - 1 ; p >= sp ; p-- )
922 {
923 if( *p != 0 )
924 {
925 len = len + 1 < max ? len + 1 : len;
926 *p = *p < len * len ? *p : len * len;
927 }
928 else
929 {
930 len = 0;
931 }
932 }
933 }
934 }
935 }
936 };
937
938
940 template < int DIMENSION >
941 struct __voronoi_distance_transform__
942 {
943 template < class Array1, class Array2 >
944 static void distance_transform( Array1 &voronoi, Array2 &dist, const _distance_utility_::__range__< DIMENSION > &range, typename Array2::size_type thread_id = 0, typename Array2::size_type thread_num = 1 )
945 {
946 typedef typename Array2::size_type size_type;
947 typedef typename Array2::value_type value_type;
948 typedef typename Array1::pointer ipointer;
949 typedef typename Array2::pointer pointer;
950 typedef typename Array2::difference_type difference_type;
951 typedef _distance_utility_::__access__< DIMENSION > access;
952
953 size_type _1s = range.begin1( );
954 size_type _2s = range.begin2( );
955 size_type _3s = range.begin3( );
956 size_type _1e = range.end1( );
957 size_type _2e = range.end2( );
958 size_type _3e = range.end3( );
959
960 double as = access::aspect( dist );
961 difference_type diff = &access::at( dist, 1, 0, 0 ) - &access::at( dist, 0, 0, 0 );
962
963 difference_type *idx = new difference_type[ _1e - _1s + 2 ];
964 double *g = new double[ _1e - _1s + 2 ];
965 double *h = new double[ _1e - _1s + 2 ];
966
967 for( size_type i3 = _3s + thread_id ; i3 <= _3e ; i3 += thread_num )
968 {
969 for( size_type i2 = _2s ; i2 <= _2e ; i2++ )
970 {
971 ipointer vp = &access::at( voronoi, 0, i2, i3 );
972 pointer p = &access::at( dist, 0, i2, i3 );
973
974 difference_type l = 0;
975
976 for( size_type i1 = _1s ; i1 <= _1e ; i1++ )
977 {
978 double nd = static_cast< double >( i1 - _1s + 1 ) * as;
979 double fn = static_cast< double >( p[ i1 * diff ] );
980
981 if( l < 2 )
982 {
983 l++;
984 g[ l ] = fn;
985 h[ l ] = nd;
986 idx[ l ] = i1;
987 }
988 else
989 {
990 while( l >= 2 && remove_edt( g[ l - 1 ], g[ l ], fn, h[ l - 1 ], h[ l ], nd ) )
991 {
992 l--;
993 }
994 l++;
995 g[ l ] = fn;
996 h[ l ] = nd;
997 idx[ l ] = i1;
998 }
999 }
1000
1001 if( l == 0 )
1002 {
1003 break;
1004 }
1005
1006 difference_type ns = l;
1007 l = 1;
1008
1009 for( size_type i1 = _1s ; i1 <= _1e ; i1++ )
1010 {
1011 double nd = ( i1 - _1s + 1 ) * as;
1012 double len = h[ l ] - nd;
1013
1014 len = g[ l ] + len * len;
1015 for( ; l < ns ; l++ )
1016 {
1017 double len_ = h[ l + 1 ] - nd;
1018 len_ = g[ l + 1 ] + len_ * len_;
1019 if( len > len_ )
1020 {
1021 len = len_;
1022 }
1023 else
1024 {
1025 break;
1026 }
1027 }
1028 p[ i1 * diff ] = static_cast< value_type >( len );
1029 vp[ i1 * diff ] = vp[ idx[ l ] * diff ];
1030 }
1031 }
1032 }
1033
1034 delete [] idx;
1035 delete [] g;
1036 delete [] h;
1037 }
1038 };
1039
1040
1042 template < >
1043 struct __voronoi_distance_transform__< 1 >
1044 {
1045 template < class Array1, class Array2 >
1046 static void distance_transform( Array1 &voronoi, Array2 &dist, const _distance_utility_::__range__< 1 > &range, typename Array2::size_type thread_id = 0, typename Array2::size_type thread_num = 1 )
1047 {
1048 typedef typename Array2::size_type size_type;
1049 typedef typename Array1::value_type ivalue_type;
1050 typedef typename Array2::value_type value_type;
1051 typedef typename Array1::pointer ipointer;
1052 typedef typename Array2::pointer pointer;
1053 typedef typename Array2::difference_type difference_type;
1054
1055 const difference_type w = dist.width( );
1056 value_type max = type_limits< value_type >::maximum( );
1057 max = static_cast< value_type >( w ) < max ? static_cast< value_type >( w ) : max;
1058
1059 ivalue_type *val = new ivalue_type[ w ];
1060 value_type len;
1061
1062 size_type sx = range.begin1( );
1063 size_type sy = range.begin2( );
1064 size_type sz = range.begin3( );
1065 size_type ex = range.end1( );
1066 size_type ey = range.end2( );
1067 size_type ez = range.end3( );
1068
1069 for( size_type k = sz ; k <= ez ; k++ )
1070 {
1071 for( size_type j = sy + thread_id ; j <= ey ; j += thread_num )
1072 {
1073 ivalue_type *vp = val;
1074 ipointer ip = &voronoi( sx, j, k );
1075 pointer sp = &dist( sx, j, k );
1076 pointer ep = &dist( ex, j, k );
1077
1078 *vp++ = *ip;
1079 if( *ip++ == 0 )
1080 {
1081 len = max;
1082 sp[ 0 ] = len * len;
1083 }
1084 else
1085 {
1086 len = 0;
1087 }
1088
1089 for( pointer p = sp + 1 ; p <= ep ; p++, ip++ )
1090 {
1091 *vp++ = *ip;
1092 if( *ip == 0 )
1093 {
1094 len = len + 1 < max ? len + 1 : len;
1095 *p = len * len;
1096 *ip = *( ip - 1 );
1097 }
1098 else
1099 {
1100 len = 0;
1101 }
1102 }
1103
1104 vp--;
1105 ip -= 2;
1106
1107 if( *vp-- == 0 )
1108 {
1109 len = max;
1110 value_type len2 = len * len;
1111
1112 if( ep[ 0 ] > len2 )
1113 {
1114 ep[ 0 ] = len2;
1115 }
1116 }
1117 else
1118 {
1119 len = 0;
1120 }
1121
1122 for( pointer p = ep - 1 ; p >= sp ; p--, ip--, vp-- )
1123 {
1124 if( *vp == 0 )
1125 {
1126 len = len + 1 < max ? len + 1 : len;
1127 value_type len2 = len * len;
1128
1129 if( *p > len2 )
1130 {
1131 *p = len2;
1132 *ip = *( ip + 1 );
1133 }
1134 }
1135 else
1136 {
1137 len = 0;
1138 }
1139 }
1140 }
1141 }
1142
1143 delete [] val;
1144 }
1145 };
1146
1147
1148 template < class Array, int DIMENSION >
1149 bool compute_object_range( const Array &in, _distance_utility_::__range__< DIMENSION > &range )
1150 {
1151 typedef typename Array::size_type size_type;
1152 typedef typename Array::difference_type difference_type;
1153
1154 difference_type w = in.width( );
1155 difference_type h = in.height( );
1156 difference_type d = in.depth( );
1157
1158 difference_type sx, ex, sy, ey, sz, ez, i, j, k;
1159
1160 sx = sy = sz = 0;
1161 ex = w - 1;
1162 ey = h - 1;
1163 ez = d - 1;
1164
1165 for( k = 0 ; k < d ; k++ )
1166 {
1167 bool found = false;
1168 for( j = 0 ; j < h ; j++ )
1169 {
1170 for( i = 0 ; i < w ; i++ )
1171 {
1172 if( in( i, j, k ) != 0 )
1173 {
1174 found = true;
1175 break;
1176 }
1177 }
1178
1179 if( found )
1180 {
1181 break;
1182 }
1183 }
1184
1185 if( found )
1186 {
1187 break;
1188 }
1189 }
1190
1191 if( k >= d )
1192 {
1193 // 何も見つからなかった・・・
1194 return( false );
1195 }
1196
1197 sz = k;
1198
1199 for( k = d - 1 ; k > sz ; k-- )
1200 {
1201 bool found = false;
1202 for( j = 0 ; j < h ; j++ )
1203 {
1204 for( i = 0 ; i < w ; i++ )
1205 {
1206 if( in( i, j, k ) != 0 )
1207 {
1208 found = true;
1209 break;
1210 }
1211 }
1212
1213 if( found )
1214 {
1215 break;
1216 }
1217 }
1218
1219 if( found )
1220 {
1221 break;
1222 }
1223 }
1224
1225 ez = k;
1226
1227 for( j = 0 ; j < h ; j++ )
1228 {
1229 bool found = false;
1230 for( k = sz ; k <= ez ; k++ )
1231 {
1232 for( i = 0 ; i < w ; i++ )
1233 {
1234 if( in( i, j, k ) != 0 )
1235 {
1236 found = true;
1237 break;
1238 }
1239 }
1240
1241 if( found )
1242 {
1243 break;
1244 }
1245 }
1246
1247 if( found )
1248 {
1249 break;
1250 }
1251 }
1252
1253 sy = j;
1254
1255 for( j = h - 1 ; j > sy ; j-- )
1256 {
1257 bool found = false;
1258 for( k = sz ; k <= ez ; k++ )
1259 {
1260 for( i = 0 ; i < w ; i++ )
1261 {
1262 if( in( i, j, k ) != 0 )
1263 {
1264 found = true;
1265 break;
1266 }
1267 }
1268
1269 if( found )
1270 {
1271 break;
1272 }
1273 }
1274
1275 if( found )
1276 {
1277 break;
1278 }
1279 }
1280
1281 ey = j;
1282
1283 for( i = 0 ; i < w ; i++ )
1284 {
1285 bool found = false;
1286 for( k = sz ; k <= ez ; k++ )
1287 {
1288 for( j = sy ; j <= ey ; j++ )
1289 {
1290 if( in( i, j, k ) != 0 )
1291 {
1292 found = true;
1293 break;
1294 }
1295 }
1296
1297 if( found )
1298 {
1299 break;
1300 }
1301 }
1302
1303 if( found )
1304 {
1305 break;
1306 }
1307 }
1308
1309 sx = i;
1310
1311 for( i = w - 1 ; i > sx ; i-- )
1312 {
1313 bool found = false;
1314 for( k = sz ; k <= ez ; k++ )
1315 {
1316 for( j = sy ; j <= ey ; j++ )
1317 {
1318 if( in( i, j, k ) != 0 )
1319 {
1320 found = true;
1321 break;
1322 }
1323 }
1324
1325 if( found )
1326 {
1327 break;
1328 }
1329 }
1330
1331 if( found )
1332 {
1333 break;
1334 }
1335 }
1336
1337 ex = i;
1338
1339 // 一周り大きい範囲を設定する
1340 range.sx = sx > 0 ? sx - 1 : 0;
1341 range.ex = ex < w - 1 ? ex + 1 : w - 1;
1342 range.sy = sy > 0 ? sy - 1 : 0;
1343 range.ey = ey < h - 1 ? ey + 1 : h - 1;
1344 range.sz = sz > 0 ? sz - 1 : 0;
1345 range.ez = ez < d - 1 ? ez + 1 : d - 1;
1346
1347 return( true );
1348 }
1349
1350 template < class T >
1351 class calvin_distance_transform_thread : public mist::thread< calvin_distance_transform_thread< T > >
1352 {
1353 public:
1354 typedef mist::thread< calvin_distance_transform_thread< T > > base;
1355 typedef typename base::thread_exit_type thread_exit_type;
1356 typedef typename T::size_type size_type;
1357 typedef typename T::value_type value_type;
1358
1359 protected:
1360 size_t thread_id_;
1361 size_t thread_num_;
1362
1363 // 入出力用の画像へのポインタ
1364 T *in_;
1365 size_type axis_;
1366 _distance_utility_::__range__< 1 > range_;
1367
1368 public:
1369 void setup_parameters( T &in, const _distance_utility_::__range__< 1 > &range, size_type axis, size_type thread_id, size_type thread_num )
1370 {
1371 in_ = &in;
1372 range_ = range;
1373 axis_ = axis;
1374 thread_id_ = thread_id;
1375 thread_num_ = thread_num;
1376 }
1377
1378 void setup_axis( size_type axis )
1379 {
1380 axis_ = axis;
1381 }
1382
1383 const calvin_distance_transform_thread& operator =( const calvin_distance_transform_thread &p )
1384 {
1385 if( &p != this )
1386 {
1387 base::operator =( p );
1388 thread_id_ = p.thread_id_;
1389 thread_num_ = p.thread_num_;
1390 in_ = p.in_;
1391 range_ = p.range_;
1392 axis_ = p.axis_;
1393 }
1394 return( *this );
1395 }
1396
1397 calvin_distance_transform_thread( size_type id = 0, size_type num = 1 )
1398 : thread_id_( id ), thread_num_( num ), in_( NULL ), axis_( 0 )
1399 {
1400 }
1401 calvin_distance_transform_thread( const calvin_distance_transform_thread &p )
1402 : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ), in_( NULL ), axis_( p.axis_ )
1403 {
1404 }
1405
1406 protected:
1407 // 継承した先で必ず実装されるスレッド関数
1408 virtual thread_exit_type thread_function( )
1409 {
1410 switch( axis_ )
1411 {
1412 case 0:
1413 __distance_transform__< 1 >::distance_transform( *in_, _distance_utility_::__range__< 1 >( range_ ), thread_id_, thread_num_ );
1414 break;
1415
1416 case 1:
1417 __distance_transform__< 2 >::distance_transform( *in_, _distance_utility_::__range__< 2 >( range_ ), thread_id_, thread_num_ );
1418 break;
1419
1420 case 2:
1421 default:
1422 __distance_transform__< 3 >::distance_transform( *in_, _distance_utility_::__range__< 3 >( range_ ), thread_id_, thread_num_ );
1423 break;
1424 }
1425 return( true );
1426 }
1427 };
1428
1429
1430 template < class T1, class T2 >
1431 class calvin_voronoi_distance_transform_thread : public mist::thread< calvin_voronoi_distance_transform_thread< T1, T2 > >
1432 {
1433 public:
1434 typedef mist::thread< calvin_voronoi_distance_transform_thread< T1, T2 > > base;
1435 typedef typename base::thread_exit_type thread_exit_type;
1436 typedef typename T1::size_type size_type;
1437 typedef typename T1::value_type value_type;
1438
1439 private:
1440 size_t thread_id_;
1441 size_t thread_num_;
1442
1443 // 入出力用の画像へのポインタ
1444 T1 *voronoi_;
1445 T2 *dist_;
1446 size_type axis_;
1447 _distance_utility_::__range__< 1 > range_;
1448
1449 public:
1450 void setup_parameters( T1 &voronoi, T2 &dist, const _distance_utility_::__range__< 1 > &range, size_type axis, size_type thread_id, size_type thread_num )
1451 {
1452 voronoi_ = &voronoi;
1453 dist_ = &dist;
1454 range_ = range;
1455 axis_ = axis;
1456 thread_id_ = thread_id;
1457 thread_num_ = thread_num;
1458 }
1459
1460 void setup_axis( size_type axis )
1461 {
1462 axis_ = axis;
1463 }
1464
1465 const calvin_voronoi_distance_transform_thread& operator =( const calvin_voronoi_distance_transform_thread &p )
1466 {
1467 if( &p != this )
1468 {
1469 base::operator =( p );
1470 thread_id_ = p.thread_id_;
1471 thread_num_ = p.thread_num_;
1472 voronoi_ = p.voronoi_;
1473 dist_ = p.dist_;
1474 range_ = p.range_;
1475 axis_ = p.axis_;
1476 }
1477 return( *this );
1478 }
1479
1480 calvin_voronoi_distance_transform_thread( size_type id = 0, size_type num = 1 )
1481 : thread_id_( id ), thread_num_( num ), voronoi_( NULL ), dist_( NULL ), axis_( 0 )
1482 {
1483 }
1484 calvin_voronoi_distance_transform_thread( const calvin_voronoi_distance_transform_thread &p )
1485 : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ), voronoi_( NULL ), dist_( NULL ), axis_( p.axis_ )
1486 {
1487 }
1488
1489 protected:
1490 // 継承した先で必ず実装されるスレッド関数
1491 virtual thread_exit_type thread_function( )
1492 {
1493 switch( axis_ )
1494 {
1495 case 0:
1496 __voronoi_distance_transform__< 1 >::distance_transform( *voronoi_, *dist_, _distance_utility_::__range__< 1 >( range_ ), thread_id_, thread_num_ );
1497 break;
1498
1499 case 1:
1500 __voronoi_distance_transform__< 2 >::distance_transform( *voronoi_, *dist_, _distance_utility_::__range__< 2 >( range_ ), thread_id_, thread_num_ );
1501 break;
1502
1503 case 2:
1504 default:
1505 __voronoi_distance_transform__< 3 >::distance_transform( *voronoi_, *dist_, _distance_utility_::__range__< 3 >( range_ ), thread_id_, thread_num_ );
1506 break;
1507 }
1508 return( true );
1509 }
1510 };
1511 }
1512
1513
1514 namespace _meijster_distance_transform_
1515 {
1516 struct MDT
1517 {
1518 template < class T, class Val >
1519 static T Func( const T &x, const T &i, const Val &gi )
1520 {
1521 return( static_cast< T >( std::abs( static_cast< int >( x - i ) ) + gi ) );
1522 }
1523
1524 template < class T, class Val >
1525 static T Sep( const T &i, const T &u, const Val &gi, const Val &gu )
1526 {
1527 if( gu >= gi + u - i )
1528 {
1529 return( type_limits< T >::maximum( ) );
1530 }
1531 else if( gi > gu + u - i )
1532 {
1533 return( type_limits< T >::minimum( ) );
1534 }
1535 else
1536 {
1537 return( ( gu - gi + u + i ) / 2 );
1538 }
1539 }
1540 };
1541
1542 struct CDT
1543 {
1544 template < class T, class Val >
1545 static T Func( const T &x, const T &i, const Val &gi )
1546 {
1547 T v = std::abs( static_cast< int >( x - i ) );
1548 return( static_cast< T >( v > gi ? v : gi ) );
1549 }
1550
1551 template < class T, class Val >
1552 static T Sep( const T &i, const T &u, const Val &gi, const Val &gu )
1553 {
1554 if( gi <= gu )
1555 {
1556 T igu = static_cast< T >( i + gu );
1557 T iu = static_cast< T >( ( i + u ) / 2 );
1558 return( igu > iu ? igu : iu );
1559 }
1560 else
1561 {
1562 T ugi = static_cast< T >( u - gi );
1563 T iu = static_cast< T >( ( i + u ) / 2 );
1564 return( ugi < iu ? ugi : iu );
1565 }
1566 }
1567 };
1568
1569 template < class Array >
1570 void distance_transform_x( Array &in, typename Array::size_type thread_id = 0, typename Array::size_type thread_num = 1 )
1571 {
1572 typedef typename Array::difference_type difference_type;
1573 typedef typename Array::value_type value_type;
1574
1575 difference_type i, j, k;
1576
1577 const difference_type w = in.width( );
1578 const difference_type h = in.height( );
1579 const difference_type d = in.depth( );
1580
1581 value_type max = type_limits< value_type >::maximum( );
1582
1583 for( k = 0 ; k < d ; k++ )
1584 {
1585 for( i = thread_id ; i < w ; i += thread_num )
1586 {
1587 if( in( i, 0, k ) != 0 )
1588 {
1589 in( i, 0, k ) = max;
1590 }
1591
1592 for( j = 1 ; j < h ; j++ )
1593 {
1594 if( in( i, j, k ) != 0 )
1595 {
1596 in( i, j, k ) = in( i, j - 1, k ) + 1;
1597 }
1598 }
1599
1600 for( j = h - 2 ; j >= 0 ; j-- )
1601 {
1602 if( in( i, j + 1, k ) < in( i, j, k ) )
1603 {
1604 in( i, j, k ) = in( i, j + 1, k ) + 1;
1605 }
1606 }
1607 }
1608 }
1609 }
1610
1611 template < class Array, class Metric >
1612 void distance_transform_y( Array &in, Metric __dmy__, typename Array::size_type thread_id = 0, typename Array::size_type thread_num = 1 )
1613 {
1614 typedef typename Array::pointer pointer;
1615 typedef typename Array::difference_type difference_type;
1616 typedef typename Array::value_type value_type;
1617
1618 const difference_type w = in.width( );
1619 const difference_type h = in.height( );
1620 const difference_type d = in.depth( );
1621
1622 difference_type i, j, k, q;
1623 difference_type *s = new difference_type[ w ];
1624 difference_type *t = new difference_type[ w ];
1625
1626 for( k = 0 ; k < d ; k++ )
1627 {
1628 for( j = thread_id ; j < h ; j += thread_num )
1629 {
1630 q = s[ 0 ] = t[ 0 ] = 0;
1631
1632 pointer pin = &in( 0, j, k );
1633
1634 for( i = 1 ; i < w ; i++ )
1635 {
1636 value_type val = in( i, j, k );
1637
1638 while( q >= 0 && Metric::Func( t[ q ], s[ q ], pin[ s[ q ] ] ) > Metric::Func( t[ q ], i, val ) )
1639 {
1640 q--;
1641 }
1642
1643 if( q < 0 )
1644 {
1645 q = 0;
1646 s[ 0 ] = i;
1647 }
1648 else
1649 {
1650 difference_type tmp = Metric::Sep( s[ q ], i, pin[ s[ q ] ], val );
1651
1652 if( tmp < w - 1 )
1653 {
1654 q++;
1655 s[ q ] = i;
1656 t[ q ] = tmp + 1;
1657 }
1658 }
1659 }
1660
1661 for( i = w - 1 ; i >= 0 ; i-- )
1662 {
1663 pin[ i ] = static_cast< value_type >( Metric::Func( i, s[ q ], pin[ s[ q ] ] ) );
1664 if( i == t[ q ] )
1665 {
1666 q--;
1667 }
1668 }
1669 }
1670 }
1671
1672 delete [] s;
1673 delete [] t;
1674 }
1675
1676 template < class Array, class Metric >
1677 void distance_transform_z( Array &in, Metric __dmy__, typename Array::size_type thread_id = 0, typename Array::size_type thread_num = 1 )
1678 {
1679 typedef typename Array::pointer pointer;
1680 typedef typename Array::difference_type difference_type;
1681 typedef typename Array::value_type value_type;
1682
1683 const difference_type w = in.width( );
1684 const difference_type h = in.height( );
1685 const difference_type d = in.depth( );
1686
1687 difference_type i, j, k, q;
1688 difference_type *s = new difference_type[ d ];
1689 difference_type *t = new difference_type[ d ];
1690 difference_type diff = &in( 0, 0, 1 ) - &in( 0, 0, 0 );
1691
1692 for( j = 0 ; j < h ; j++ )
1693 {
1694 for( i = thread_id ; i < w ; i += thread_num )
1695 {
1696 q = s[ 0 ] = t[ 0 ] = 0;
1697
1698 pointer pin = &in( i, j, 0 );
1699
1700 for( k = 1 ; k < d ; k++ )
1701 {
1702 value_type val = in( i, j, k );
1703
1704 while( q >= 0 && Metric::Func( t[ q ], s[ q ], pin[ s[ q ] * diff ] ) > Metric::Func( t[ q ], k, val ) )
1705 {
1706 q--;
1707 }
1708
1709 if( q < 0 )
1710 {
1711 q = 0;
1712 s[ 0 ] = k;
1713 }
1714 else
1715 {
1716 difference_type tmp = Metric::Sep( s[ q ], k, pin[ s[ q ] * diff ], val );
1717
1718 if( tmp < d - 1 )
1719 {
1720 q++;
1721 s[ q ] = k;
1722 t[ q ] = tmp + 1;
1723 }
1724 }
1725 }
1726
1727 for( k = d - 1 ; k >= 0 ; k-- )
1728 {
1729 pin[ k * diff ] = static_cast< value_type >( Metric::Func( k, s[ q ], pin[ s[ q ] * diff ] ) );
1730 if( k == t[ q ] )
1731 {
1732 q--;
1733 }
1734 }
1735 }
1736 }
1737
1738 delete [] s;
1739 delete [] t;
1740 }
1741
1742
1743 template < class T, class Metric >
1744 class distance_transform_thread : public mist::thread< distance_transform_thread< T, Metric > >
1745 {
1746 public:
1747 typedef mist::thread< distance_transform_thread< T, Metric > > base;
1748 typedef typename base::thread_exit_type thread_exit_type;
1749 typedef typename T::size_type size_type;
1750 typedef typename T::value_type value_type;
1751
1752 private:
1753 size_t thread_id_;
1754 size_t thread_num_;
1755
1756 // 入出力用の画像へのポインタ
1757 T *in_;
1758 size_type axis_;
1759 Metric dmy_;
1760
1761 public:
1762 void setup_parameters( T &in, size_type axis, size_type thread_id, size_type thread_num )
1763 {
1764 in_ = &in;
1765 axis_ = axis;
1766 thread_id_ = thread_id;
1767 thread_num_ = thread_num;
1768 }
1769
1770 void setup_axis( size_type axis )
1771 {
1772 axis_ = axis;
1773 }
1774
1775 const distance_transform_thread& operator =( const distance_transform_thread &p )
1776 {
1777 if( &p != this )
1778 {
1779 base::operator =( p );
1780 thread_id_ = p.thread_id_;
1781 thread_num_ = p.thread_num_;
1782 in_ = p.in_;
1783 axis_ = p.axis_;
1784 }
1785 return( *this );
1786 }
1787
1788 distance_transform_thread( size_type id = 0, size_type num = 1 )
1789 : thread_id_( id ), thread_num_( num ), in_( NULL ), axis_( 0 )
1790 {
1791 }
1792 distance_transform_thread( const distance_transform_thread &p )
1793 : base( p ), thread_id_( p.thread_id_ ), thread_num_( p.thread_num_ ), in_( NULL ), axis_( p.axis_ )
1794 {
1795 }
1796
1797 protected:
1798 // 継承した先で必ず実装されるスレッド関数
1799 virtual thread_exit_type thread_function( )
1800 {
1801 switch( axis_ )
1802 {
1803 case 1:
1804 distance_transform_y( *in_, dmy_, thread_id_, thread_num_ );
1805 break;
1806
1807 case 2:
1808 distance_transform_z( *in_, dmy_, thread_id_, thread_num_ );
1809 break;
1810
1811 case 0:
1812 default:
1813 distance_transform_x( *in_, thread_id_, thread_num_ );
1814 break;
1815 }
1816 return( true );
1817 }
1818 };
1819 }
1820
1821
1822
1830
1831
1839
1840
1841
1843 namespace saito
1844 {
1867 template < class Array1, class Array2 >
1868 void distance_transform( const Array1 &in, Array2 &out, double max_length = -1.0, typename Array1::size_type thread_num = 0 )
1869 {
1870 typedef typename Array2::size_type size_type;
1871 typedef typename Array2::value_type value_type;
1872 typedef __saito__::saito_distance_transform_thread< Array2 > saito_distance_transform_thread;
1873
1874 if( thread_num == 0 )
1875 {
1876 thread_num = static_cast< size_type >( get_cpu_num( ) );
1877 }
1878
1879 out.resize( in.size1( ), in.size2( ), in.size3( ) );
1880 out.reso1( in.reso1( ) );
1881 out.reso2( in.reso2( ) );
1882 out.reso3( in.reso3( ) );
1883
1884 size_type i;
1885
1886 for( i = 0 ; i < in.size( ) ; i++ )
1887 {
1888 out[ i ] = static_cast< value_type >( in[ i ] != 0 ? 1 : 0 );
1889 }
1890
1891 saito_distance_transform_thread *thread = new saito_distance_transform_thread[ thread_num ];
1892
1893 if( in.width( ) > 1 )
1894 {
1895 // X軸方向の処理
1896 for( i = 0 ; i < thread_num ; i++ )
1897 {
1898 thread[ i ].setup_parameters( out, max_length, 0, i, thread_num );
1899 }
1900
1901 do_threads_( thread, thread_num );
1902 }
1903
1904 if( in.height( ) > 1 )
1905 {
1906 // Y軸方向の処理
1907 for( i = 0 ; i < thread_num ; i++ )
1908 {
1909 thread[ i ].setup_parameters( out, max_length, 1, i, thread_num );
1910 }
1911
1912 do_threads_( thread, thread_num );
1913 }
1914
1915 if( in.depth( ) > 1 )
1916 {
1917 // Y軸方向の処理
1918 for( i = 0 ; i < thread_num ; i++ )
1919 {
1920 thread[ i ].setup_parameters( out, max_length, 2, i, thread_num );
1921 }
1922
1923 do_threads_( thread, thread_num );
1924 }
1925
1926 delete [] thread;
1927 }
1928
1929
1948 template < class Array1, class Array2 >
1949 void inverse_distance_transform( const Array1 &in, Array2 &out, typename Array1::size_type thread_num = 0 )
1950 {
1951 typedef typename Array2::size_type size_type;
1952 typedef typename Array2::difference_type difference_type;
1953 typedef typename Array2::value_type value_type;
1954 typedef __saito__::saito_inverse_distance_transform_thread< Array2 > saito_inverse_distance_transform_thread;
1955
1956 if( thread_num == 0 )
1957 {
1958 thread_num = static_cast< size_type >( get_cpu_num( ) );
1959 }
1960
1961 out = in;
1962 //out.resize( in.size1( ), in.size2( ), in.size3( ) );
1963 //out.reso1( in.reso1( ) );
1964 //out.reso2( in.reso2( ) );
1965 //out.reso3( in.reso3( ) );
1966
1967 //for( size_type i = 0 ; i < in.size( ) ; i++ )
1968 //{
1969 // out[ i ] = static_cast< value_type >( in[ i ] );
1970 //}
1971
1972 _distance_utility_::__range__< 0 > object_range;
1973 object_range.sx = 0;
1974 object_range.ex = in.size1( ) - 1;
1975 object_range.sy = 0;
1976 object_range.ey = in.size2( ) - 1;
1977 object_range.sz = 0;
1978 object_range.ez = in.size3( ) - 1;
1979
1980 saito_inverse_distance_transform_thread *thread = new saito_inverse_distance_transform_thread[ thread_num ];
1981
1982 if( in.width( ) > 1 )
1983 {
1984 // X軸方向の処理
1985 for( size_type i = 0 ; i < thread_num ; i++ )
1986 {
1987 thread[ i ].setup_parameters( out, object_range, 0, i, thread_num );
1988 }
1989
1990 do_threads_( thread, thread_num );
1991 }
1992
1993 if( in.height( ) > 1 )
1994 {
1995 // Y軸方向の処理
1996 for( size_type i = 0 ; i < thread_num ; i++ )
1997 {
1998 thread[ i ].setup_parameters( out, object_range, 1, i, thread_num );
1999 }
2000
2001 do_threads_( thread, thread_num );
2002 }
2003
2004 if( in.depth( ) > 1 )
2005 {
2006 // Z軸方向の処理
2007 for( size_type i = 0 ; i < thread_num ; i++ )
2008 {
2009 thread[ i ].setup_parameters( out, object_range, 2, i, thread_num );
2010 }
2011
2012 do_threads_( thread, thread_num );
2013 }
2014
2015 delete [] thread;
2016
2017 // 逆距離変換の結果得られる距離値は正しいユークリッド2乗距離ではないため，結果をもとに2値図形に直す．
2018 for( size_type i = 0 ; i < out.size( ) ; i++ )
2019 {
2020 out[ i ] = out[ i ] != 0 ? 1 : 0;
2021 }
2022 }
2023 }
2024
2025
2027 namespace calvin
2028 {
2050 template < class Array1, class Array2 >
2051 void distance_transform( const Array1 &in, Array2 &out, typename Array1::size_type thread_num = 0 )
2052 {
2053 typedef typename Array2::size_type size_type;
2054 typedef typename Array2::value_type value_type;
2055 typedef __calvin__::calvin_distance_transform_thread< Array2 > calvin_distance_transform_thread;
2056
2057 if( thread_num == 0 )
2058 {
2059 thread_num = static_cast< size_type >( get_cpu_num( ) );
2060 }
2061
2062 out.resize( in.size1( ), in.size2( ), in.size3( ) );
2063 out.reso1( in.reso1( ) );
2064 out.reso2( in.reso2( ) );
2065 out.reso3( in.reso3( ) );
2066
2067 size_type i;
2068
2069 for( i = 0 ; i < in.size( ) ; i++ )
2070 {
2071 out[ i ] = static_cast< value_type >( in[ i ] != 0 ? 1 : 0 );
2072 }
2073
2074 _distance_utility_::__range__< 0 > object_range;
2075 if( !__calvin__::compute_object_range( out, object_range ) )
2076 {
2077 // １つも１画素が見つからなかったので終了する
2078 return;
2079 }
2080
2081 calvin_distance_transform_thread *thread = new calvin_distance_transform_thread[ thread_num ];
2082
2083 if( in.width( ) > 1 )
2084 {
2085 // X軸方向の処理
2086 for( i = 0 ; i < thread_num ; i++ )
2087 {
2088 thread[ i ].setup_parameters( out, object_range, 0, i, thread_num );
2089 }
2090
2091 do_threads_( thread, thread_num );
2092 }
2093
2094 if( in.height( ) > 1 )
2095 {
2096 // Y軸方向の処理
2097 for( i = 0 ; i < thread_num ; i++ )
2098 {
2099 thread[ i ].setup_parameters( out, object_range, 1, i, thread_num );
2100 }
2101
2102 do_threads_( thread, thread_num );
2103 }
2104
2105 if( in.depth( ) > 1 )
2106 {
2107 // Z軸方向の処理
2108 for( i = 0 ; i < thread_num ; i++ )
2109 {
2110 thread[ i ].setup_parameters( out, object_range, 2, i, thread_num );
2111 }
2112
2113 do_threads_( thread, thread_num );
2114 }
2115
2116 delete [] thread;
2117 }
2118
2119
2139 template < class Array1, class Array2 >
2140 void voronoi_distance_transform( Array1 &voronoi, Array2 &dist, typename Array1::size_type thread_num = 0 )
2141 {
2142 typedef typename Array2::size_type size_type;
2143 typedef typename Array2::value_type value_type;
2144 typedef __calvin__::calvin_voronoi_distance_transform_thread< Array1, Array2 > calvin_voronoi_distance_transform_thread;
2145
2146 if( thread_num == 0 )
2147 {
2148 thread_num = static_cast< size_type >( get_cpu_num( ) );
2149 }
2150
2151 dist.resize( voronoi.size1( ), voronoi.size2( ), voronoi.size3( ) );
2152 dist.reso1( voronoi.reso1( ) );
2153 dist.reso2( voronoi.reso2( ) );
2154 dist.reso3( voronoi.reso3( ) );
2155
2156 size_type i;
2157
2158 _distance_utility_::__range__< 0 > object_range;
2159 object_range.sx = 0;
2160 object_range.ex = voronoi.size1( ) - 1;
2161 object_range.sy = 0;
2162 object_range.ey = voronoi.size2( ) - 1;
2163 object_range.sz = 0;
2164 object_range.ez = voronoi.size3( ) - 1;
2165
2166 calvin_voronoi_distance_transform_thread *thread = new calvin_voronoi_distance_transform_thread[ thread_num ];
2167
2168 if( voronoi.width( ) > 1 )
2169 {
2170 // X軸方向の処理
2171 for( i = 0 ; i < thread_num ; i++ )
2172 {
2173 thread[ i ].setup_parameters( voronoi, dist, object_range, 0, i, thread_num );
2174 }
2175
2176 do_threads_( thread, thread_num );
2177 }
2178
2179 if( voronoi.height( ) > 1 )
2180 {
2181 // Y軸方向の処理
2182 for( i = 0 ; i < thread_num ; i++ )
2183 {
2184 thread[ i ].setup_parameters( voronoi, dist, object_range, 1, i, thread_num );
2185 }
2186
2187 do_threads_( thread, thread_num );
2188 }
2189
2190 if( voronoi.depth( ) > 1 )
2191 {
2192 // Z軸方向の処理
2193 for( i = 0 ; i < thread_num ; i++ )
2194 {
2195 thread[ i ].setup_parameters( voronoi, dist, object_range, 2, i, thread_num );
2196 }
2197
2198 do_threads_( thread, thread_num );
2199 }
2200
2201 delete [] thread;
2202 }
2203
2221 template < class Array >
2222 void voronoi_transform( Array &voronoi, typename Array::size_type thread_num = 0 )
2223 {
2224 double ax = voronoi.reso1( );
2225 double ay = voronoi.reso2( );
2226 double az = voronoi.reso3( );
2227 if( ax == ay && ay == az )
2228 {
2229 typename Array::template rebind< unsigned int >::other dist;
2230 voronoi_distance_transform( voronoi, dist, thread_num );
2231 }
2232 else
2233 {
2234 typename Array::template rebind< float >::other dist;
2235 voronoi_distance_transform( voronoi, dist, thread_num );
2236 }
2237 }
2238 }
2239
2240
2242 namespace euclidean
2243 {
2265 template < class Array1, class Array2 >
2266 void distance_transform( const Array1 &in, Array2 &out, typename Array1::size_type thread_num = 0 )
2267 {
2268 calvin::distance_transform( in, out, thread_num );
2269 }
2270
2271
2290 template < class Array1, class Array2 >
2291 void inverse_distance_transform( const Array1 &in, Array2 &out, typename Array1::size_type thread_num = 0 )
2292 {
2293 saito::inverse_distance_transform( in, out, thread_num );
2294 }
2295
2296
2304
2305
2325 template < class Array1, class Array2 >
2326 void voronoi_distance_transform( Array1 &voronoi, Array2 &dist, typename Array1::size_type thread_num = 0 )
2327 {
2328 calvin::voronoi_distance_transform( voronoi, dist, thread_num );
2329 }
2330
2348 template < class Array >
2349 void voronoi_transform( Array &voronoi, typename Array::size_type thread_num = 0 )
2350 {
2351 calvin::voronoi_transform( voronoi, thread_num );
2352 }
2353
2355 // ボロノイ分割グループの終わり
2356 }
2357
2359 // ユークリッド距離変換グループの終わり
2360
2361
2369
2371 namespace meijster
2372 {
2388 template < class Array1, class Array2, class Metric >
2389 void distance_transform( const Array1 &in, Array2 &out, Metric __metric__, typename Array1::size_type thread_num = 0 )
2390 {
2391 typedef typename Array2::size_type size_type;
2392 typedef typename Array2::value_type value_type;
2393 typedef _meijster_distance_transform_::distance_transform_thread< Array2, Metric > meijster_distance_transform_thread;
2394
2395 if( thread_num == 0 )
2396 {
2397 thread_num = static_cast< size_type >( get_cpu_num( ) );
2398 }
2399
2400 out.resize( in.size1( ), in.size2( ), in.size3( ) );
2401 out.reso1( in.reso1( ) );
2402 out.reso2( in.reso2( ) );
2403 out.reso3( in.reso3( ) );
2404
2405 size_type i;
2406
2407 for( i = 0 ; i < in.size( ) ; i++ )
2408 {
2409 out[ i ] = static_cast< value_type >( in[ i ] != 0 ? 1 : 0 );
2410 }
2411
2412 meijster_distance_transform_thread *thread = new meijster_distance_transform_thread[ thread_num ];
2413
2414 if( in.width( ) > 1 )
2415 {
2416 // X軸方向の処理
2417 for( i = 0 ; i < thread_num ; i++ )
2418 {
2419 thread[ i ].setup_parameters( out, 0, i, thread_num );
2420 }
2421
2422 do_threads_( thread, thread_num );
2423 }
2424
2425 if( in.height( ) > 1 )
2426 {
2427 // Y軸方向の処理
2428 for( i = 0 ; i < thread_num ; i++ )
2429 {
2430 thread[ i ].setup_parameters( out, 1, i, thread_num );
2431 }
2432
2433 do_threads_( thread, thread_num );
2434 }
2435
2436 if( in.depth( ) > 1 )
2437 {
2438 // Z軸方向の処理
2439 for( i = 0 ; i < thread_num ; i++ )
2440 {
2441 thread[ i ].setup_parameters( out, 2, i, thread_num );
2442 }
2443
2444 do_threads_( thread, thread_num );
2445 }
2446
2447 delete [] thread;
2448 }
2449 }
2450
2452 // A. Meijster による距離変換グループの終わり
2453
2454
2455
2463
2464
2466 namespace manhattan
2467 {
2482 template < class Array1, class Array2 >
2483 void distance_transform( const Array1 &in, Array2 &out, typename Array1::size_type thread_num = 0 )
2484 {
2485 _meijster_distance_transform_::MDT metric;
2486 meijster::distance_transform( in, out, metric, thread_num );
2487 }
2488 }
2489
2490
2492 namespace cityblock
2493 {
2508 template < class Array1, class Array2 >
2509 void distance_transform( const Array1 &in, Array2 &out, typename Array1::size_type thread_num = 0 )
2510 {
2511 _meijster_distance_transform_::MDT metric;
2512 meijster::distance_transform( in, out, metric, thread_num );
2513 }
2514 }
2515
2517 // マンハッタン距離（シティーブロック距離）変換
2518
2519
2520
2528
2530 namespace chessboard
2531 {
2546 template < class Array1, class Array2 >
2547 void distance_transform( const Array1 &in, Array2 &out, typename Array1::size_type thread_num = 0 )
2548 {
2549 _meijster_distance_transform_::CDT metric;
2550 meijster::distance_transform( in, out, metric, thread_num );
2551 }
2552 }
2553
2555 // チェスボード距離距離変換
2556
2557
2558
2566
2574 template < class Array1, class Array2 >
2575 void skeleton( const Array1 &in, Array2 &out )
2576 {
2577 typedef typename Array1::size_type size_type;
2578 typedef typename Array1::difference_type difference_type;
2579 typedef typename Array2::value_type value_type;
2580
2581 out.resize( in.size1( ), in.size2( ), in.size3( ) );
2582 out.reso1( in.reso1( ) );
2583 out.reso2( in.reso2( ) );
2584 out.reso3( in.reso3( ) );
2585
2586 // 作業用配列
2587 array3< difference_type > temp( in.width( ), in.height( ), in.depth( ) );
2588
2589 for( size_type i = 0 ; i < in.size( ) ; i++ )
2590 {
2591 temp[ i ] = -1;
2592 }
2593
2594 difference_type w = in.width( );
2595 difference_type h = in.height( );
2596 difference_type d = in.depth( );
2597 for( difference_type k = 0 ; k < d ; k++ )
2598 {
2599 for( difference_type j = 0 ; j < h ; j++ )
2600 {
2601 for( difference_type i = 0 ; i < w ; i++ )
2602 {
2603 difference_type index = &in( i, j, k ) - &in[ 0 ];
2604 difference_type rr = in[ index ];
2605
2606 if( rr == 0 )
2607 {
2608 continue;
2609 }
2610
2611 difference_type rx = static_cast< difference_type >( std::ceil( std::sqrt( static_cast< double >( rr ) ) ) );
2612 difference_type ry = rx;
2613 difference_type rz = rx;
2614
2615 for( difference_type z = -rz ; z <= rz ; z++ )
2616 {
2617 difference_type pz = k + z;
2618 if( pz < 0 || d <= pz )
2619 {
2620 continue;
2621 }
2622
2623 difference_type zz = z * z;
2624
2625 for( difference_type y = -ry ; y <= ry ; y++ )
2626 {
2627 difference_type py = j + y;
2628 if( py < 0 || h <= py )
2629 {
2630 continue;
2631 }
2632
2633 difference_type yy = y * y;
2634
2635 for( difference_type x = -rx ; x <= rx ; x++ )
2636 {
2637 difference_type px = i + x;
2638 if( px < 0 || w <= px )
2639 {
2640 continue;
2641 }
2642
2643 difference_type xx = x * x;
2644
2645 if( xx + yy + zz < rr )
2646 {
2647 difference_type &p = temp( px, py, pz );
2648 if( p < 0 || in[ p ] < rr )
2649 {
2650 p = index;
2651 }
2652 }
2653 }
2654 }
2655 }
2656 }
2657 }
2658 }
2659
2660 for( size_type i = 0 ; i < in.size( ) ; i++ )
2661 {
2662 if( temp[ i ] > 0 )
2663 {
2664 temp[ temp[ i ] ] = -2;
2665 }
2666 }
2667
2668 for( size_type i = 0 ; i < in.size( ) ; i++ )
2669 {
2670 if( temp[ i ] != -2 )
2671 {
2672 out[ i ] = 0;
2673 }
2674 else
2675 {
2676 out[ i ] = static_cast< value_type >( in[ i ] );
2677 }
2678 }
2679 }
2680
2682 // スケルトン抽出グループの終わり
2683
2684
2686 // 距離変換グループの終わり
2687
2688
2689
2690
2691 // mist名前空間の終わり
2692 _MIST_END
2693
2694
2695 #endif // __INCLUDE_MIST_DISTANCE_TRANSFORM__
2696