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

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

region_growing.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_REGIONGROWING__
35 #define __INCLUDE_MIST_REGIONGROWING__
36
37
38 #ifndef __INCLUDE_MIST_H__
39 #include "../mist.h"
40 #endif
41
42 #ifndef __INCLUDE_MIST_LIMITS__
43 #include "../limits.h"
44 #endif
45
46 #ifndef __INCLUDE_MIST_VECTOR__
47 #include "../vector.h"
48 #endif
49
50 #include <deque>
51 #include <vector>
52 #include <cmath>
53
54
55 // mist名前空間の始まり
56 _MIST_BEGIN
57
58
59
82
83
88 namespace region_growing_utility
89 {
91 enum expand_mode_type
92 {
93 NC4,
94 NC8,
95 NC6,
96 NC18,
97 NC26,
98 ALL
99 };
100
101
103 struct pixel
104 {
105 typedef size_t size_type;
106 typedef ptrdiff_t difference_type;
107
108 size_type width( ) const { return( 1 ); }
109 size_type height( ) const { return( 1 ); }
110 size_type depth( ) const { return( 1 ); }
111
112 expand_mode_type __expand_mode__;
113
115 bool operator ()( difference_type /* i */, difference_type /* j */, difference_type /* k */ ) const { return( true ); }
116
118 expand_mode_type expand_mode( ) const { return( __expand_mode__ ); }
119
121 pixel( ) : __expand_mode__( NC4 ){ }
122
124 pixel( expand_mode_type mode ) : __expand_mode__( mode ){ }
125 };
126
127
129 struct voxel
130 {
131 typedef size_t size_type;
132 typedef ptrdiff_t difference_type;
133
134 size_type width( ) const { return( 1 ); }
135 size_type height( ) const { return( 1 ); }
136 size_type depth( ) const { return( 1 ); }
137
138 expand_mode_type __expand_mode__;
139
141 bool operator ()( difference_type /* i */, difference_type /* j */, difference_type /* k */ ) const { return( true ); }
142
144 expand_mode_type expand_mode( ) const { return( __expand_mode__ ); }
145
147 voxel( ) : __expand_mode__( NC6 ){ }
148
150 voxel( expand_mode_type mode ) : __expand_mode__( mode ){ }
151 };
152
153
155 class circle
156 {
157 public:
158 typedef size_t size_type;
159 typedef ptrdiff_t difference_type;
160
161 private:
162 double radius_;
163 double resoX_;
164 double resoY_;
165 size_type width_;
166 size_type height_;
167
168 public:
169 size_type width( ) const { return( width_ ); }
170 size_type height( ) const { return( height_ ); }
171 size_type depth( ) const { return( 1 ); }
172
174 bool operator ()( difference_type i, difference_type j, difference_type k ) const
175 {
176 if( k != 0 )
177 {
178 return( false );
179 }
180 else
181 {
182 return( i * i * resoX_ * resoX_ + j * j * resoY_ * resoY_ <= radius_ * radius_ );
183 }
184 }
185
187 expand_mode_type expand_mode( ) const { return( NC4 ); }
188
190 circle( double radius, double resoX = 1.0, double resoY = 1.0, bool radiusInPhysicalCoords = false ) : radius_( radius )
191 {
192 if( radiusInPhysicalCoords )
193 {
194 resoX_ = resoX;
195 resoY_ = resoY;
196 }
197 else
198 {
199 double max_reso = resoX > resoY ? resoX: resoY;
200
201 resoX_ = resoX / max_reso;
202 resoY_ = resoY / max_reso;
203 }
204 width_ = static_cast< size_type >( std::ceil( radius_ / resoX_ ) ) * 2 + 1;
205 height_ = static_cast< size_type >( std::ceil( radius_ / resoY_ ) ) * 2 + 1;
206 }
207 };
208
209
211 class sphere
212 {
213 public:
214 typedef size_t size_type;
215 typedef ptrdiff_t difference_type;
216
217 private:
218 double radius_;
219 double resoX_;
220 double resoY_;
221 double resoZ_;
222 size_type width_;
223 size_type height_;
224 size_type depth_;
225
226 public:
227 size_type width( ) const { return( width_ ); }
228 size_type height( ) const { return( height_ ); }
229 size_type depth( ) const { return( depth_ ); }
230
232 bool operator ()( difference_type i, difference_type j, difference_type k ) const
233 {
234 return( i * i * resoX_ * resoX_ + j * j * resoY_ * resoY_+ k * k * resoZ_ * resoZ_ <= radius_ * radius_ );
235 }
236
238 expand_mode_type expand_mode( ) const { return( NC6 ); }
239
241 sphere( double radius, double resoX = 1.0, double resoY = 1.0, double resoZ = 1.0, bool radiusInPhysicalCoords = false ) : radius_( radius )
242 {
243 if (radiusInPhysicalCoords)
244 {
245 resoX_ = resoX;
246 resoY_ = resoY;
247 resoZ_ = resoZ;
248 }
249 else
250 {
251 double max_reso = resoX > resoY ? resoX: resoY;
252 max_reso = max_reso > resoZ ? max_reso : resoZ;
253
254 resoX_ = resoX / max_reso;
255 resoY_ = resoY / max_reso;
256 resoZ_ = resoZ / max_reso;
257 }
258 width_ = static_cast< size_type >( std::ceil( radius_ / resoX_ ) ) * 2 + 1;
259 height_ = static_cast< size_type >( std::ceil( radius_ / resoY_ ) ) * 2 + 1;
260 depth_ = static_cast< size_type >( std::ceil( radius_ / resoZ_ ) ) * 2 + 1;
261 }
262 };
263
264
269 template < class T >
270 class less
271 {
272 private:
273 typedef size_t size_type;
274 typedef ptrdiff_t difference_type;
275 typedef T value_type;
276
277 protected:
278 value_type th_;
279
280 public:
282 template < class VerifyList >
283 bool operator()( const VerifyList &elements, size_type num ) const
284 {
285 typedef typename VerifyList::value_type verify_value_type;
286
287 for( size_type i = 0 ; i < num ; i++ )
288 {
289 const verify_value_type &v = elements[ i ];
290 if( !( v < th_ ) )
291 {
292 return ( false );
293 }
294 }
295
296 return ( true );
297 }
298
300 bool require_all_elements( ) const { return( false ); }
301
303 less( value_type threshold ) : th_( threshold ){ }
304 };
305
306
311 template < class T >
312 class greater
313 {
314 private:
315 typedef size_t size_type;
316 typedef ptrdiff_t difference_type;
317 typedef T value_type;
318
319 protected:
320 value_type th_;
321
322 public:
324 template < class VerifyList >
325 bool operator()( const VerifyList &elements, size_type num ) const
326 {
327 typedef typename VerifyList::value_type verify_value_type;
328
329 for( size_type i = 0 ; i < num ; i++ )
330 {
331 const verify_value_type &v = elements[ i ];
332 if( v < th_ )
333 {
334 return ( false );
335 }
336 }
337
338 return ( true );
339 }
340
342 bool require_all_elements( ) const { return( false ); }
343
345 greater( value_type threshold ) : th_( threshold ){ }
346 };
347
348
353 template < class T >
354 class equal
355 {
356 private:
357 typedef size_t size_type;
358 typedef ptrdiff_t difference_type;
359 typedef T value_type;
360
361 protected:
362 value_type th_;
363
364 public:
366 template < class VerifyList >
367 bool operator()( const VerifyList &elements, size_type num ) const
368 {
369 typedef typename VerifyList::value_type verify_value_type;
370
371 for( size_type i = 0 ; i < num ; i++ )
372 {
373 const verify_value_type &v = elements[ i ];
374 if( v < th_ || th_ < v )
375 {
376 return ( false );
377 }
378 }
379
380 return ( true );
381 }
382
384 bool require_all_elements( ) const { return( false ); }
385
387 equal( value_type threshold ) : th_( threshold ){ }
388 };
389
390
395 template < class T >
396 class range
397 {
398 private:
399 typedef size_t size_type;
400 typedef ptrdiff_t difference_type;
401 typedef T value_type;
402
403 protected:
404 value_type min_;
405 value_type max_;
406
407 public:
409 template < class VerifyList >
410 bool operator()( const VerifyList &elements, size_type num ) const
411 {
412 typedef typename VerifyList::value_type verify_value_type;
413
414 for( size_type i = 0 ; i < num ; i++ )
415 {
416 const verify_value_type &v = elements[ i ];
417 if( v < min_ || max_ < v )
418 {
419 return ( false );
420 }
421 }
422
423 return ( true );
424 }
425
427 bool require_all_elements( ) const { return( false ); }
428
430 range( value_type min, value_type max ) : min_( min ), max_( max ){ }
431 };
432 }
433
434
435 // 領域拡張法の内部で利用する型など
436 namespace __region_growing_utility__
437 {
438 struct pointer_diff2
439 {
440 typedef ptrdiff_t difference_type;
441
442 ptrdiff_t diff1;
443 ptrdiff_t diff2;
444
445 pointer_diff2( ptrdiff_t d1 = 0, ptrdiff_t d2 = 0 ) : diff1( d1 ), diff2( d2 ){ }
446 };
447
448 typedef vector3< ptrdiff_t > point_type;
449 typedef std::vector< point_type > point_list_type;
450 typedef std::vector< pointer_diff2 > ptrdiff_list_type;
451
452
453 struct no_mask
454 {
455 typedef size_t size_type;
456 typedef ptrdiff_t difference_type;
457
458 bool empty( ) const { return( true ); }
459 bool operator ()( difference_type /* i */, difference_type /* j */, difference_type /* k */ ) const { return( false ); }
460 };
461
462
463 template < class Array1, class Array2, class Component >
464 ptrdiff_list_type create_component_list( const Array1 &in1, const Array2 &in2, const Component &components )
465 {
466 typedef typename Array1::const_pointer const_pointer1;
467 typedef typename Array2::const_pointer const_pointer2;
468 typedef typename ptrdiff_list_type::size_type size_type;
469 typedef typename ptrdiff_list_type::difference_type difference_type;
470 typedef typename ptrdiff_list_type::value_type value_type;
471
472 ptrdiff_list_type list;
473
474 difference_type w = components.width( );
475 difference_type h = components.height( );
476 difference_type d = components.depth( );
477 difference_type cx = w / 2;
478 difference_type cy = h / 2;
479 difference_type cz = d / 2;
480 difference_type ox1 = in1.width( ) / 2;
481 difference_type oy1 = in1.height( ) / 2;
482 difference_type oz1 = in1.depth( ) / 2;
483 difference_type ox2 = in2.width( ) / 2;
484 difference_type oy2 = in2.height( ) / 2;
485 difference_type oz2 = in2.depth( ) / 2;
486
487 const_pointer1 p1 = &in1( ox1, oy1, oz1 );
488 const_pointer2 p2 = &in2( ox2, oy2, oz2 );
489
490 for( difference_type k = 0 ; k < d ; k++ )
491 {
492 for( difference_type j = 0 ; j < h ; j++ )
493 {
494 for( difference_type i = 0 ; i < w ; i++ )
495 {
496 if( components( i - cx, j - cy, k - cz ) )
497 {
498 ptrdiff_t d1 = &in1( ox1 + i - cx, oy1 + j - cy, oz1 + k - cz ) - p1;
499 ptrdiff_t d2 = &in2( ox2 + i - cx, oy2 + j - cy, oz2 + k - cz ) - p2;
500 list.push_back( pointer_diff2( d1, d2 ) );
501 }
502 }
503 }
504 }
505
506 return( list );
507 }
508
509
510 template < class Array1, class Array2, class Component >
511 ptrdiff_list_type create_update_list( const Array1 &in1, const Array2 &in2, const Component &components )
512 {
513 typedef typename Array1::const_pointer const_pointer1;
514 typedef typename Array2::const_pointer const_pointer2;
515 typedef typename ptrdiff_list_type::size_type size_type;
516 typedef typename ptrdiff_list_type::difference_type difference_type;
517 typedef typename ptrdiff_list_type::value_type value_type;
518
519 ptrdiff_list_type list;
520
521 difference_type ox1 = in1.width( ) / 2;
522 difference_type oy1 = in1.height( ) / 2;
523 difference_type oz1 = in1.depth( ) / 2;
524 difference_type ox2 = in2.width( ) / 2;
525 difference_type oy2 = in2.height( ) / 2;
526 difference_type oz2 = in2.depth( ) / 2;
527
528 const_pointer1 p1 = &in1( ox1, oy1, oz1 );
529 const_pointer2 p2 = &in2( ox2, oy2, oz2 );
530
531 switch( components.expand_mode( ) )
532 {
533 case region_growing_utility::NC4:
534 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1, oz1 ) - p1, &in2( ox2 - 1, oy2, oz2 ) - p2 ) );
535 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1, oz1 ) - p1, &in2( ox2 + 1, oy2, oz2 ) - p2 ) );
536 list.push_back( pointer_diff2( &in1( ox1, oy1 - 1, oz1 ) - p1, &in2( ox2, oy2 - 1, oz2 ) - p2 ) );
537 list.push_back( pointer_diff2( &in1( ox1, oy1 + 1, oz1 ) - p1, &in2( ox2, oy2 + 1, oz2 ) - p2 ) );
538 break;
539
540 case region_growing_utility::NC8:
541 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 - 1, oz1 ) - p1, &in2( ox2 - 1, oy2 - 1, oz2 ) - p2 ) );
542 list.push_back( pointer_diff2( &in1( ox1 , oy1 - 1, oz1 ) - p1, &in2( ox2 , oy2 - 1, oz2 ) - p2 ) );
543 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 - 1, oz1 ) - p1, &in2( ox2 + 1, oy2 - 1, oz2 ) - p2 ) );
544 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 , oz1 ) - p1, &in2( ox2 - 1, oy2 , oz2 ) - p2 ) );
545 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 , oz1 ) - p1, &in2( ox2 + 1, oy2 , oz2 ) - p2 ) );
546 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 + 1, oz1 ) - p1, &in2( ox2 - 1, oy2 + 1, oz2 ) - p2 ) );
547 list.push_back( pointer_diff2( &in1( ox1 , oy1 + 1, oz1 ) - p1, &in2( ox2 , oy2 + 1, oz2 ) - p2 ) );
548 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 + 1, oz1 ) - p1, &in2( ox2 + 1, oy2 + 1, oz2 ) - p2 ) );
549 break;
550
551 case region_growing_utility::NC6:
552 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1, oz1 ) - p1, &in2( ox2 - 1, oy2, oz2 ) - p2 ) );
553 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1, oz1 ) - p1, &in2( ox2 + 1, oy2, oz2 ) - p2 ) );
554 list.push_back( pointer_diff2( &in1( ox1, oy1 - 1, oz1 ) - p1, &in2( ox2, oy2 - 1, oz2 ) - p2 ) );
555 list.push_back( pointer_diff2( &in1( ox1, oy1 + 1, oz1 ) - p1, &in2( ox2, oy2 + 1, oz2 ) - p2 ) );
556 list.push_back( pointer_diff2( &in1( ox1, oy1, oz1 - 1 ) - p1, &in2( ox2, oy2, oz2 - 1 ) - p2 ) );
557 list.push_back( pointer_diff2( &in1( ox1, oy1, oz1 + 1 ) - p1, &in2( ox2, oy2, oz2 + 1 ) - p2 ) );
558 break;
559
560 case region_growing_utility::NC18:
561 list.push_back( pointer_diff2( &in1( ox1 , oy1 - 1, oz1 - 1 ) - p1, &in2( ox2 , oy2 - 1, oz2 - 1 ) - p2 ) );
562 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 , oz1 - 1 ) - p1, &in2( ox2 - 1, oy2 , oz2 - 1 ) - p2 ) );
563 list.push_back( pointer_diff2( &in1( ox1 , oy1 , oz1 - 1 ) - p1, &in2( ox2 , oy2 , oz2 - 1 ) - p2 ) );
564 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 , oz1 - 1 ) - p1, &in2( ox2 + 1, oy2 , oz2 - 1 ) - p2 ) );
565 list.push_back( pointer_diff2( &in1( ox1 , oy1 + 1, oz1 - 1 ) - p1, &in2( ox2 , oy2 + 1, oz2 - 1 ) - p2 ) );
566 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 - 1, oz1 ) - p1, &in2( ox2 - 1, oy2 - 1, oz2 ) - p2 ) );
567 list.push_back( pointer_diff2( &in1( ox1 , oy1 - 1, oz1 ) - p1, &in2( ox2 , oy2 - 1, oz2 ) - p2 ) );
568 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 - 1, oz1 ) - p1, &in2( ox2 + 1, oy2 - 1, oz2 ) - p2 ) );
569 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 , oz1 ) - p1, &in2( ox2 - 1, oy2 , oz2 ) - p2 ) );
570 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 , oz1 ) - p1, &in2( ox2 + 1, oy2 , oz2 ) - p2 ) );
571 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 + 1, oz1 ) - p1, &in2( ox2 - 1, oy2 + 1, oz2 ) - p2 ) );
572 list.push_back( pointer_diff2( &in1( ox1 , oy1 + 1, oz1 ) - p1, &in2( ox2 , oy2 + 1, oz2 ) - p2 ) );
573 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 + 1, oz1 ) - p1, &in2( ox2 + 1, oy2 + 1, oz2 ) - p2 ) );
574 list.push_back( pointer_diff2( &in1( ox1 , oy1 - 1, oz1 + 1 ) - p1, &in2( ox2 , oy2 - 1, oz2 + 1 ) - p2 ) );
575 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 , oz1 + 1 ) - p1, &in2( ox2 - 1, oy2 , oz2 + 1 ) - p2 ) );
576 list.push_back( pointer_diff2( &in1( ox1 , oy1 , oz1 + 1 ) - p1, &in2( ox2 , oy2 , oz2 + 1 ) - p2 ) );
577 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 , oz1 + 1 ) - p1, &in2( ox2 + 1, oy2 , oz2 + 1 ) - p2 ) );
578 list.push_back( pointer_diff2( &in1( ox1 , oy1 + 1, oz1 + 1 ) - p1, &in2( ox2 , oy2 + 1, oz2 + 1 ) - p2 ) );
579 break;
580
581 case region_growing_utility::NC26:
582 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 - 1, oz1 - 1 ) - p1, &in2( ox2 - 1, oy2 - 1, oz2 - 1 ) - p2 ) );
583 list.push_back( pointer_diff2( &in1( ox1 , oy1 - 1, oz1 - 1 ) - p1, &in2( ox2 , oy2 - 1, oz2 - 1 ) - p2 ) );
584 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 - 1, oz1 - 1 ) - p1, &in2( ox2 + 1, oy2 - 1, oz2 - 1 ) - p2 ) );
585 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 , oz1 - 1 ) - p1, &in2( ox2 - 1, oy2 , oz2 - 1 ) - p2 ) );
586 list.push_back( pointer_diff2( &in1( ox1 , oy1 , oz1 - 1 ) - p1, &in2( ox2 , oy2 , oz2 - 1 ) - p2 ) );
587 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 , oz1 - 1 ) - p1, &in2( ox2 + 1, oy2 , oz2 - 1 ) - p2 ) );
588 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 + 1, oz1 - 1 ) - p1, &in2( ox2 - 1, oy2 + 1, oz2 - 1 ) - p2 ) );
589 list.push_back( pointer_diff2( &in1( ox1 , oy1 + 1, oz1 - 1 ) - p1, &in2( ox2 , oy2 + 1, oz2 - 1 ) - p2 ) );
590 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 + 1, oz1 - 1 ) - p1, &in2( ox2 + 1, oy2 + 1, oz2 - 1 ) - p2 ) );
591 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 - 1, oz1 ) - p1, &in2( ox2 - 1, oy2 - 1, oz2 ) - p2 ) );
592 list.push_back( pointer_diff2( &in1( ox1 , oy1 - 1, oz1 ) - p1, &in2( ox2 , oy2 - 1, oz2 ) - p2 ) );
593 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 - 1, oz1 ) - p1, &in2( ox2 + 1, oy2 - 1, oz2 ) - p2 ) );
594 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 , oz1 ) - p1, &in2( ox2 - 1, oy2 , oz2 ) - p2 ) );
595 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 , oz1 ) - p1, &in2( ox2 + 1, oy2 , oz2 ) - p2 ) );
596 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 + 1, oz1 ) - p1, &in2( ox2 - 1, oy2 + 1, oz2 ) - p2 ) );
597 list.push_back( pointer_diff2( &in1( ox1 , oy1 + 1, oz1 ) - p1, &in2( ox2 , oy2 + 1, oz2 ) - p2 ) );
598 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 + 1, oz1 ) - p1, &in2( ox2 + 1, oy2 + 1, oz2 ) - p2 ) );
599 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 - 1, oz1 + 1 ) - p1, &in2( ox2 - 1, oy2 - 1, oz2 + 1 ) - p2 ) );
600 list.push_back( pointer_diff2( &in1( ox1 , oy1 - 1, oz1 + 1 ) - p1, &in2( ox2 , oy2 - 1, oz2 + 1 ) - p2 ) );
601 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 - 1, oz1 + 1 ) - p1, &in2( ox2 + 1, oy2 - 1, oz2 + 1 ) - p2 ) );
602 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 , oz1 + 1 ) - p1, &in2( ox2 - 1, oy2 , oz2 + 1 ) - p2 ) );
603 list.push_back( pointer_diff2( &in1( ox1 , oy1 , oz1 + 1 ) - p1, &in2( ox2 , oy2 , oz2 + 1 ) - p2 ) );
604 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 , oz1 + 1 ) - p1, &in2( ox2 + 1, oy2 , oz2 + 1 ) - p2 ) );
605 list.push_back( pointer_diff2( &in1( ox1 - 1, oy1 + 1, oz1 + 1 ) - p1, &in2( ox2 - 1, oy2 + 1, oz2 + 1 ) - p2 ) );
606 list.push_back( pointer_diff2( &in1( ox1 , oy1 + 1, oz1 + 1 ) - p1, &in2( ox2 , oy2 + 1, oz2 + 1 ) - p2 ) );
607 list.push_back( pointer_diff2( &in1( ox1 + 1, oy1 + 1, oz1 + 1 ) - p1, &in2( ox2 + 1, oy2 + 1, oz2 + 1 ) - p2 ) );
608 break;
609
610 case region_growing_utility::ALL:
611 default:
612 list = create_component_list( in1, in2, components );
613 break;
614 }
615
616 return( list );
617 }
618
619
620 template < class T, class Allocator, class Array, class PointType >
621 bool check_component_list( const array2< T, Allocator > &in1, const Array &in2, const PointType &pt, pointer_diff2 &ptr )
622 {
623 typedef typename array2< T, Allocator >::const_pointer const_pointer1;
624 typedef typename Array::const_pointer const_pointer2;
625 typedef typename ptrdiff_list_type::difference_type difference_type;
626
627 difference_type x = static_cast< difference_type >( pt.x );
628 difference_type y = static_cast< difference_type >( pt.y );
629 difference_type w = in1.width( );
630 difference_type h = in1.height( );
631
632 if( x < 0 || x >= w || y < 0 || y >= h )
633 {
634 return( false );
635 }
636 else if( in2( x, y ) != 0 )
637 {
638 return( false );
639 }
640 else
641 {
642 ptr.diff1 = &in1( x, y) - &in1[ 0 ];
643 ptr.diff2 = &in2( x, y ) - &in2[ 0 ];
644 return( true );
645 }
646 }
647
648
649 template < class T, class Allocator, class Array, class PointType >
650 bool check_component_list( const array3< T, Allocator > &in1, const Array &in2, const PointType &pt, pointer_diff2 &ptr )
651 {
652 typedef typename array3< T, Allocator >::const_pointer const_pointer1;
653 typedef typename Array::const_pointer const_pointer2;
654 typedef typename ptrdiff_list_type::difference_type difference_type;
655
656 difference_type x = static_cast< difference_type >( pt.x );
657 difference_type y = static_cast< difference_type >( pt.y );
658 difference_type z = static_cast< difference_type >( pt.z );
659 difference_type w = in1.width( );
660 difference_type h = in1.height( );
661 difference_type d = in1.depth( );
662
663 if( x < 0 || x >= w || y < 0 || y >= h || z < 0 || z >= d )
664 {
665 return( false );
666 }
667 else if( in2( x, y, z ) != 0 )
668 {
669 return( false );
670 }
671 else
672 {
673 ptr.diff1 = &in1( x, y, z ) - &in1[ 0 ];
674 ptr.diff2 = &in2( x, y, z ) - &in2[ 0 ];
675 return( true );
676 }
677 }
678
679
680 template < class PointType >
681 struct point_list_converter
682 {
683 typedef std::vector< PointType > point_list_type;
684
685 static point_list_type create_point_list( const PointType &pt )
686 {
687 point_list_type list( 1 );
688 list[ 0 ] = pt;
689 return( list );
690 }
691 };
692
693
694 template < class PointType >
695 struct point_list_converter< std::vector< PointType > >
696 {
697 typedef std::vector< PointType > point_list_type;
698
699 static point_list_type create_point_list( const point_list_type &list )
700 {
701 return( list );
702 }
703 };
704
705 template < class PointType >
706 struct point_list_converter< std::deque< PointType > >
707 {
708 typedef std::deque< PointType > point_list_type;
709
710 static point_list_type create_point_list( const point_list_type &list )
711 {
712 return( list );
713 }
714 };
715
716
717 template < class T >
718 inline const T maximum( const T &v0, const typename type_trait< T >::value_type &v1 )
719 {
720 return( v0 > v1 ? v0 : v1 );
721 }
722
723 template < class T >
724 inline const T maximum( const T &v0, const typename type_trait< T >::value_type &v1, const typename type_trait< T >::value_type &v2 )
725 {
726 return( v0 > v1 ? ( v0 > v2 ? v0 : v2 ) : ( v1 > v2 ? v1 : v2 ) );
727 }
728
729 template < class T >
730 inline const T maximum( const T &v0, const typename type_trait< T >::value_type &v1, const typename type_trait< T >::value_type &v2, const typename type_trait< T >::value_type &v3 )
731 {
732 return( maximum( maximum( v0, v1 ), maximum( v2, v3 ) ) );
733 }
734 }
735
736
752 template < class Array1, class Array2, class MaskType, class PointList, class Component, class Condition >
753 typename Array1::difference_type region_growing( const Array1 &in, Array2 &out, const MaskType &mask, const PointList &start_points, typename Array2::value_type output_value,
754 const Component &components, const Condition &condition, typename Array1::size_type max_paint )
755 {
756 if( in.empty( ) || is_same_object( in, out ) )
757 {
758 return( -1 );
759 }
760
761 typedef typename Array1::template rebind< unsigned char >::other mask_type;
762 typedef __region_growing_utility__::pointer_diff2 pointer_diff_type;
763
764 typedef typename Array1::size_type size_type;
765 typedef typename Array1::value_type value_type;
766 typedef typename Array1::difference_type difference_type;
767 typedef typename Array2::value_type out_value_type;
768
769 typedef typename Array1::const_pointer const_pointer;
770 typedef typename mask_type::pointer work_pointer;
771 typedef typename Array2::pointer output_pointer;
772
773 size_type rx = components.width( ) / 2;
774 size_type ry = components.height( ) / 2;
775 size_type rz = components.depth( ) / 2;
776
777 // 注目点が範囲外に出ないことを監視するマスク
778 // 最低でも外側に１画素の安全領域を作成する
779 marray< mask_type > work( in, __region_growing_utility__::maximum( rx, ry, rz, 1 ) );
780
781 __region_growing_utility__::ptrdiff_list_type clist = __region_growing_utility__::create_component_list( in, work, components ); // 構造要素内の画素のリスト
782 __region_growing_utility__::ptrdiff_list_type ulist = __region_growing_utility__::create_update_list( in, work, components ); // ある注目点の次に注目点になる画素のリスト
783
784 std::deque< pointer_diff_type > que; // 注目点画素のリスト
785 std::vector< value_type > element( clist.size( ) ); // 拡張条件の判定に用いる画素の配列
786 __region_growing_utility__::ptrdiff_list_type elist( clist.size( ) ); // 拡張条件の判定に用いる画素へのポインタを覚える配列
787
788 // 画像の外側に拡張しないようにマスクを設定する
789 work.fill( );
790 work.fill_margin( 255 );
791
792 // マスクの内容を注目点のデータに反映させる
793 if( !mask.empty( ) )
794 {
795 for( size_type k = 0 ; k < in.depth( ) ; k++ )
796 {
797 for( size_type j = 0 ; j < in.height( ) ; j++ )
798 {
799 for( size_type i = 0 ; i < in.width( ) ; i++ )
800 {
801 if( mask( i, j, k ) )
802 {
803 work( i, j, k ) = 255;
804 }
805 }
806 }
807 }
808 }
809
810 typedef __region_growing_utility__::point_list_converter< PointList > start_point_list_converter;
811 typedef typename start_point_list_converter::point_list_type start_point_list_type;
812 start_point_list_type sps = start_point_list_converter::create_point_list( start_points );
813
814 // 拡張開始点が画像内に存在しているかどうかをチェックする
815 for( typename start_point_list_type::const_iterator ite = sps.begin( ) ; ite != sps.end( ) ; ++ite )
816 {
817 __region_growing_utility__::pointer_diff2 ptr;
818 if( __region_growing_utility__::check_component_list( in, work, *ite, ptr ) )
819 {
820 que.push_back( ptr );
821 }
822 }
823
824 // 出力画像の大きさをチェックする
825 if( in.size( ) != out.size( ) )
826 {
827 out.resize( in.width( ), in.height( ), in.depth( ) );
828 out.reso1( in.reso1( ) );
829 out.reso2( in.reso2( ) );
830 out.reso3( in.reso3( ) );
831 }
832
833
834 size_type num_painted = 0; // 塗りつぶした画素数のカウンタ
835
836
837 // 領域拡張の条件によって，構造要素内の全要素を判定に用いるかどうかを決める
838 // require_all_elements が false の場合には，すでに判定済みの画素は2度目の判定を行わない
839 if( condition.require_all_elements( ) )
840 {
841 // 構造要素内のすべての点を用いて評価を行う
842 while( !que.empty( ) )
843 {
844 if( num_painted >= max_paint )
845 {
846 return( num_painted );
847 }
848
849 // リストの先頭画素を取り出す
850 pointer_diff_type cur = que.front( );
851 que.pop_front( );
852
853 const_pointer pi = &in[ 0 ] + cur.diff1;
854 work_pointer pw = &work[ 0 ] + cur.diff2;
855 output_pointer po = &out[ 0 ] + cur.diff1;
856
857 // 以降の処理で注目点とならないように，現在の注目点をマスクする
858 pw[ 0 ] |= 0xf0;
859
860 // 拡張条件判定に用いる画素を列挙する
861 // ここで，マスクされている領域は範囲に含めない（画像外など）
862 size_type num = 0;
863 for( size_type i = 0 ; i < clist.size( ) ; i++ )
864 {
865 const pointer_diff_type &d = clist[ i ];
866 if( ( pw[ d.diff2 ] & 0x0f ) == 0 )
867 {
868 element[ num ] = pi[ d.diff1 ];
869 elist[ num ] = d;
870 num++;
871 }
872 }
873
874 // 拡張条件の判定を行う
875 if( num != 0 && condition( element, num ) )
876 {
877 // 条件を満たした構造要素内の画素すべてを塗りつぶす
878 for( size_type i = 0 ; i < num ; i++ )
879 {
880 const pointer_diff_type &d = elist[ i ];
881 num_painted += po[ d.diff1 ] != output_value ? 1 : 0;
882 po[ d.diff1 ] = output_value;
883 }
884
885 // 構造要素によって次の注目点を決定し，リストの最後に追加する
886 for( size_type i = 0 ; i < ulist.size( ) ; i++ )
887 {
888 const pointer_diff_type &d = ulist[ i ];
889 if( ( pw[ d.diff2 ] & 0xf0 ) == 0 )
890 {
891 pw[ d.diff2 ] |= 0xf0;
892 que.push_back( pointer_diff_type( cur.diff1 + d.diff1, cur.diff2 + d.diff2 ) );
893 }
894 }
895 }
896 }
897 }
898 else
899 {
900 // すでに判定された画素については評価を行わない
901 while( !que.empty( ) )
902 {
903 if( num_painted >= max_paint )
904 {
905 return( num_painted );
906 }
907
908 // リストの先頭画素を取り出す
909 pointer_diff_type cur = que.front( );
910 que.pop_front( );
911
912 const_pointer pi = &in[ 0 ] + cur.diff1;
913 work_pointer pw = &work[ 0 ] + cur.diff2;
914 output_pointer po = &out[ 0 ] + cur.diff1;
915
916 // 以降の処理で注目点とならないように，現在の注目点をマスクする
917 pw[ 0 ] |= 0xf0;
918
919 // 拡張条件判定に用いる画素を列挙する
920 // ここで，マスクされている領域は範囲に含めない（画像外など）
921 // また，すでに判定済みの画素も除外する
922 // 判定画素と更新部分のリストを作成する
923 size_type num = 0;
924 for( size_type i = 0 ; i < clist.size( ) ; i++ )
925 {
926 const pointer_diff_type &d = clist[ i ];
927 if( ( pw[ d.diff2 ] & 0x0f ) == 0 )
928 {
929 element[ num ] = pi[ d.diff1 ];
930 elist[ num ] = d;
931 num++;
932 }
933 }
934
935 // 拡張条件の判定を行う
936 if( condition( element, num ) )
937 {
938 // 条件を満たした構造要素内の画素すべてを塗りつぶす
939 for( size_type i = 0 ; i < num ; i++ )
940 {
941 const pointer_diff_type &d = elist[ i ];
942 num_painted += po[ d.diff1 ] != output_value ? 1 : 0;
943 pw[ d.diff2 ] |= 1;
944 po[ d.diff1 ] = output_value;
945 }
946
947 // 構造要素によって次の注目点を決定し，リストの最後に追加する
948 for( size_type i = 0 ; i < ulist.size( ) ; i++ )
949 {
950 const pointer_diff_type &d = ulist[ i ];
951 if( ( pw[ d.diff2 ] & 0xf0 ) == 0 )
952 {
953 pw[ d.diff2 ] |= 0xf0;
954 que.push_back( pointer_diff_type( cur.diff1 + d.diff1, cur.diff2 + d.diff2 ) );
955 }
956 }
957 }
958 }
959 }
960
961 return( num_painted );
962 }
963
964
965
982 template < class Array1, class Array2, class PointList, class Component, class Condition >
983 typename Array1::difference_type region_growing( const Array1 &in, Array2 &out, const PointList &start_points, typename Array2::value_type output_value,
984 const Component &components, const Condition &condition, typename Array1::size_type max_paint = type_limits< typename Array1::size_type >::maximum( ) )
985 {
986 return( region_growing( in, out, __region_growing_utility__::no_mask( ), start_points, output_value, components, condition, max_paint ) );
987 }
988
989
991 // 領域拡張法グループの終わり
992
993
994 // mist名前空間の終わり
995 _MIST_END
996
997
998 #endif // __INCLUDE_MIST_REGIONGROWING__
999