thinning.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 
39 
40 #ifndef __INCLUDE_MIST_THINNING__
41 #define __INCLUDE_MIST_THINNING__
42 
43 
44 #ifndef __INCLUDE_MIST_H__
45 #include "../mist.h"
46 #endif
47 
48 #ifndef __INCLUDE_MIST_LIMITS__
49 #include "../limits.h"
50 #endif
51 
52 #ifndef __INCLUDE_MIST_MATRIX__
53 #include "../matrix.h"
54 #endif
55 
56 #ifndef __INCLUDE_MIST_DISTANCE_TRANSFORM__
57 #include "distance.h"
58 #endif
59 
60 #include <vector>
61 
62 
63 // mist名前空間の始まり
64 _MIST_BEGIN
65 
66 
67 namespace __thinning_controller__
68 {
69  template < class T, class Allocator >
70  inline void val9( const array2< T, Allocator > &in, typename array2< T, Allocator >::value_type p[ 9 ],
71  typename array2< T, Allocator >::size_type i, typename array2< T, Allocator >::size_type j )
72  {
73  p[ 0 ] = in( i , j );
74  p[ 1 ] = in( i + 1, j );
75  p[ 2 ] = in( i + 1, j - 1 );
76  p[ 3 ] = in( i , j - 1 );
77  p[ 4 ] = in( i - 1, j - 1 );
78  p[ 5 ] = in( i - 1, j );
79  p[ 6 ] = in( i - 1, j + 1 );
80  p[ 7 ] = in( i , j + 1 );
81  p[ 8 ] = in( i + 1, j + 1 );
82  }
83 
84  template < class T > inline T pixel( T p )
85  {
86  return ( ( p > 0 ) ? 1 : 0 );
87  }
88 
89  template < class T > inline T pixel( T p[ 9 ], int index )
90  {
91  if( index <= 0 )
92  {
93  return ( pixel( p[ index + 8 ] ) );
94  }
95  else if( index >= 9 )
96  {
97  return ( pixel( p[ index - 8 ] ) );
98  }
99  else
100  {
101  return ( pixel( p[ index ] ) );
102  }
103  }
104 
105  template < class T > inline T pixel_( T p[ 9 ], int index )
106  {
107  return ( 1 - pixel( p, index ) );
108  }
109 
110  template < class T > inline int Nc4( T p[ 9 ] )
111  {
112  T ret = 0;
113 
114  for( int i = 1 ; i < 9 ; i += 2 )
115  {
116  ret += pixel( p, i ) - pixel( p, i ) * pixel( p, i + 1 ) * pixel( p, i + 2 );
117  }
118 
119  return ( ret );
120  }
121 
122  template < class T > inline int Nc8( T p[ 9 ] )
123  {
124  T ret = 0;
125 
126  for( int i = 1 ; i < 9 ; i += 2 )
127  {
128  ret += pixel_( p, i ) - pixel_( p, i ) * pixel_( p, i + 1 ) * pixel_( p, i + 2 );
129  }
130 
131  return ( ret );
132  }
133 
134  // C. J. Hilditch, ``Linear Skeleton from Square Cupboards,'' In: Machine Intelligence 6, B. Meltzer and D. Michie eds., Edinburgh Univ. Press, pp.403?420, 1969
135  // 鳥脇純一郎, ``画像理解のためのディジタル画像処理〔II〕，'' 昭晃堂，pp.56-59，1988
136  template < class T, class Allocator >
137  void thinning( array2< T, Allocator > &ia )
138  {
139  typedef typename array2< T, Allocator >::size_type size_type;
140  typedef typename array2< T, Allocator >::value_type value_type;
141 
142  bool flag;
143  size_type w = ia.width( );
144  size_type h = ia.height( );
145  size_type i, j;
146  char p[ 9 ], q[ 9 ], value;
147  array2< char > ic( w, h );
148  array2< char > id( w, h );
149 
150  for( i = 0 ; i < ia.size( ) ; i++ )
151  {
152  id[ i ] = ia[ i ] == 0 ? 0 : 1;
153  }
154 
155  flag = true;
156 
157  while( flag )
158  {
159  flag = false;
160  ic = id;
161 
162  for( j = 1 ; j < h - 1 ; j++ )
163  {
164  for( i = 1 ; i < w - 1 ; i++ )
165  {
166  val9( ic, p, i, j );
167  val9( id, q, i, j );
168 
169  if( pixel( p[ 0 ] ) != 1 )
170  {
171  continue;
172  }
173 
174  // 条件1 = 入力画像Fにおいて，x0の4近傍に少なくとも1個の0画素が存在する
175  if( p[ 1 ] == 1 && p[ 3 ] == 1 && p[ 5 ] == 1 && p[ 7 ] == 1 )
176  {
177  continue;
178  }
179 
180  // 条件2 = 入力画像Fにおいて，x0の8近傍に2個以上の1画素がある(x0が端点でない)
181  value = pixel( p[ 1 ] ) + pixel( p[ 2 ] ) + pixel( p[ 3 ] ) + pixel( p[ 4 ] ) + pixel( p[ 5 ] ) + pixel( p[ 6 ] ) + pixel( p[ 7 ] ) + pixel( p[ 8 ] );
182 
183  if( value <= 1 )
184  {
185  continue;
186  }
187 
188  // 条件3 = 作業画像F'において，-1の画素も0画素とみなした時，x0が孤立点ではない
189  value = pixel( q[ 1 ] ) + pixel( q[ 2 ] ) + pixel( q[ 3 ] ) + pixel( q[ 4 ] ) + pixel( q[ 5 ] ) + pixel( q[ 6 ] ) + pixel( q[ 7 ] ) + pixel( q[ 8 ] );
190 
191  if( value == 0 )
192  {
193  continue;
194  }
195 
196  // 条件4 = 入力画像Fにおいて，x0におけるNc[8]の連結数が1である
197  if( Nc8( p ) != 1 )
198  {
199  continue;
200  }
201 
202  // 条件5 = i=3,5に対して次の事柄が成り立つ
203  // xiがちょうどこの反復で1から0に変えられてはいない(xiが-1になっていない)か，
204  // または，それが-1になっていても，-1を全て0とみなして計算したx0の連結数Nc[8]が1に等しい
205  if( q[ 3 ] == -1 || q[ 5 ] == -1 )
206  {
207  if( Nc8( q ) != 1 )
208  {
209  continue;
210  }
211  }
212 
213  // 全ての条件を満たした画素は-1にする
214  id( i, j ) = -1;
215 
216  flag = true;
217  }
218  }
219 
220  for( i = 0 ; i < id.size( ) ; i++ )
221  {
222  id[ i ] = id[ i ] == -1 ? 0 : id[ i ];
223  }
224  }
225 
226  for( i = 0 ; i < ia.size( ) ; i++ )
227  {
228  ia[ i ] = id[ i ] > 0 ? 1 : 0;
229  }
230  }
231 }
232 
233 
234 
235 namespace __euclidean_utility__
236 {
237  template < class T, class Allocator >
238  inline void val9_1( const array2< T, Allocator > &in, typename array2< T, Allocator >::value_type p[ 9 ],
239  typename array2< T, Allocator >::size_type i, typename array2< T, Allocator >::size_type j )
240  {
241  p[ 0 ] = in( i , j );
242  p[ 1 ] = in( i + 1, j );
243  p[ 2 ] = in( i + 1, j - 1 );
244  p[ 3 ] = in( i , j - 1 );
245  p[ 4 ] = in( i - 1, j - 1 );
246  p[ 5 ] = in( i - 1, j );
247  p[ 6 ] = in( i - 1, j + 1 );
248  p[ 7 ] = in( i , j + 1 );
249  p[ 8 ] = in( i + 1, j + 1 );
250  }
251 
252  template < class T, class Allocator >
253  inline void val9_2( const array2< T, Allocator > &in, typename array2< T, Allocator >::value_type p[ 9 ],
254  typename array2< T, Allocator >::size_type i, typename array2< T, Allocator >::size_type j )
255  {
256  p[ 0 ] = in( i , j );
257  p[ 1 ] = in( i - 1, j );
258  p[ 2 ] = in( i - 1, j + 1 );
259  p[ 3 ] = in( i , j + 1 );
260  p[ 4 ] = in( i + 1, j + 1 );
261  p[ 5 ] = in( i + 1, j );
262  p[ 6 ] = in( i + 1, j - 1 );
263  p[ 7 ] = in( i , j - 1 );
264  p[ 8 ] = in( i - 1, j - 1 );
265  }
266 
267  template < class T, class Allocator >
268  inline void val9_3( const array2< T, Allocator > &in, typename array2< T, Allocator >::value_type p[ 9 ],
269  typename array2< T, Allocator >::size_type i, typename array2< T, Allocator >::size_type j )
270  {
271  p[ 0 ] = in( i , j );
272  p[ 1 ] = in( i + 1, j );
273  p[ 2 ] = in( i + 1, j + 1 );
274  p[ 3 ] = in( i , j + 1 );
275  p[ 4 ] = in( i - 1, j + 1 );
276  p[ 5 ] = in( i - 1, j );
277  p[ 6 ] = in( i - 1, j - 1 );
278  p[ 7 ] = in( i , j - 1 );
279  p[ 8 ] = in( i + 1, j - 1 );
280  }
281 
282  template < class T > inline T max_pixel( T p[ 9 ] )
283  {
284  T max = 0;
285 
286  for( int i = 1 ; i < 9 ; i++ )
287  {
288  if( std::abs( p[ i ] ) > max )
289  {
290  max = std::abs( p[ i ] );
291  }
292  }
293 
294  return ( max );
295  }
296 
297  template < class T > inline T plus_pixel_num( T p[ 9 ] )
298  {
299  int num = 0;
300 
301  for( int i = 1 ; i < 9 ; i++ )
302  {
303  if( p[ i ] > 0 )
304  {
305  num++;
306  }
307  }
308 
309  return ( num );
310  }
311 
312  template < class T > inline T pixel( T p )
313  {
314  return ( ( p > 0 ) ? 1 : 0 );
315  }
316 
317  template < class T > inline T pixel( T p[ 9 ], int index )
318  {
319  if( index <= 0 )
320  {
321  return ( pixel( p[ index + 8 ] ) );
322  }
323  else if( index >= 9 )
324  {
325  return ( pixel( p[ index - 8 ] ) );
326  }
327  else
328  {
329  return ( pixel( p[ index ] ) );
330  }
331  }
332 
333  template < class T > inline T pixel_( T p[ 9 ], int index )
334  {
335  return ( 1 - pixel( p, index ) );
336  }
337 
338  template < class T > inline int Nc4( T p[ 9 ] )
339  {
340  T ret = 0;
341 
342  for( int i = 1 ; i < 9 ; i += 2 )
343  {
344  ret += pixel( p, i ) - pixel( p, i ) * pixel( p, i + 1 ) * pixel( p, i + 2 );
345  }
346 
347  return ( ret );
348  }
349 
350  template < class T > inline T Nc8( T p[ 9 ] )
351  {
352  T ret = 0;
353 
354  for( int i = 1 ; i < 9 ; i += 2 )
355  {
356  ret += pixel_( p, i ) - pixel_( p, i ) * pixel_( p, i + 1 ) * pixel_( p, i + 2 );
357  }
358 
359  return ( ret );
360  }
361 
362  // ユークリッド距離を用いた細線化の実装 by 林さん
363  //
364  // - 参考文献
365  // - 鈴木智, 阿部圭一, ``距離変換の結果を利用した二値画像の逐次細線化,'' 電子情報通信学会論文誌D, vol.68-D, no.4, pp.473-480, 1985.
366  //
367  template < class T, class Allocator >
368  void thinning( array2< T, Allocator > &ia )
369  {
370  typedef typename array2< T, Allocator >::size_type size_type;
371  typedef typename array2< T, Allocator >::value_type value_type;
372 
373  size_type w = ia.width( );
374  size_type h = ia.height( );
375  size_type i, j;
376  double p[ 9 ], q[ 9 ];
377  array2< double > ic( w, h, ia.reso1( ), ia.reso2( ) );
378  array2< double > id( w, h, ia.reso1( ), ia.reso2( ) );
379 
380  for( i = 0 ; i < ia.size( ) ; i++ )
381  {
382  id[ i ] = ia[ i ] == 0 ? 0 : 1;
383  }
384 
385  //距離変換
386  euclidean::distance_transform( id, id );
387 
388  //細め処理　1回目
389  ic = id;
390  for( j = 1 ; j < h - 1 ; j++ )
391  {
392  for( i = 1 ; i < w - 1 ; i++ )
393  {
394  val9_1( ic, p, i, j );
395  val9_1( id, q, i, j );
396 
397  if( q[ 0 ] <= 0 )
398  {
399  continue;
400  }
401 
402  if( p[ 3 ] > 0 && p[ 5 ] > 0 )
403  {
404  continue;
405  }
406 
407  if( Nc8( p ) != 1 )
408  {
409  continue;
410  }
411 
412  if( max_pixel( p ) <= std::abs( q[ 0 ] ) )
413  {
414  continue;
415  }
416 
417  if( ( q[ 0 ] > q[ 1 ] && q[ 1 ] > 0 ) && ( ( p[ 2 ] > 0 ) || ( p[ 3 ] > 0 ) ) && ( ( q[ 7 ] > 0 ) || ( q[ 8 ] > 0 ) ) )
418  {
419  continue;
420  }
421 
422  if( ( q[ 0 ] > q[ 7 ] && q[ 7 ] > 0 ) && ( ( p[ 5 ] > 0 ) || ( q[ 6 ] > 0 ) ) && ( ( q[ 8 ] > 0 ) || ( q[ 1 ] > 0 ) ) )
423  {
424  continue;
425  }
426 
427  ic( i, j ) = - id( i, j );
428  }
429  }
430 
431  //細め処理　2回目
432  id = ic;
433  for( j = h - 2 ; j > 0 ; j-- )
434  {
435  for( i = w - 2 ; i > 0 ; i-- )
436  {
437  val9_2( id, p, i, j );
438  val9_2( ic, q, i, j );
439 
440  if( q[ 0 ] <= 0 )
441  {
442  continue;
443  }
444 
445  if( p[ 3 ] > 0 && p[ 5 ] > 0 )
446  {
447  continue;
448  }
449 
450  val9_1( id, p, i, j );
451 
452  if( Nc8( p ) != 1 )
453  {
454  continue;
455  }
456 
457  val9_2( id, p, i, j );
458 
459  if( max_pixel( p ) <= std::abs( q[ 0 ] ) )
460  {
461  continue;
462  }
463 
464  if( ( q[ 0 ] > q[ 1 ] && q[ 1 ] > 0 ) && ( ( p[ 2 ] > 0 ) || ( p[ 3 ] > 0 ) ) && ( ( q[ 7 ] > 0 ) || ( q[ 8 ] > 0 ) ) )
465  {
466  continue;
467  }
468 
469  if( ( q[ 0 ] > q[ 7 ] && q[ 7 ] > 0 ) && ( ( p[ 5 ] > 0 ) || ( q[ 6 ] >= q[ 0 ] ) ) && ( ( q[ 8 ] > 0 ) || ( q[ 1 ] > 0 ) ) )
470  {
471  continue;
472  }
473 
474  id( i, j ) = - ic( i, j );
475  }
476  }
477 
478  //細め処理　3回目
479  ic = id;
480  for( j = h - 2 ; j > 0 ; j-- )
481  {
482  for( i = 1 ; i < w - 1 ; i++ )
483  {
484  val9_3( ic, p, i, j );
485  val9_3( id, q, i, j );
486 
487  if( q[ 0 ] <= 0 )
488  {
489  continue;
490  }
491 
492  if( p[ 3 ] > 0 && p[ 5 ] > 0 )
493  {
494  continue;
495  }
496 
497  val9_1( ic, p, i, j );
498 
499  if( Nc8( p ) != 1 )
500  {
501  continue;
502  }
503 
504  val9_3( ic, p, i, j );
505 
506  if( max_pixel( p ) <= std::abs( q[ 0 ] ) )
507  {
508  continue;
509  }
510 
511  ic( i, j ) = - id( i, j );
512  }
513  }
514 
515  //後処理
516  id = ic;
517  for( j = 1 ; j < h - 1 ; j++ )
518  {
519  for( i = 1 ; i < w - 1 ; i++ )
520  {
521  val9_1( id, p, i, j );
522  val9_1( ic, q, i, j );
523 
524  if( q[ 0 ] <= 0 )
525  {
526  id( i, j ) = 0;
527  continue;
528  }
529 
530  if( ( p[ 1 ] > 0 ) && ( p[ 3 ] > 0 ) && ( p[ 5 ] > 0 ) && ( p[ 7 ] > 0 ) )
531  {
532  continue;
533  }
534 
535  if( Nc8( p ) != 1 )
536  {
537  continue;
538  }
539 
540  if( plus_pixel_num( p ) < 2 )
541  {
542  continue;
543  }
544 
545  if( plus_pixel_num( p ) == 2 && ( q[ 4 ] > 0 && p[ 4 ] == 0 ) )
546  {
547  p[ 4 ] = 1;
548 
549  if( Nc8( p ) != 1 )
550  {
551  continue;
552  }
553  }
554 
555  id( i, j ) = 0;
556  }
557  }
558 
559  for( i = 0 ; i < ia.size( ) ; i++ )
560  {
561  ia[ i ] = id[ i ] != 0 ? 1 : 0;
562  }
563  }
564 
565  // 消去可能性判定
566  // - 参考文献
567  // - 米倉達広, 横井茂樹, 鳥脇純一郎, 福村晃夫, ``3次元ディジタル画像における1-要素の消去可能性と図形収縮について,'' 電子情報通信学会論文誌D, vol.65-D, no.12, pp.1543-1549, 1982
568  //
569  template < size_t Nc >
570  struct neighbor
571  {
572  typedef size_t size_type;
573  typedef ptrdiff_t difference_type;
574 
575  template < class T >
576  static T P1( T *p[ 4 ], difference_type i, difference_type j )
577  {
578  if( j >= 9 )
579  {
580  return( p[ i ][ j - 8 ] );
581  }
582  else
583  {
584  return( p[ i ][ j ] );
585  }
586  }
587 
588  template < class T >
589  static T P2( T *p[ 4 ], difference_type i, difference_type j )
590  {
591  if( j >= 9 )
592  {
593  return( p[ i ][ j - 8 ] );
594  }
595  else if( j <= 0 )
596  {
597  return( p[ i ][ j + 8 ] );
598  }
599  else
600  {
601  return( p[ i ][ j ] );
602  }
603  }
604 
605  //連結数 6連結
606  template < class T >
607  static difference_type Nc6( T *p[ 4 ] )
608  {
609  static size_type S0[] = { 1, 3 };
610  static size_type S1[] = { 1, 3, 5, 7 };
611 
612  difference_type ret = 0;
613 
614  for( size_type h = 0 ; h < 2 ; h++ )
615  {
616  difference_type tmp = 0;
617  for( size_type k = 0 ; k < 4 ; k++ )
618  {
619  tmp += P1( p, S0[ h ], S1[ k ] ) * P1( p, 2, S1[ k ] );
620  }
621  ret += P1( p, S0[ h ], 0 ) * ( 1 - tmp );
622  }
623 
624  for( size_type k = 0 ; k < 4 ; k++ )
625  {
626  difference_type tmp = 0;
627  for( size_type h = 0 ; h < 2 ; h++ )
628  {
629  tmp += P1( p, S0[ h ], 0 ) * P1( p, S0[ h ], S1[ k ] ) * P1( p, S0[ h ], S1[ k ] + 1 ) * P1( p, S0[ h ], S1[ k ] + 2 );
630  }
631  ret += P1( p, 2, S1[ k ] ) * ( 1 - P1( p, 2, S1[ k ] + 1 ) * P1( p, 2, S1[ k ] + 2 ) * ( 1 - tmp ) );
632  }
633 
634  return( ret );
635  }
636 
637  //連結数 6連結
638  template < class T >
639  static difference_type Sx( T *p[ 4 ] )
640  {
641  static size_type S0[] = { 1, 3 };
642  static size_type S1[] = { 1, 3, 5, 7 };
643 
644  difference_type sx = 0;
645 
646  for( size_type k = 0 ; k < 4 ; k++ )
647  {
648  for( size_type h = 0 ; h < 2 ; h++ )
649  {
650  sx += ( 1 - P1( p, S0[ h ], S1[ k ] + 1 ) ) * P1( p, S0[ h ], S1[ k ] ) * P1( p, S0[ h ], S1[ k ] + 2 ) * P1( p, S0[ h ], 0 ) * P1( p, 2, S1[ k ] ) * P1( p, 2, S1[ k ] + 1 ) * P1( p, 2, S1[ k ] + 2 );
651  }
652  }
653 
654  return( sx );
655  }
656 
657  static bool is_deletable( const int v[ 27 ] )
658  {
659  // 近傍情報を設定
660  const int *p[ 4 ] = { NULL, v + 0, v + 9, v + 18 };
661 
662  //１．連結数
663  if( Nc6( p ) != 1 )
664  {
665  return( false );
666  }
667 
668  //２．6近傍にある1-要素の個数
669  difference_type num = p[ 1 ][ 0 ] + p[ 2 ][ 1 ] + p[ 2 ][ 3 ] + p[ 2 ][ 5 ] + p[ 2 ][ 7 ] + p[ 3 ][ 0 ];
670  if( num <= 3 )
671  {
672  return( true );
673  }
674  //３．
675  else if( num == 4 )
676  {
677  if( Sx( p ) == 1 )
678  {
679  return( false );
680  }
681  else
682  {
683  return( true );
684  }
685  }
686  //４．
687  else if( num == 5 )
688  {
689  if( Sx( p ) == 1 )
690  {
691  return( false );
692  }
693  else
694  {
695  //射影グラフに孤立点があるか
696  if( p[ 2 ][ 0 ] == 0 )
697  {
698  return( false );
699  }
700 
701  static size_type S1[] = { 1, 3, 5, 7 };
702  difference_type hole = 0;
703  if( p[ 1 ][ 0 ] != 0 )
704  {
705  for( size_type h = 0 ; h < 4 ; h++ )
706  {
707  hole += p[ 1 ][ S1[ h ] ] * p[ 2 ][ S1[ h ] ];
708  }
709  if( hole == 0 )
710  {
711  return( false );
712  }
713  }
714 
715  hole = 0;
716  if( p[ 3 ][ 0 ] != 0 )
717  {
718  for( size_type h = 0 ; h < 4 ; h++ )
719  {
720  hole += p[ 3 ][ S1[ h ] ] * p[ 2 ][ S1[ h ] ];
721  }
722  if( hole == 0 )
723  {
724  return( false );
725  }
726  }
727 
728  for( size_type h = 0 ; h < 4 ; h++ )
729  {
730  hole = 0;
731  if( p[ 2 ][ S1[ h ] ] != 0 )
732  {
733  hole = p[ 1 ][ 0 ] * p[ 1 ][ S1[ h ] ] + p[ 3 ][ 0 ] * p[ 3 ][ S1[ h ] ] + P2( p, 2, S1[ h ] - 1 ) * P2( p, 2, S1[ h ] - 2 ) + P2( p, 2, S1[ h ] + 1 ) * P2( p, 2, S1[ h ] + 2 );
734  if( hole == 0 )
735  {
736  return( false );
737  }
738  }
739  }
740  return( true );
741  }
742  }
743  //５．
744  else if( num == 6 )
745  {
746  typedef matrix< difference_type > matrix_type;
747 
748  //隣接行列
749  matrix_type M( 6, 6 );
750 
751  M( 0, 0 ) = 0 ; M( 0, 1 ) = p[ 1 ][ 1 ]; M( 0, 2 ) = p[ 1 ][ 3 ]; M( 0, 3 ) = p[ 1 ][ 7 ]; M( 0, 4 ) = p[ 1 ][ 5 ]; M( 0, 5 ) = 0;
752  M( 1, 0 ) = p[ 1 ][ 1 ]; M( 1, 1 ) = 0 ; M( 1, 2 ) = p[ 2 ][ 2 ]; M( 1, 3 ) = p[ 2 ][ 8 ]; M( 1, 4 ) = 0 ; M( 1, 5 ) = p[ 3 ][ 1 ];
753  M( 2, 0 ) = p[ 1 ][ 3 ]; M( 2, 1 ) = p[ 2 ][ 2 ]; M( 2, 2 ) = 0 ; M( 2, 3 ) = 0 ; M( 2, 4 ) = p[ 2 ][ 4 ]; M( 2, 5 ) = p[ 3 ][ 3 ];
754  M( 3, 0 ) = p[ 1 ][ 7 ]; M( 3, 1 ) = p[ 2 ][ 8 ]; M( 3, 2 ) = 0 ; M( 3, 3 ) = 0 ; M( 3, 4 ) = p[ 2 ][ 6 ]; M( 3, 5 ) = p[ 3 ][ 7 ];
755  M( 4, 0 ) = p[ 1 ][ 5 ]; M( 4, 1 ) = 0 ; M( 4, 2 ) = p[ 2 ][ 4 ]; M( 4, 3 ) = p[ 2 ][ 6 ]; M( 4, 4 ) = 0 ; M( 4, 5 ) = p[ 3 ][ 5 ];
756  M( 5, 0 ) = 0 ; M( 5, 1 ) = p[ 3 ][ 1 ]; M( 5, 2 ) = p[ 3 ][ 3 ]; M( 5, 3 ) = p[ 3 ][ 7 ]; M( 5, 4 ) = p[ 3 ][ 5 ]; M( 5, 5 ) = 0;
757 
758  matrix_type I = matrix_type::identity( 6, 6 );
759  matrix_type N = M * matrix_type( I + M * matrix_type( I + M * matrix_type( I + M * ( I + M ) ) ));
760 
761  for( size_type r = 0 ; r < N.size( ) ; r++ )
762  {
763  if( N[ r ] == 0 )
764  {
765  return( false );
766  }
767  }
768  return( true );
769  }
770  else
771  {
772  std::cout << "Error" << std::endl;
773  }
774 
775  return( true );
776  }
777 
778  //３次元単体のチェック
779  static bool is_3d_simplex( const int v[ 27 ] )
780  {
781  // 近傍情報を設定
782  int p[ 3 ][ 3 ][ 3 ];
783 
784  p[ 0 ][ 0 ][ 0 ] = v[ 4 ];
785  p[ 1 ][ 0 ][ 0 ] = v[ 5 ];
786  p[ 2 ][ 0 ][ 0 ] = v[ 6 ];
787  p[ 0 ][ 1 ][ 0 ] = v[ 3 ];
788  p[ 1 ][ 1 ][ 0 ] = v[ 0 ];
789  p[ 2 ][ 1 ][ 0 ] = v[ 7 ];
790  p[ 0 ][ 2 ][ 0 ] = v[ 2 ];
791  p[ 1 ][ 2 ][ 0 ] = v[ 1 ];
792  p[ 2 ][ 2 ][ 0 ] = v[ 8 ];
793  p[ 0 ][ 0 ][ 1 ] = v[ 13 ];
794  p[ 1 ][ 0 ][ 1 ] = v[ 14 ];
795  p[ 2 ][ 0 ][ 1 ] = v[ 15 ];
796  p[ 0 ][ 1 ][ 1 ] = v[ 12 ];
797  p[ 1 ][ 1 ][ 1 ] = v[ 9 ];
798  p[ 2 ][ 1 ][ 1 ] = v[ 16 ];
799  p[ 0 ][ 2 ][ 1 ] = v[ 11 ];
800  p[ 1 ][ 2 ][ 1 ] = v[ 10 ];
801  p[ 2 ][ 2 ][ 1 ] = v[ 17 ];
802  p[ 0 ][ 0 ][ 2 ] = v[ 22 ];
803  p[ 1 ][ 0 ][ 2 ] = v[ 23 ];
804  p[ 2 ][ 0 ][ 2 ] = v[ 24 ];
805  p[ 0 ][ 1 ][ 2 ] = v[ 21 ];
806  p[ 1 ][ 1 ][ 2 ] = v[ 18 ];
807  p[ 2 ][ 1 ][ 2 ] = v[ 25 ];
808  p[ 0 ][ 2 ][ 2 ] = v[ 20 ];
809  p[ 1 ][ 2 ][ 2 ] = v[ 19 ];
810  p[ 2 ][ 2 ][ 2 ] = v[ 26 ];
811 
812  if( p[ 1 ][ 1 ][ 1 ] == 0 ) return( false );
813 
814  for( size_type k = 1 ; k < 3 ; k++ )
815  {
816  for( size_type j = 1 ; j < 3 ; j++ )
817  {
818  for( size_type i = 1 ; i < 3 ; i++ )
819  {
820  if( p[ i - 1 ][ j ][ k ] * p[ i ][ j - 1 ][ k ] * p[ i ][ j ][ k - 1 ] *
821  p[ i - 1 ][ j - 1 ][ k ] * p[ i ][ j - 1 ][ k - 1 ] * p[ i - 1 ][ j ][ k - 1 ] *
822  p[ i ][ j ][ k ] * p[ i - 1 ][ j - 1 ][ k - 1 ] != 0 )
823  {
824  return( true );
825  }
826  }
827  }
828  }
829 
830  return( false );
831  }
832  };
833 
834  template < >
835  struct neighbor< 26 >
836  {
837  typedef size_t size_type;
838  typedef ptrdiff_t difference_type;
839 
840  template < class T >
841  static T P1( T *p[ 4 ], difference_type i, difference_type j )
842  {
843  if( j >= 9 )
844  {
845  return( p[ i ][ j - 8 ] );
846  }
847  else
848  {
849  return( p[ i ][ j ] );
850  }
851  }
852 
853  template < class T >
854  static T P2( T *p[ 4 ], difference_type i, difference_type j )
855  {
856  if( j >= 9 )
857  {
858  return( p[ i ][ j - 8 ] );
859  }
860  else if( j <= 0 )
861  {
862  return( p[ i ][ j + 8 ] );
863  }
864  else
865  {
866  return( p[ i ][ j ] );
867  }
868  }
869 
870  template < class T >
871  static T P1_( T *p[ 4 ], difference_type i, difference_type j )
872  {
873  if( j >= 9 )
874  {
875  return( 1 - p[ i ][ j - 8 ] );
876  }
877  else
878  {
879  return( 1 - p[ i ][ j ] );
880  }
881  }
882 
883  template < class T >
884  static T P2_( T *p[ 4 ], difference_type i, difference_type j )
885  {
886  if( j >= 9 )
887  {
888  return( 1 - p[ i ][ j - 8 ] );
889  }
890  else if( j <= 0 )
891  {
892  return( 1 - p[ i ][ j + 8 ] );
893  }
894  else
895  {
896  return( 1 - p[ i ][ j ] );
897  }
898  }
899 
900 
901  // 連結数 26連結
902  template < class T >
903  static difference_type Nc26( T *p[ 4 ] )
904  {
905  static size_type S0[] = { 1, 3 };
906  static size_type S1[] = { 1, 3, 5, 7 };
907 
908  difference_type ret = 0;
909 
910  for( size_type h = 0 ; h < 2 ; h++ )
911  {
912  difference_type tmp = 0;
913  for( size_type k = 0 ; k < 4 ; k++ )
914  {
915  tmp += P1_( p, S0[ h ], S1[ k ] ) * P1_( p, 2, S1[ k ] );
916  }
917  ret += P1_( p, S0[ h ], 0 ) * ( 1 - tmp );
918  }
919 
920  for( size_type k = 0 ; k < 4 ; k++ )
921  {
922  difference_type tmp = 0;
923  for( size_type h = 0 ; h < 2 ; h++ )
924  {
925  tmp += P1_( p, S0[ h ], 0 ) * P1_( p, S0[ h ], S1[ k ] ) * P1_( p, S0[ h ], S1[ k ] + 1 ) * P1_( p, S0[ h ], S1[ k ] + 2 );
926  }
927  ret += P1_( p, 2, S1[ k ] ) * ( 1 - P1_( p, 2, S1[ k ] + 1 ) * P1_( p, 2, S1[ k ] + 2 ) * ( 1 - tmp ) );
928  }
929 
930  ret = 2 - ret;
931 
932  return( ret );
933  }
934 
935  //
936  template < class T >
937  static difference_type Sx( T *p[ 4 ] )
938  {
939  static size_type S0[] = { 1, 3 };
940  static size_type S1[] = { 1, 3, 5, 7 };
941 
942  difference_type sx = 0;
943 
944  for( size_type k = 0 ; k < 4 ; k++ )
945  {
946  for( size_type h = 0 ; h < 2 ; h++ )
947  {
948  sx += ( 1 - P1( p, S0[ h ], S1[ k ] + 1 ) ) * P1( p, S0[ h ], S1[ k ] ) * P1( p, S0[ h ], S1[ k ] + 2 ) * P1( p, S0[ h ], 0 ) * P1( p, 2, S1[ k ] ) * P1( p, 2, S1[ k ] + 1 ) * P1( p, 2, S1[ k ] + 2 );
949  }
950  }
951 
952  return( sx );
953  }
954 
955 
956  static bool is_deletable( const int v[ 27 ] )
957  {
958  // 近傍情報を設定
959  int pp[ 27 ];
960  const int *p[ 4 ] = { NULL, v + 0, v + 9, v + 18 };
961  memcpy( pp, v, sizeof( int ) * 27 );
962 
963  //１．連結数
964  if( Nc26( p ) != 1 )
965  {
966  return( false );
967  }
968 
969  // 近傍情報を取得 0-1反転
970  for( size_t i = 0 ; i < 9 ; i++ )
971  {
972  pp[ i ] = v[ i ] == 0 ? 1 : 0;
973  }
974  pp[ 9 ] = v[ 9 ];
975  for( size_t i = 10 ; i < 27 ; i++ )
976  {
977  pp[ i ] = v[ i ] == 0 ? 1 : 0;
978  }
979 
980  p[ 1 ] = pp;
981  p[ 2 ] = pp + 9;
982  p[ 3 ] = pp + 18;
983 
984  //２．6近傍にある1-要素の個数
985  difference_type num = p[ 1 ][ 0 ] + p[ 2 ][ 1 ] + p[ 2 ][ 3 ] + p[ 2 ][ 5 ] + p[ 2 ][ 7 ] + p[ 3 ][ 0 ];
986  if( num <= 3 )
987  {
988  return( true );
989  }
990  //３．
991  else if( num == 4 )
992  {
993  if( Sx( p ) == 1 )
994  {
995  return( false );
996  }
997  else
998  {
999  return( true );
1000  }
1001  }
1002  //４．
1003  else if( num == 5 )
1004  {
1005  if( Sx( p ) == 1 )
1006  {
1007  return( false );
1008  }
1009  else
1010  {
1011  //射影グラフに孤立点があるか
1012  if( p[ 2 ][ 0 ] == 0 )
1013  {
1014  return( false );
1015  }
1016 
1017  static size_type S1[] = { 1, 3, 5, 7 };
1018  difference_type hole = 0;
1019  if( p[ 1 ][ 0 ] != 0 )
1020  {
1021  for( size_type h = 0 ; h < 4 ; h++ )
1022  {
1023  hole += p[ 1 ][ S1[ h ] ] * p[ 2 ][ S1[ h ] ];
1024  }
1025  if( hole == 0 )
1026  {
1027  return( false );
1028  }
1029  }
1030 
1031  hole = 0;
1032  if( p[ 3 ][ 0 ] != 0 )
1033  {
1034  for( size_type h = 0 ; h < 4 ; h++ )
1035  {
1036  hole += p[ 3 ][ S1[ h ] ] * p[ 2 ][ S1[ h ] ];
1037  }
1038  if( hole == 0 )
1039  {
1040  return( false );
1041  }
1042  }
1043 
1044  for( size_type h = 0 ; h < 4 ; h++ )
1045  {
1046  hole = 0;
1047  if( p[ 2 ][ S1[ h ] ] != 0 )
1048  {
1049  hole = p[ 1 ][ 0 ] * p[ 1 ][ S1[ h ] ] + p[ 3 ][ 0 ] * p[ 3 ][ S1[ h ] ] + P2( p, 2, S1[ h ] - 1 ) * P2( p, 2, S1[ h ] - 2 ) + P2( p, 2, S1[ h ] + 1 ) * P2( p, 2, S1[ h ] + 2 );
1050  if( hole == 0 )
1051  {
1052  return( false );
1053  }
1054  }
1055  }
1056  return( true );
1057  }
1058  }
1059  //５．
1060  else if( num == 6 )
1061  {
1062  typedef matrix< int > matrix_type;
1063 
1064  //隣接行列
1065  matrix_type M( 6, 6 );
1066 
1067  M( 0, 0 ) = 0 ; M( 0, 1 ) = p[ 1 ][ 1 ]; M( 0, 2 ) = p[ 1 ][ 3 ]; M( 0, 3 ) = p[ 1 ][ 7 ]; M( 0, 4 ) = p[ 1 ][ 5 ]; M( 0, 5 ) = 0;
1068  M( 1, 0 ) = p[ 1 ][ 1 ]; M( 1, 1 ) = 0 ; M( 1, 2 ) = p[ 2 ][ 2 ]; M( 1, 3 ) = p[ 2 ][ 8 ]; M( 1, 4 ) = 0 ; M( 1, 5 ) = p[ 3 ][ 1 ];
1069  M( 2, 0 ) = p[ 1 ][ 3 ]; M( 2, 1 ) = p[ 2 ][ 2 ]; M( 2, 2 ) = 0 ; M( 2, 3 ) = 0 ; M( 2, 4 ) = p[ 2 ][ 4 ]; M( 2, 5 ) = p[ 3 ][ 3 ];
1070  M( 3, 0 ) = p[ 1 ][ 7 ]; M( 3, 1 ) = p[ 2 ][ 8 ]; M( 3, 2 ) = 0 ; M( 3, 3 ) = 0 ; M( 3, 4 ) = p[ 2 ][ 6 ]; M( 3, 5 ) = p[ 3 ][ 7 ];
1071  M( 4, 0 ) = p[ 1 ][ 5 ]; M( 4, 1 ) = 0 ; M( 4, 2 ) = p[ 2 ][ 4 ]; M( 4, 3 ) = p[ 2 ][ 6 ]; M( 4, 4 ) = 0 ; M( 4, 5 ) = p[ 3 ][ 5 ];
1072  M( 5, 0 ) = 0 ; M( 5, 1 ) = p[ 3 ][ 1 ]; M( 5, 2 ) = p[ 3 ][ 3 ]; M( 5, 3 ) = p[ 3 ][ 7 ]; M( 5, 4 ) = p[ 3 ][ 5 ]; M( 5, 5 ) = 0;
1073 
1074  matrix_type I = matrix_type::identity( 6, 6 );
1075  matrix_type N = M * matrix_type( I + M * matrix_type( I + M * matrix_type( I + M * ( I + M ) ) ));
1076 
1077  for( size_type r = 0 ; r < N.size( ) ; r++ )
1078  {
1079  if( N[ r ] == 0 )
1080  {
1081  return( false );
1082  }
1083  }
1084  return( true );
1085  }
1086  else
1087  {
1088  std::cout << "Error" << std::endl;
1089  }
1090 
1091  return( true );
1092  }
1093 
1094  //３次元単体のチェック
1095  static bool is_3d_simplex( const int v[ 27 ] )
1096  {
1097  // 近傍情報を設定
1098  int p[ 3 ][ 3 ][ 3 ];
1099 
1100  p[ 0 ][ 0 ][ 0 ] = v[ 4 ];
1101  p[ 1 ][ 0 ][ 0 ] = v[ 5 ];
1102  p[ 2 ][ 0 ][ 0 ] = v[ 6 ];
1103  p[ 0 ][ 1 ][ 0 ] = v[ 3 ];
1104  p[ 1 ][ 1 ][ 0 ] = v[ 0 ];
1105  p[ 2 ][ 1 ][ 0 ] = v[ 7 ];
1106  p[ 0 ][ 2 ][ 0 ] = v[ 2 ];
1107  p[ 1 ][ 2 ][ 0 ] = v[ 1 ];
1108  p[ 2 ][ 2 ][ 0 ] = v[ 8 ];
1109  p[ 0 ][ 0 ][ 1 ] = v[ 13 ];
1110  p[ 1 ][ 0 ][ 1 ] = v[ 14 ];
1111  p[ 2 ][ 0 ][ 1 ] = v[ 15 ];
1112  p[ 0 ][ 1 ][ 1 ] = v[ 12 ];
1113  p[ 1 ][ 1 ][ 1 ] = v[ 9 ];
1114  p[ 2 ][ 1 ][ 1 ] = v[ 16 ];
1115  p[ 0 ][ 2 ][ 1 ] = v[ 11 ];
1116  p[ 1 ][ 2 ][ 1 ] = v[ 10 ];
1117  p[ 2 ][ 2 ][ 1 ] = v[ 17 ];
1118  p[ 0 ][ 0 ][ 2 ] = v[ 22 ];
1119  p[ 1 ][ 0 ][ 2 ] = v[ 23 ];
1120  p[ 2 ][ 0 ][ 2 ] = v[ 24 ];
1121  p[ 0 ][ 1 ][ 2 ] = v[ 21 ];
1122  p[ 1 ][ 1 ][ 2 ] = v[ 18 ];
1123  p[ 2 ][ 1 ][ 2 ] = v[ 25 ];
1124  p[ 0 ][ 2 ][ 2 ] = v[ 20 ];
1125  p[ 1 ][ 2 ][ 2 ] = v[ 19 ];
1126  p[ 2 ][ 2 ][ 2 ] = v[ 26 ];
1127 
1128  if( p[ 1 ][ 1 ][ 1 ] == 0 ) return( false );
1129 
1130  for( size_type k = 0 ; k < 3 ; k += 2 )
1131  {
1132  for( size_type j = 0 ; j < 3 ; j += 2 )
1133  {
1134  for( size_type i = 0 ; i < 3 ; i += 2 )
1135  {
1136  int num = 0;
1137 
1138  if( p[ i ][ j ][ 1 ] != 0 ) num++;
1139  if( p[ i ][ j ][ k ] != 0 ) num++;
1140  if( p[ i ][ 1 ][ 1 ] != 0 ) num++;
1141  if( p[ i ][ 1 ][ k ] != 0 ) num++;
1142  if( p[ 1 ][ j ][ 1 ] != 0 ) num++;
1143  if( p[ 1 ][ j ][ k ] != 0 ) num++;
1144  if( p[ 1 ][ 1 ][ k ] != 0 ) num++;
1145 
1146  if( num > 3 ) return( true );
1147  if( num < 3 ) continue;
1148 
1149  if( p[ 1 ][ 1 ][ k ] * ( p[ 1 ][ 1 ][ k ] * p[ i ][ 1 ][ k ] +
1150  p[ 1 ][ j ][ 1 ] * p[ 1 ][ j ][ k ] +
1151  p[ i ][ j ][ 1 ] * p[ i ][ j ][ k ] ) != 0 )
1152  {
1153  continue;
1154  }
1155  if( p[ i ][ 1 ][ 1 ] * ( p[ 1 ][ j ][ 1 ] * p[ i ][ j ][ 1 ] +
1156  p[ 1 ][ j ][ k ] * p[ i ][ j ][ k ] ) != 0 )
1157  {
1158  continue;
1159  }
1160  if( p[ 1 ][ j ][ 1 ] * ( p[ i ][ 1 ][ k ] * p[ i ][ j ][ k ] ) != 0 )
1161  {
1162  continue;
1163  }
1164 
1165  return( true );
1166 
1167  }
1168  }
1169  }
1170  return( false );
1171  }
1172  };
1173 
1174 
1175  template < class Tin, class Tout, class Difference >
1176  void create_neighbor_list( const Tin *pin, Tout *pout, Difference p[ 27 ] )
1177  {
1178  for( size_t i = 0 ; i < 27 ; i++ )
1179  {
1180  pout[ i ] = pin[ p[ i ] ] > 0 ? 1 : 0;
1181  }
1182  }
1183 
1184  template < class T, class Allocator, class Neighbor >
1185  void shrink_skelton( array3< T, Allocator > &in, Neighbor __dmy__ )
1186  {
1187  typedef typename array3< T, Allocator >::size_type size_type;
1188  typedef typename array3< T, Allocator >::difference_type difference_type;
1189  typedef typename array3< T, Allocator >::const_pointer const_pointer;
1190  typedef typename array3< T, Allocator >::value_type value_type;
1191  typedef Neighbor neighbor_type;
1192 
1193  difference_type diff[ 27 ];
1194  int val[ 27 ];
1195 
1196  {
1197  difference_type ox = in.width( ) / 2;
1198  difference_type oy = in.height( ) / 2;
1199  difference_type oz = in.depth( ) / 2;
1200 
1201  const_pointer p0 = &in( ox, oy, oz );
1202 
1203  diff[ 0 ] = &in( ox , oy , oz - 1 ) - p0;
1204  diff[ 1 ] = &in( ox , oy + 1, oz - 1 ) - p0;
1205  diff[ 2 ] = &in( ox - 1, oy + 1, oz - 1 ) - p0;
1206  diff[ 3 ] = &in( ox - 1, oy , oz - 1 ) - p0;
1207  diff[ 4 ] = &in( ox - 1, oy - 1, oz - 1 ) - p0;
1208  diff[ 5 ] = &in( ox , oy - 1, oz - 1 ) - p0;
1209  diff[ 6 ] = &in( ox + 1, oy - 1, oz - 1 ) - p0;
1210  diff[ 7 ] = &in( ox + 1, oy , oz - 1 ) - p0;
1211  diff[ 8 ] = &in( ox + 1, oy + 1, oz - 1 ) - p0;
1212 
1213  diff[ 9 ] = &in( ox , oy , oz ) - p0;
1214  diff[ 10 ] = &in( ox , oy + 1, oz ) - p0;
1215  diff[ 11 ] = &in( ox - 1, oy + 1, oz ) - p0;
1216  diff[ 12 ] = &in( ox - 1, oy , oz ) - p0;
1217  diff[ 13 ] = &in( ox - 1, oy - 1, oz ) - p0;
1218  diff[ 14 ] = &in( ox , oy - 1, oz ) - p0;
1219  diff[ 15 ] = &in( ox + 1, oy - 1, oz ) - p0;
1220  diff[ 16 ] = &in( ox + 1, oy , oz ) - p0;
1221  diff[ 17 ] = &in( ox + 1, oy + 1, oz ) - p0;
1222 
1223  diff[ 18 ] = &in( ox , oy , oz + 1 ) - p0;
1224  diff[ 19 ] = &in( ox , oy + 1, oz + 1 ) - p0;
1225  diff[ 20 ] = &in( ox - 1, oy + 1, oz + 1 ) - p0;
1226  diff[ 21 ] = &in( ox - 1, oy , oz + 1 ) - p0;
1227  diff[ 22 ] = &in( ox - 1, oy - 1, oz + 1 ) - p0;
1228  diff[ 23 ] = &in( ox , oy - 1, oz + 1 ) - p0;
1229  diff[ 24 ] = &in( ox + 1, oy - 1, oz + 1 ) - p0;
1230  diff[ 25 ] = &in( ox + 1, oy , oz + 1 ) - p0;
1231  diff[ 26 ] = &in( ox + 1, oy + 1, oz + 1 ) - p0;
1232  }
1233 
1234  // 図形の端は0
1235  for( size_type j = 0 ; j < in.height( ) ; j++ )
1236  {
1237  for( size_type i = 0 ; i < in.width( ) ; i++ )
1238  {
1239  in( i, j, 0 ) = 0;
1240  in( i, j, in.depth( ) - 1 ) = 0;
1241  }
1242  }
1243 
1244  for( size_type k = 0 ; k < in.depth( ) ; k++ )
1245  {
1246  for( size_type j = 0 ; j < in.height( ) ; j++ )
1247  {
1248  in( 0, j, k ) = 0;
1249  in( in.width( ) - 1, j, k ) = 0;
1250  }
1251  }
1252 
1253  for( size_type k = 0 ; k < in.depth( ) ; k++ )
1254  {
1255  for( size_type i = 0 ; i < in.width( ) ; i++ )
1256  {
1257  in( i, 0, k ) = 0;
1258  in( i, in.height( ) - 1, k ) = 0;
1259  }
1260  }
1261 
1262  bool loop;
1263  do
1264  {
1265  loop = false;
1266  for( size_type k = 1 ; k < in.depth( ) - 1 ; k++ )
1267  {
1268  for( size_type j = 1 ; j < in.height( ) - 1 ; j++ )
1269  {
1270  for( size_type i = 1 ; i < in.width( ) - 1 ; i++ )
1271  {
1272  value_type &v = in( i, j, k );
1273  if( v != 0 )
1274  {
1275  create_neighbor_list( &v, val, diff );
1276 
1277  if( neighbor_type::is_deletable( val ) )
1278  {
1279  v = 0;
1280  loop = true;
1281  }
1282  }
1283  }
1284  }
1285  }
1286  } while( loop );
1287  }
1288 
1289  //　ユークリッド距離変換を用いた3次元ディジタル画像の細線化と薄面化
1290  // - 参考文献
1291  // - 斉藤豊文, 森健策, 鳥脇純一郎, ``ユークリッド距離変換を用いた3次元ディジタル画像の薄面化および細線化の逐次型アルゴリズムとその諸性質,'' 電子情報通信学会論文誌D-II, vol.J79-D-II, no.10, pp.1675-1685, 1996
1292  //
1293  template< class T >
1294  struct border
1295  {
1296  typedef ptrdiff_t difference_type;
1297 
1298  difference_type diff;
1299  T value;
1300  T distance;
1301 
1302  border( difference_type d, T val ) : diff( d ), value( val ), distance( val ) {}
1303  };
1304 
1305  template < class T >
1306  inline std::ostream &operator <<( std::ostream &out, const border< T > &b )
1307  {
1308  out << "( ";
1309  out << b.i << ", ";
1310  out << b.j << ", ";
1311  out << b.k;
1312  out << " ) = ( ";
1313  out << b.value << ", ";
1314  out << b.distance;
1315  out << " )";
1316  return( out );
1317  }
1318 
1319  template < class T, class Allocator, class Neighbor >
1320  void thinning( array3< T, Allocator > &in, Neighbor __dmy__ )
1321  {
1322  typedef typename array3< T, Allocator >::size_type size_type;
1323  typedef typename array3< T, Allocator >::difference_type difference_type;
1324  typedef typename array3< T, Allocator >::const_pointer const_pointer;
1325  typedef typename array3< T, Allocator >::pointer pointer;
1326  typedef T value_type;
1327  typedef Neighbor neighbor_type;
1328  typedef border< T > border_type;
1329  typedef std::vector< border_type > border_list_type;
1330 
1331  // Step1 距離変換
1332  euclidean::distance_transform( in, in );
1333 
1334  // 図形の端は0
1335  for( size_type j = 0 ; j < in.height( ) ; j++ )
1336  {
1337  for( size_type i = 0 ; i < in.width( ) ; i++ )
1338  {
1339  in( i, j, 0 ) = 0;
1340  in( i, j, in.depth( ) - 1 ) = 0;
1341  }
1342  }
1343 
1344  for( size_type k = 0 ; k < in.depth( ) ; k++ )
1345  {
1346  for( size_type j = 0 ; j < in.height( ) ; j++ )
1347  {
1348  in( 0, j, k ) = 0;
1349  in( in.width( ) - 1, j, k ) = 0;
1350  }
1351  }
1352 
1353  for( size_type k = 0 ; k < in.depth( ) ; k++ )
1354  {
1355  for( size_type i = 0 ; i < in.width( ) ; i++ )
1356  {
1357  in( i, 0, k ) = 0;
1358  in( i, in.height( ) - 1, k ) = 0;
1359  }
1360  }
1361 
1362  value_type min = type_limits< value_type >::maximum( ), max = type_limits< value_type >::minimum( );
1363  for( size_type i = 0 ; i < in.size( ) ; i++ )
1364  {
1365  value_type &v = in[ i ];
1366  if( v != 0 )
1367  {
1368  v += 20;
1369  if( v > max ) max = v;
1370  if( v < min ) min = v;
1371  }
1372  }
1373 
1374  border_list_type blist, slist[ 9 ];
1375  difference_type diff[ 27 ], nc6[ 6 ];
1376  int val[ 27 ];
1377 
1378  {
1379  difference_type ox = in.width( ) / 2;
1380  difference_type oy = in.height( ) / 2;
1381  difference_type oz = in.depth( ) / 2;
1382 
1383  const_pointer p0 = &in( ox, oy, oz );
1384 
1385 
1386  diff[ 0 ] = &in( ox , oy , oz - 1 ) - p0;
1387  diff[ 1 ] = &in( ox , oy + 1, oz - 1 ) - p0;
1388  diff[ 2 ] = &in( ox - 1, oy + 1, oz - 1 ) - p0;
1389  diff[ 3 ] = &in( ox - 1, oy , oz - 1 ) - p0;
1390  diff[ 4 ] = &in( ox - 1, oy - 1, oz - 1 ) - p0;
1391  diff[ 5 ] = &in( ox , oy - 1, oz - 1 ) - p0;
1392  diff[ 6 ] = &in( ox + 1, oy - 1, oz - 1 ) - p0;
1393  diff[ 7 ] = &in( ox + 1, oy , oz - 1 ) - p0;
1394  diff[ 8 ] = &in( ox + 1, oy + 1, oz - 1 ) - p0;
1395 
1396  diff[ 9 ] = &in( ox , oy , oz ) - p0;
1397  diff[ 10 ] = &in( ox , oy + 1, oz ) - p0;
1398  diff[ 11 ] = &in( ox - 1, oy + 1, oz ) - p0;
1399  diff[ 12 ] = &in( ox - 1, oy , oz ) - p0;
1400  diff[ 13 ] = &in( ox - 1, oy - 1, oz ) - p0;
1401  diff[ 14 ] = &in( ox , oy - 1, oz ) - p0;
1402  diff[ 15 ] = &in( ox + 1, oy - 1, oz ) - p0;
1403  diff[ 16 ] = &in( ox + 1, oy , oz ) - p0;
1404  diff[ 17 ] = &in( ox + 1, oy + 1, oz ) - p0;
1405 
1406  diff[ 18 ] = &in( ox , oy , oz + 1 ) - p0;
1407  diff[ 19 ] = &in( ox , oy + 1, oz + 1 ) - p0;
1408  diff[ 20 ] = &in( ox - 1, oy + 1, oz + 1 ) - p0;
1409  diff[ 21 ] = &in( ox - 1, oy , oz + 1 ) - p0;
1410  diff[ 22 ] = &in( ox - 1, oy - 1, oz + 1 ) - p0;
1411  diff[ 23 ] = &in( ox , oy - 1, oz + 1 ) - p0;
1412  diff[ 24 ] = &in( ox + 1, oy - 1, oz + 1 ) - p0;
1413  diff[ 25 ] = &in( ox + 1, oy , oz + 1 ) - p0;
1414  diff[ 26 ] = &in( ox + 1, oy + 1, oz + 1 ) - p0;
1415 
1416  nc6[ 0 ] = diff[ 12 ];
1417  nc6[ 1 ] = diff[ 16 ];
1418  nc6[ 2 ] = diff[ 14 ];
1419  nc6[ 3 ] = diff[ 10 ];
1420  nc6[ 4 ] = diff[ 0 ];
1421  nc6[ 5 ] = diff[ 18 ];
1422  }
1423 
1424  {
1425  const_pointer p0 = &in[ 0 ];
1426 
1427  // Step2 境界画素の検出
1428  for( size_type i = 0 ; i < in.size( ) ; i++ )
1429  {
1430  value_type &v = in[ i ];
1431  if( v > 20 )
1432  {
1433  const_pointer p = &v;
1434  if( p[ nc6[ 0 ] ] == 0 || p[ nc6[ 1 ] ] == 0 ||
1435  p[ nc6[ 2 ] ] == 0 || p[ nc6[ 3 ] ] == 0 ||
1436  p[ nc6[ 4 ] ] == 0 || p[ nc6[ 5 ] ] == 0 )
1437  {
1438  blist.push_back( border_type( p - p0, v ) );
1439  v = 1;
1440  }
1441  }
1442  }
1443  }
1444 
1445  // 以降の処理のためにメモリの予備を確保する
1446  blist.reserve( blist.size( ) * 2 + 1 );
1447  slist[ 0 ].reserve( blist.size( ) );
1448  slist[ 1 ].reserve( blist.size( ) );
1449  slist[ 2 ].reserve( blist.size( ) );
1450  slist[ 3 ].reserve( blist.size( ) );
1451  slist[ 4 ].reserve( blist.size( ) );
1452  slist[ 5 ].reserve( blist.size( ) );
1453  slist[ 6 ].reserve( blist.size( ) );
1454  slist[ 7 ].reserve( blist.size( ) );
1455  slist[ 8 ].reserve( blist.size( ) );
1456 
1457  size_type count = 0, loop = 0;
1458 
1459  do
1460  {
1461  size_type __length__ = 0;
1462 
1463  //Step3 メインサイクル
1464  for( size_type l = 0 ; l < blist.size( ) ; l++ )
1465  {
1466  border_type &b = blist[ l ];
1467  pointer p = &in[ b.diff ];
1468 
1469  if( b.value <= min )
1470  {
1471  create_neighbor_list( p, val, diff );
1472 
1473  // 消去不可能なら一時保存点
1474  if( !neighbor_type::is_deletable( val ) )
1475  {
1476  b.value = 16;
1477  }
1478  else
1479  {
1480  // 自分自身を除く
1481  size_type num = val[ 0 ];
1482  num += val[ 1 ];
1483  num += val[ 2 ];
1484  num += val[ 3 ];
1485  num += val[ 4 ];
1486  num += val[ 5 ];
1487  num += val[ 6 ];
1488  num += val[ 7 ];
1489  num += val[ 8 ];
1490  num += val[ 10 ];
1491  num += val[ 11 ];
1492  num += val[ 12 ];
1493  num += val[ 13 ];
1494  num += val[ 14 ];
1495  num += val[ 15 ];
1496  num += val[ 16 ];
1497  num += val[ 17 ];
1498  num += val[ 18 ];
1499  num += val[ 19 ];
1500  num += val[ 20 ];
1501  num += val[ 21 ];
1502  num += val[ 22 ];
1503  num += val[ 23 ];
1504  num += val[ 24 ];
1505  num += val[ 25 ];
1506  num += val[ 26 ];
1507 
1508  // 端点なら永久保存点
1509  if( num == 1 )
1510  {
1511  continue;
1512  }
1513  else
1514  {
1515  difference_type val = num / 3 + 7;
1516  b.value = static_cast< value_type >( val );
1517 
1518  if( 7 <= val && val <= 15 )
1519  {
1520  slist[ val - 7 ].push_back( b );
1521  continue;
1522  }
1523  }
1524  }
1525  }
1526 
1527  blist[ __length__++ ] = b;
1528  }
1529 
1530  blist.erase( blist.begin( ) + __length__, blist.end( ) );
1531 
1532  // Step4 サブサイクル
1533  for( size_type ll = 0 ; ll < 9 ; ll++ )
1534  {
1535  border_list_type &list = slist[ ll ];
1536 
1537  for( size_type l = 0 ; l < list.size( ) ; l++ )
1538  {
1539  border_type &b = list[ l ];
1540  pointer p = &in[ b.diff ];
1541 
1542  create_neighbor_list( p, val, diff );
1543 
1544  //消去不可能なら一時保存点
1545  if( !neighbor_type::is_deletable( val ) )
1546  {
1547  b.value = 16;
1548  blist.push_back( b );
1549  }
1550  else
1551  {
1552  size_type num = 0, i;
1553 
1554  // 自分自身を除く周りを調べる
1555  for( i = 0 ; i < 9 && num < 2 ; i++ )
1556  {
1557  if( val[ i ] > 0 )
1558  {
1559  num++;
1560  }
1561  }
1562  for( i = 10 ; i < 27 && num < 2 ; i++ )
1563  {
1564  if( val[ i ] > 0 )
1565  {
1566  num++;
1567  }
1568  }
1569 
1570  //端点なら永久保存点
1571  if( num > 1 )
1572  {
1573  //画素の消去
1574  p[ 0 ] = 0;
1575 
1576  if( p[ nc6[ 0 ] ] > 20 )
1577  {
1578  blist.push_back( border_type( b.diff + nc6[ 0 ], p[ nc6[ 0 ] ] ) );
1579  p[ nc6[ 0 ] ] = 1;
1580  }
1581  if( p[ nc6[ 1 ] ] > 20 )
1582  {
1583  blist.push_back( border_type( b.diff + nc6[ 1 ], p[ nc6[ 1 ] ] ) );
1584  p[ nc6[ 1 ] ] = 1;
1585  }
1586  if( p[ nc6[ 2 ] ] > 20 )
1587  {
1588  blist.push_back( border_type( b.diff + nc6[ 2 ], p[ nc6[ 2 ] ] ) );
1589  p[ nc6[ 2 ] ] = 1;
1590  }
1591  if( p[ nc6[ 3 ] ] > 20 )
1592  {
1593  blist.push_back( border_type( b.diff + nc6[ 3 ], p[ nc6[ 3 ] ] ) );
1594  p[ nc6[ 3 ] ] = 1;
1595  }
1596  if( p[ nc6[ 4 ] ] > 20 )
1597  {
1598  blist.push_back( border_type( b.diff + nc6[ 4 ], p[ nc6[ 4 ] ] ) );
1599  p[ nc6[ 4 ] ] = 1;
1600  }
1601  if( p[ nc6[ 5 ] ] > 20 )
1602  {
1603  blist.push_back( border_type( b.diff + nc6[ 5 ], p[ nc6[ 5 ] ] ) );
1604  p[ nc6[ 5 ] ] = 1;
1605  }
1606  }
1607  }
1608  }
1609 
1610  // 次以降の処理のためにリスト空にする（内部で使用するメモリ容量は変化しない）
1611  list.clear( );
1612  }
1613 
1614  for( size_type l = 0 ; l < blist.size( ) ; l++ )
1615  {
1616  border_type &b = blist[ l ];
1617  pointer p = &in[ b.diff ];
1618 
1619  if( b.value == 16 )
1620  {
1621  create_neighbor_list( p, val, diff );
1622 
1623  if( neighbor_type::is_deletable( val ) )
1624  {
1625  b.value = b.distance;
1626  }
1627  }
1628  }
1629 
1630  //Step5 終了判定
1631  min = max;
1632  count = 0;
1633  for( size_type l = 0 ; l < blist.size( ) ; l++ )
1634  {
1635  border_type &b = blist[ l ];
1636  if( b.value < min && b.value > 20 )
1637  {
1638  min = b.value;
1639  }
1640 
1641  if( b.value == 16 )
1642  {
1643  count++;
1644  }
1645  }
1646 
1647  std::cout << loop++ << " \r";
1648 
1649  } while( min < max || count != blist.size( ) );
1650 
1651  std::cout << std::endl;
1652 
1653  // Step6 後処理
1654  for( size_type i = 0 ; i < in.size( ) ; i++ )
1655  {
1656  if( in[ i ] != 0 )
1657  {
1658  in[ i ] = 1;
1659  }
1660  }
1661  }
1662 
1663  //薄面化
1664  template < class T, class Allocator, class Neighbor >
1665  void surface_thinning( array3< T, Allocator > &in, Neighbor __dmy__ )
1666  {
1667  typedef typename array3< T, Allocator >::size_type size_type;
1668  typedef typename array3< T, Allocator >::difference_type difference_type;
1669  typedef typename array3< T, Allocator >::const_pointer const_pointer;
1670  typedef typename array3< T, Allocator >::pointer pointer;
1671  typedef T value_type;
1672  typedef Neighbor neighbor_type;
1673  typedef border< T > border_type;
1674  typedef std::vector< border_type > border_list_type;
1675 
1676  // 図形の端は０
1677  for( size_type j = 0 ; j < in.height( ) ; j++ )
1678  {
1679  for( size_type i = 0 ; i < in.width( ) ; i++ )
1680  {
1681  in( i, j, 0 ) = 0;
1682  in( i, j, in.depth( ) - 1 ) = 0;
1683  }
1684  }
1685 
1686  for( size_type k = 0 ; k < in.depth( ) ; k++ )
1687  {
1688  for( size_type j = 0 ; j < in.height( ) ; j++ )
1689  {
1690  in( 0, j, k ) = 0;
1691  in( in.width( ) - 1, j, k ) = 0;
1692  }
1693  }
1694 
1695  for( size_type k = 0 ; k < in.depth( ) ; k++ )
1696  {
1697  for( size_type i = 0 ; i < in.width( ) ; i++ )
1698  {
1699  in( i, 0, k ) = 0;
1700  in( i, in.height( ) - 1, k ) = 0;
1701  }
1702  }
1703 
1704  // Step1 距離変換
1705  euclidean::distance_transform( in, in );
1706 
1707  value_type min = type_limits< value_type >::maximum( ), max = type_limits< value_type >::minimum( );
1708  for( size_type i = 0 ; i < in.size( ) ; i++ )
1709  {
1710  value_type &v = in[ i ];
1711  if( v != 0 )
1712  {
1713  v += 20;
1714  if( v > max ) max = v;
1715  if( v < min ) min = v;
1716  }
1717  }
1718 
1719  border_list_type blist, slist[ 9 ];
1720  difference_type diff[ 27 ], nc6[ 6 ];
1721  int val[ 27 ];
1722 
1723  {
1724  difference_type ox = in.width( ) / 2;
1725  difference_type oy = in.height( ) / 2;
1726  difference_type oz = in.depth( ) / 2;
1727 
1728  const_pointer p0 = &in( ox, oy, oz );
1729 
1730 
1731  diff[ 0 ] = &in( ox , oy , oz - 1 ) - p0;
1732  diff[ 1 ] = &in( ox , oy + 1, oz - 1 ) - p0;
1733  diff[ 2 ] = &in( ox - 1, oy + 1, oz - 1 ) - p0;
1734  diff[ 3 ] = &in( ox - 1, oy , oz - 1 ) - p0;
1735  diff[ 4 ] = &in( ox - 1, oy - 1, oz - 1 ) - p0;
1736  diff[ 5 ] = &in( ox , oy - 1, oz - 1 ) - p0;
1737  diff[ 6 ] = &in( ox + 1, oy - 1, oz - 1 ) - p0;
1738  diff[ 7 ] = &in( ox + 1, oy , oz - 1 ) - p0;
1739  diff[ 8 ] = &in( ox + 1, oy + 1, oz - 1 ) - p0;
1740 
1741  diff[ 9 ] = &in( ox , oy , oz ) - p0;
1742  diff[ 10 ] = &in( ox , oy + 1, oz ) - p0;
1743  diff[ 11 ] = &in( ox - 1, oy + 1, oz ) - p0;
1744  diff[ 12 ] = &in( ox - 1, oy , oz ) - p0;
1745  diff[ 13 ] = &in( ox - 1, oy - 1, oz ) - p0;
1746  diff[ 14 ] = &in( ox , oy - 1, oz ) - p0;
1747  diff[ 15 ] = &in( ox + 1, oy - 1, oz ) - p0;
1748  diff[ 16 ] = &in( ox + 1, oy , oz ) - p0;
1749  diff[ 17 ] = &in( ox + 1, oy + 1, oz ) - p0;
1750 
1751  diff[ 18 ] = &in( ox , oy , oz + 1 ) - p0;
1752  diff[ 19 ] = &in( ox , oy + 1, oz + 1 ) - p0;
1753  diff[ 20 ] = &in( ox - 1, oy + 1, oz + 1 ) - p0;
1754  diff[ 21 ] = &in( ox - 1, oy , oz + 1 ) - p0;
1755  diff[ 22 ] = &in( ox - 1, oy - 1, oz + 1 ) - p0;
1756  diff[ 23 ] = &in( ox , oy - 1, oz + 1 ) - p0;
1757  diff[ 24 ] = &in( ox + 1, oy - 1, oz + 1 ) - p0;
1758  diff[ 25 ] = &in( ox + 1, oy , oz + 1 ) - p0;
1759  diff[ 26 ] = &in( ox + 1, oy + 1, oz + 1 ) - p0;
1760 
1761  nc6[ 0 ] = diff[ 12 ];
1762  nc6[ 1 ] = diff[ 16 ];
1763  nc6[ 2 ] = diff[ 14 ];
1764  nc6[ 3 ] = diff[ 10 ];
1765  nc6[ 4 ] = diff[ 0 ];
1766  nc6[ 5 ] = diff[ 18 ];
1767  }
1768 
1769  {
1770  const_pointer p0 = &in[ 0 ];
1771 
1772  // Step2 境界画素の検出
1773  for( size_type i = 0 ; i < in.size( ) ; i++ )
1774  {
1775  value_type &v = in[ i ];
1776  if( v > 20 )
1777  {
1778  const_pointer p = &v;
1779  if( p[ nc6[ 0 ] ] == 0 || p[ nc6[ 1 ] ] == 0 ||
1780  p[ nc6[ 2 ] ] == 0 || p[ nc6[ 3 ] ] == 0 ||
1781  p[ nc6[ 4 ] ] == 0 || p[ nc6[ 5 ] ] == 0 )
1782  {
1783  blist.push_back( border_type( p - p0, v ) );
1784  v = 1;
1785  }
1786  }
1787  }
1788  }
1789 
1790  // 以降の処理のためにメモリの予備を確保する
1791  blist.reserve( blist.size( ) * 2 + 1 );
1792  slist[ 0 ].reserve( blist.size( ) );
1793  slist[ 1 ].reserve( blist.size( ) );
1794  slist[ 2 ].reserve( blist.size( ) );
1795  slist[ 3 ].reserve( blist.size( ) );
1796  slist[ 4 ].reserve( blist.size( ) );
1797  slist[ 5 ].reserve( blist.size( ) );
1798  slist[ 6 ].reserve( blist.size( ) );
1799  slist[ 7 ].reserve( blist.size( ) );
1800  slist[ 8 ].reserve( blist.size( ) );
1801 
1802  size_type count = 0, loop = 0;
1803 
1804  do
1805  {
1806  size_type __length__ = 0;
1807 
1808  //Step3 メインサイクル
1809  for( size_type l = 0 ; l < blist.size( ) ; l++ )
1810  {
1811  border_type &b = blist[ l ];
1812  pointer p = &in[ b.diff ];
1813 
1814  if( b.value <= min )
1815  {
1816  create_neighbor_list( p, val, diff );
1817 
1818  // 消去不可能なら永久保存点
1819  if( !neighbor_type::is_deletable( val ) )
1820  {
1821  continue;
1822  }
1823  else
1824  {
1825  // 図形の厚さが１なら永久保存点
1826  if( !neighbor_type::is_3d_simplex( val ) )
1827  {
1828  continue;
1829  }
1830  else
1831  {
1832  // 自分自身を除く
1833  size_type num = val[ 0 ];
1834  num += val[ 1 ];
1835  num += val[ 2 ];
1836  num += val[ 3 ];
1837  num += val[ 4 ];
1838  num += val[ 5 ];
1839  num += val[ 6 ];
1840  num += val[ 7 ];
1841  num += val[ 8 ];
1842  num += val[ 10 ];
1843  num += val[ 11 ];
1844  num += val[ 12 ];
1845  num += val[ 13 ];
1846  num += val[ 14 ];
1847  num += val[ 15 ];
1848  num += val[ 16 ];
1849  num += val[ 17 ];
1850  num += val[ 18 ];
1851  num += val[ 19 ];
1852  num += val[ 20 ];
1853  num += val[ 21 ];
1854  num += val[ 22 ];
1855  num += val[ 23 ];
1856  num += val[ 24 ];
1857  num += val[ 25 ];
1858  num += val[ 26 ];
1859 
1860  difference_type val = num / 3 + 7;
1861  b.value = static_cast< value_type >( val );
1862 
1863  if( 7 <= val && val <= 15 )
1864  {
1865  slist[ val - 7 ].push_back( b );
1866  continue;
1867  }
1868  }
1869  }
1870  }
1871 
1872  blist[ __length__++ ] = b;
1873  }
1874 
1875  blist.erase( blist.begin( ) + __length__, blist.end( ) );
1876 
1877  // Step4 サブサイクル
1878  for( size_type ll = 0 ; ll < 9 ; ll++ )
1879  {
1880  border_list_type &list = slist[ ll ];
1881 
1882  for( size_type l = 0 ; l < list.size( ) ; l++ )
1883  {
1884  border_type &b = list[ l ];
1885  pointer p = &in[ b.diff ];
1886 
1887  create_neighbor_list( p, val, diff );
1888 
1889  //消去不可能なら永久保存点
1890  if( !neighbor_type::is_deletable( val ) )
1891  {
1892  continue;
1893  }
1894  else
1895  {
1896  //図形の厚さが１なら永久保存点
1897  if( !neighbor_type::is_3d_simplex( val ) )
1898  {
1899  continue;
1900  }
1901  else
1902  {
1903  //画素の消去
1904  p[ 0 ] = 0;
1905 
1906  if( p[ nc6[ 0 ] ] > 20 )
1907  {
1908  blist.push_back( border_type( b.diff + nc6[ 0 ], p[ nc6[ 0 ] ] ) );
1909  p[ nc6[ 0 ] ] = 1;
1910  }
1911  if( p[ nc6[ 1 ] ] > 20 )
1912  {
1913  blist.push_back( border_type( b.diff + nc6[ 1 ], p[ nc6[ 1 ] ] ) );
1914  p[ nc6[ 1 ] ] = 1;
1915  }
1916  if( p[ nc6[ 2 ] ] > 20 )
1917  {
1918  blist.push_back( border_type( b.diff + nc6[ 2 ], p[ nc6[ 2 ] ] ) );
1919  p[ nc6[ 2 ] ] = 1;
1920  }
1921  if( p[ nc6[ 3 ] ] > 20 )
1922  {
1923  blist.push_back( border_type( b.diff + nc6[ 3 ], p[ nc6[ 3 ] ] ) );
1924  p[ nc6[ 3 ] ] = 1;
1925  }
1926  if( p[ nc6[ 4 ] ] > 20 )
1927  {
1928  blist.push_back( border_type( b.diff + nc6[ 4 ], p[ nc6[ 4 ] ] ) );
1929  p[ nc6[ 4 ] ] = 1;
1930  }
1931  if( p[ nc6[ 5 ] ] > 20 )
1932  {
1933  blist.push_back( border_type( b.diff + nc6[ 5 ], p[ nc6[ 5 ] ] ) );
1934  p[ nc6[ 5 ] ] = 1;
1935  }
1936  }
1937  }
1938  }
1939 
1940  // 次以降の処理のためにリスト空にする（内部で使用するメモリ容量は変化しない）
1941  list.clear( );
1942  }
1943 
1944  //Step5 終了判定
1945  min = max;
1946  count = 0;
1947  for( size_type l = 0 ; l < blist.size( ) ; l++ )
1948  {
1949  border_type &b = blist[ l ];
1950  if( b.value < min && b.value > 20 )
1951  {
1952  min = b.value;
1953  }
1954  }
1955 
1956  std::cout << loop++ << " \r";
1957 
1958  } while( min < max || blist.size( ) > 0 );
1959 
1960  std::cout << std::endl;
1961 
1962  // Step6 後処理
1963  for( size_type i = 0 ; i < in.size( ) ; i++ )
1964  {
1965  if( in[ i ] != 0 )
1966  {
1967  in[ i ] = 1;
1968  }
1969  }
1970  }
1971 }
1972 
1973 
1974 
1976 namespace hilditch
1977 {
1985 
2002  template < class T1, class T2, class Allocator1, class Allocator2 >
2003  void thinning( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out )
2004  {
2005  typedef typename array2< T2, Allocator2 >::size_type size_type;
2006  typedef typename array2< T2, Allocator2 >::value_type value_type;
2007 
2008  out.resize( in.size1( ), in.size2( ) );
2009  out.reso1( in.reso1( ) );
2010  out.reso2( in.reso2( ) );
2011 
2012  for( size_type i = 0 ; i < in.size( ) ; i++ )
2013  {
2014  out[ i ] = in[ i ] > 0 ? 1 : 0;
2015  }
2016 
2017  __thinning_controller__::thinning( out );
2018  }
2019 
2021  // 細線化グループの終わり
2022 }
2023 
2024 
2026 namespace euclidean
2027 {
2030 
2052  template < class T1, class T2, class Allocator1, class Allocator2 >
2053  void thinning8( const array2< T1, Allocator1 > &in, array2< T2, Allocator2 > &out )
2054  {
2055  typedef typename array2< T2, Allocator2 >::size_type size_type;
2056  typedef typename array2< T2, Allocator2 >::value_type value_type;
2057 
2058  out.resize( in.size1( ), in.size2( ) );
2059 
2060  if( in.reso1( ) <= in.reso2( ) || is_float< value_type >::value )
2061  {
2062  out.reso1( in.reso1( ) );
2063  out.reso2( in.reso2( ) );
2064 
2065  for( size_type i = 0 ; i < in.size( ) ; i++ )
2066  {
2067  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2068  }
2069  __euclidean_utility__::thinning( out );
2070  }
2071  else
2072  {
2073  double r1 = in.reso1( );
2074  double r2 = in.reso2( );
2075 
2076  out.reso1( 1.0 );
2077  out.reso2( 1.0 );
2078 
2079  for( size_type i = 0 ; i < in.size( ) ; i++ )
2080  {
2081  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2082  }
2083  __euclidean_utility__::thinning( out );
2084 
2085  out.reso1( r1 );
2086  out.reso2( r2 );
2087  }
2088  }
2089 
2090 
2101  template < class T1, class T2, class Allocator1, class Allocator2 >
2102  void shrink_skelton6( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
2103  {
2104  typedef typename array3< T2, Allocator2 >::size_type size_type;
2105  typedef typename array3< T2, Allocator2 >::value_type value_type;
2106 
2107  out.resize( in.size1( ), in.size2( ), in.size3( ) );
2108 
2109  out.reso1( in.reso1( ) );
2110  out.reso2( in.reso2( ) );
2111  out.reso3( in.reso3( ) );
2112 
2113  for( size_type i = 0 ; i < in.size( ) ; i++ )
2114  {
2115  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2116  }
2117  __euclidean_utility__::shrink_skelton( out, __euclidean_utility__::neighbor< 6 >( ) );
2118  }
2119 
2120 
2131  template < class T1, class T2, class Allocator1, class Allocator2 >
2132  void shrink_skelton26( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
2133  {
2134  typedef typename array3< T2, Allocator2 >::size_type size_type;
2135  typedef typename array3< T2, Allocator2 >::value_type value_type;
2136 
2137  out.resize( in.size1( ), in.size2( ), in.size3( ) );
2138 
2139  out.reso1( in.reso1( ) );
2140  out.reso2( in.reso2( ) );
2141  out.reso3( in.reso3( ) );
2142 
2143  for( size_type i = 0 ; i < in.size( ) ; i++ )
2144  {
2145  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2146  }
2147  __euclidean_utility__::shrink_skelton( out, __euclidean_utility__::neighbor< 26 >( ) );
2148  }
2149 
2150 
2171  template < class T1, class T2, class Allocator1, class Allocator2 >
2172  void thinning6( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
2173  {
2174  typedef typename array3< T2, Allocator2 >::size_type size_type;
2175  typedef typename array3< T2, Allocator2 >::value_type value_type;
2176 
2177  out.resize( in.size1( ), in.size2( ), in.size3( ) );
2178 
2179  if( ( in.reso1( ) <= in.reso2( ) && in.reso1( ) <= in.reso3( ) ) || is_float< value_type >::value )
2180  {
2181  out.reso1( in.reso1( ) );
2182  out.reso2( in.reso2( ) );
2183  out.reso3( in.reso3( ) );
2184 
2185  for( size_type i = 0 ; i < in.size( ) ; i++ )
2186  {
2187  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2188  }
2189  __euclidean_utility__::thinning( out, __euclidean_utility__::neighbor< 6 >( ) );
2190  }
2191  else
2192  {
2193  double r1 = in.reso1( );
2194  double r2 = in.reso2( );
2195  double r3 = in.reso3( );
2196 
2197  out.reso1( 1.0 );
2198  out.reso2( 1.0 );
2199  out.reso3( 1.0 );
2200 
2201  for( size_type i = 0 ; i < in.size( ) ; i++ )
2202  {
2203  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2204  }
2205  __euclidean_utility__::thinning( out, __euclidean_utility__::neighbor< 6 >( ) );
2206 
2207  out.reso1( r1 );
2208  out.reso2( r2 );
2209  out.reso3( r3 );
2210  }
2211  }
2212 
2213 
2234  template < class T1, class T2, class Allocator1, class Allocator2 >
2235  void thinning26( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
2236  {
2237  typedef typename array3< T2, Allocator2 >::size_type size_type;
2238  typedef typename array3< T2, Allocator2 >::value_type value_type;
2239 
2240  out.resize( in.size1( ), in.size2( ), in.size3( ) );
2241 
2242  if( ( in.reso1( ) <= in.reso2( ) && in.reso1( ) <= in.reso3( ) ) || is_float< value_type >::value )
2243  {
2244  out.reso1( in.reso1( ) );
2245  out.reso2( in.reso2( ) );
2246  out.reso3( in.reso3( ) );
2247 
2248  for( size_type i = 0 ; i < in.size( ) ; i++ )
2249  {
2250  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2251  }
2252  __euclidean_utility__::thinning( out, __euclidean_utility__::neighbor< 26 >( ) );
2253  }
2254  else
2255  {
2256  double r1 = in.reso1( );
2257  double r2 = in.reso2( );
2258  double r3 = in.reso3( );
2259 
2260  out.reso1( 1.0 );
2261  out.reso2( 1.0 );
2262  out.reso3( 1.0 );
2263 
2264  for( size_type i = 0 ; i < in.size( ) ; i++ )
2265  {
2266  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2267  }
2268  __euclidean_utility__::thinning( out, __euclidean_utility__::neighbor< 26 >( ) );
2269 
2270  out.reso1( r1 );
2271  out.reso2( r2 );
2272  out.reso3( r3 );
2273  }
2274  }
2275 
2292  template < class T1, class T2, class Allocator1, class Allocator2 >
2293  void surface_thinning6( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
2294  {
2295  typedef typename array3< T2, Allocator2 >::size_type size_type;
2296  typedef typename array3< T2, Allocator2 >::value_type value_type;
2297 
2298  out.resize( in.size1( ), in.size2( ), in.size3( ) );
2299 
2300  if( ( in.reso1( ) <= in.reso2( ) && in.reso1( ) <= in.reso3( ) ) || is_float< value_type >::value )
2301  {
2302  out.reso1( in.reso1( ) );
2303  out.reso2( in.reso2( ) );
2304  out.reso3( in.reso3( ) );
2305 
2306  for( size_type i = 0 ; i < in.size( ) ; i++ )
2307  {
2308  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2309  }
2310  __euclidean_utility__::surface_thinning( out, __euclidean_utility__::neighbor< 6 >( ) );
2311  }
2312  else
2313  {
2314  double r1 = in.reso1( );
2315  double r2 = in.reso2( );
2316  double r3 = in.reso3( );
2317 
2318  out.reso1( 1.0 );
2319  out.reso2( 1.0 );
2320  out.reso3( 1.0 );
2321 
2322  for( size_type i = 0 ; i < in.size( ) ; i++ )
2323  {
2324  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2325  }
2326  __euclidean_utility__::surface_thinning( out, __euclidean_utility__::neighbor< 6 >( ) );
2327 
2328  out.reso1( r1 );
2329  out.reso2( r2 );
2330  out.reso3( r3 );
2331  }
2332  }
2333 
2334 
2351  template < class T1, class T2, class Allocator1, class Allocator2 >
2352  void surface_thinning26( const array3< T1, Allocator1 > &in, array3< T2, Allocator2 > &out )
2353  {
2354  typedef typename array3< T2, Allocator2 >::size_type size_type;
2355  typedef typename array3< T2, Allocator2 >::value_type value_type;
2356 
2357  out.resize( in.size1( ), in.size2( ), in.size3( ) );
2358 
2359  if( ( in.reso1( ) <= in.reso2( ) && in.reso1( ) <= in.reso3( ) ) || is_float< value_type >::value )
2360  {
2361  out.reso1( in.reso1( ) );
2362  out.reso2( in.reso2( ) );
2363  out.reso3( in.reso3( ) );
2364 
2365  for( size_type i = 0 ; i < in.size( ) ; i++ )
2366  {
2367  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2368  }
2369  __euclidean_utility__::surface_thinning( out, __euclidean_utility__::neighbor< 26 >( ) );
2370  }
2371  else
2372  {
2373  double r1 = in.reso1( );
2374  double r2 = in.reso2( );
2375  double r3 = in.reso3( );
2376 
2377  out.reso1( 1.0 );
2378  out.reso2( 1.0 );
2379  out.reso3( 1.0 );
2380 
2381  for( size_type i = 0 ; i < in.size( ) ; i++ )
2382  {
2383  out[ i ] = static_cast< value_type >( in[ i ] > 0 ? 1 : 0 );
2384  }
2385  __euclidean_utility__::surface_thinning( out, __euclidean_utility__::neighbor< 26 >( ) );
2386 
2387  out.reso1( r1 );
2388  out.reso2( r2 );
2389  out.reso3( r3 );
2390  }
2391  }
2392 
2394  // 細線化グループの終わり
2395 }
2396 
2397 
2398 
2399 // mist名前空間の終わり
2400 _MIST_END
2401 
2402 
2403 #endif // __INCLUDE_MIST_THINNING__
2404 
Generated on Wed Nov 12 2014 19:44:25 for MIST by doxygen 1.8.1.2