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