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

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

decomposition.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 #ifndef __INCLUDE_MIST_FIGURE_DECOMPOSITION__
34 #define __INCLUDE_MIST_FIGURE_DECOMPOSITION__
35
36
37 #ifndef __INCLUDE_MIST_H__
38 #include "../mist.h"
39 #endif
40
41 #ifndef __INCLUDE_MIST_LIMITS__
42 #include "../limits.h"
43 #endif
44
45 #ifndef __INCLUDE_MIST_DISTANCE_TRANSFORM__
46 #include "distance.h"
47 #endif
48
49 #ifndef __INCLUDE_MIST_LABELING__
50 #include "labeling.h"
51 #endif
52
53
54 #include <deque>
55 #include <vector>
56 #include <set>
57 #include <algorithm>
58
59
60 // mist名前空間の始まり
61 _MIST_BEGIN
62
63
64 // 図形分割用のツール
65 namespace __figure_dedomposition__
66 {
67 struct Position
68 {
69 typedef ptrdiff_t difference_type;
70 difference_type RADIUS;
71 double radius;
72 difference_type index;
73
74 Position( difference_type r = 0, double R = 0, difference_type indx = 0 ) : RADIUS( r ), radius( R - std::sqrt( static_cast< double >( r ) ) ), index( indx )
75 {
76 }
77
78 bool operator <( const Position &p ) const
79 {
80 return( RADIUS > p.RADIUS );
81 }
82 };
83
84 template < class T, class Allocator >
85 inline typename array2< T, Allocator >::size_type do_labeling( array2< T, Allocator > &in )
86 {
87 return( labeling8( in, in ) );
88 }
89
90 template < class T, class Allocator >
91 inline typename array3< T, Allocator >::size_type do_labeling( array3< T, Allocator > &in )
92 {
93 return( labeling26( in, in ) );
94 }
95
96 // 画像の端に画素を持っているかどうかを調べる
97 template < class T, class Allocator >
98 bool has_voxel_at_side( const array2< T, Allocator > &in )
99 {
100 typedef typename array2< T, Allocator >::size_type size_type;
101
102 for( size_type i = 0 ; i < in.width( ) ; i++ )
103 {
104 if( in( i, 0 ) != 0 || in( i, in.height( ) - 1 ) != 0 )
105 {
106 return( true );
107 }
108 }
109
110 for( size_type j = 0 ; j < in.height( ) ; j++ )
111 {
112 if( in( 0, j ) != 0 || in( in.width( ) - 1, j ) != 0 )
113 {
114 return( true );
115 }
116 }
117
118 return( false );
119 }
120
121 // 画像の端に画素を持っているかどうかを調べる
122 template < class T, class Allocator >
123 bool has_voxel_at_side( const array3< T, Allocator > &in )
124 {
125 typedef typename array2< T, Allocator >::size_type size_type;
126
127 if( in.depth( ) == 1 )
128 {
129 for( size_type i = 0 ; i < in.width( ) ; i++ )
130 {
131 if( in( i, 0, 0 ) != 0 || in( i, in.height( ) - 1, 0 ) != 0 )
132 {
133 return( true );
134 }
135 }
136
137 for( size_type j = 0 ; j < in.height( ) ; j++ )
138 {
139 if( in( 0, j, 0 ) != 0 || in( in.width( ) - 1, j, 0 ) != 0 )
140 {
141 return( true );
142 }
143 }
144 }
145 else
146 {
147 for( size_type j = 0 ; j < in.height( ) ; j++ )
148 {
149 for( size_type i = 0 ; i < in.width( ) ; i++ )
150 {
151 if( in( i, j, 0 ) != 0 || in( i, j, in.depth( ) - 1 ) != 0 )
152 {
153 return( true );
154 }
155 }
156 }
157
158 for( size_type k = 0 ; k < in.depth( ) ; k++ )
159 {
160 for( size_type i = 0 ; i < in.width( ) ; i++ )
161 {
162 if( in( i, 0, k ) != 0 || in( i, in.height( ) - 1, k ) != 0 )
163 {
164 return( true );
165 }
166 }
167 }
168
169 for( size_type k = 0 ; k < in.depth( ) ; k++ )
170 {
171 for( size_type j = 0 ; j < in.height( ) ; j++ )
172 {
173 if( in( 0, j, k ) != 0 || in( in.width( ) - 1, j, k ) != 0 )
174 {
175 return( true );
176 }
177 }
178 }
179 }
180
181 return( false );
182 }
183
184 struct __mist_dmy_fd_callback__
185 {
194 template < class Array1, class Array2 >
195 void operator()( size_t loop, size_t label_num, double radius, const Array1 &in, const Array2 &out ) const
196 {
197 }
198 };
199 }
200
201
209
210
211
212
224 template < class Array1, class Array2, class Functor >
225 typename Array1::size_type figure_decomposition( const Array1 &in, Array2 &out, double max_distance, Functor f )
226 {
227 if( in.empty( ) )
228 {
229 return( 0 );
230 }
231
232 typedef typename Array1::size_type size_type;
233 typedef typename Array1::difference_type difference_type;
234 typedef typename Array1::const_pointer const_pointer;
235 typedef typename Array2::value_type value_type;
236 typedef typename Array2::template rebind< size_type >::other distance_type;
237 typedef typename Array2::template rebind< double >::other mask_type;
238
239 distance_type dist;
240 dist.resize( in.width( ), in.height( ), in.depth( ) );
241 dist.reso1( 1 );
242 dist.reso2( 1 );
243 dist.reso3( 1 );
244
245 // 距離変換用の画像を作成
246 for( size_type i = 0 ; i < in.size( ) ; i++ )
247 {
248 dist[ i ] = in[ i ] != 0 ? 1 : 0;
249 }
250
251 if( __figure_dedomposition__::has_voxel_at_side( in ) )
252 {
253 // 画像の端に１画素がある場合
254 marray< distance_type > dist_tmp( dist, 1 );
255 distance_type &dist_tmp_map = dist_tmp;
256
257 // 画像の外側には０画素が存在するものとして距離変換を行う
258 calvin::distance_transform( dist_tmp_map, dist_tmp_map );
259
260 for( size_type k = 0 ; k < dist.depth( ) ; k++ )
261 {
262 for( size_type j = 0 ; j < dist.height( ) ; j++ )
263 {
264 for( size_type i = 0 ; i < dist.width( ) ; i++ )
265 {
266 dist( i, j, k ) = dist_tmp( i, j, k );
267 }
268 }
269 }
270 }
271 else
272 {
273 calvin::distance_transform( dist, dist );
274 }
275
276
277 typedef __figure_dedomposition__::Position position_type;
278
279 // 画像中に含まれる距離値のリストを作成する
280 std::map< size_type, size_type > distance_list;
281 for( size_type l = 0 ; l < dist.size( ) ; l++ )
282 {
283 distance_list[ dist[ l ] ]++;
284 }
285
286 size_type label_count = 0; // 現在のラベル数
287 size_type loop_count = 0; // 現在のラベル数
288 size_type current_label = 0; // 現在処理中のラベル
289 size_type label_max = type_limits< value_type >::maximum( ); // ラベルの最大値
290 position_type pt;
291
292 mask_type mask( in );
293 mask.fill( );
294 out.resize( in.width( ), in.height( ), in.depth( ) );
295 out.reso1( in.reso1( ) );
296 out.reso2( in.reso2( ) );
297 out.reso3( in.reso3( ) );
298 out.fill( );
299
300 typename std::map< size_type, size_type >::reverse_iterator dite = distance_list.rbegin( );
301
302 // 最大距離値より大きい部分は処理対象から省く
303 if( max_distance > 0 )
304 {
305 double md = max_distance * max_distance;
306 double Md = ( max_distance + 1 ) * ( max_distance + 1 );
307 for( size_type l = 0 ; l < dist.size( ) ; l++ )
308 {
309 if( dist[ l ] >= md )
310 {
311 out[ l ] = 1;
312 }
313 }
314 for( ; dite != distance_list.rend( ) ; ++dite )
315 {
316 if( dite->first < Md )
317 {
318 break;
319 }
320 }
321
322 label_count = __figure_dedomposition__::do_labeling( out );
323 }
324
325 f( 0, label_count, std::sqrt( static_cast< double >( dite->first ) ), in, out );
326
327 // 距離値の大きいものから順番にラベルを復元していく
328 for( ; dite != distance_list.rend( ) ; ++dite )
329 {
330 size_type rr = dite->first;
331 difference_type count = dite->second;
332
333 if( rr == 0 )
334 {
335 // ０の画素はアクセスしない
336 break;
337 }
338
339 double r = std::sqrt( static_cast< double >( rr ) );
340
341 std::vector< difference_type > list;
342 list.reserve( count );
343 for( size_type i = 0 ; i < dist.size( ) ; i++ )
344 {
345 // 現段階で最大の距離値を持つ画素を探す
346 if( dist[ i ] == rr && r >= mask[ i ] )
347 //if( dist[ i ] == rr && r >= static_cast< difference_type >( mask[ i ] + 1 ) )
348 {
349 list.push_back( i );
350 }
351 }
352
353
354 difference_type rx = static_cast< difference_type >( std::ceil( r ) );
355 difference_type ry = in.height( ) < 2 ? 0 : rx;
356 difference_type rz = in.depth( ) < 2 ? 0 : rx;
357 size_type RR = ( rx + 1 ) * ( rx + 1 );
358
359 if( rr > 1 )
360 {
361 std::vector< position_type > sphere;
362 sphere.reserve( ( 2 * rx + 1 ) * ( 2 * ry + 1 ) * ( 2 * rz + 1 ) );
363 {
364 difference_type cx = out.width( ) / 2;
365 difference_type cy = out.height( ) / 2;
366 difference_type cz = out.depth( ) / 2;
367 const_pointer cp = &out( cx, cy, cz );
368 for( difference_type k = -rz ; k <= rz ; k++ )
369 {
370 size_type kk = k * k;
371
372 for( difference_type j = -ry ; j <= ry ; j++ )
373 {
374 size_type jj = j * j;
375
376 for( difference_type i = -rx ; i <= rx ; i++ )
377 {
378 size_type ii = i * i;
379
380 size_type rrr = ii + jj + kk;
381 if( rrr < rr )
382 {
383 // インデックスとして，球の中心からの差を代入する
384 sphere.push_back( position_type( rrr, r, &out( cx + i, cy + j, cz + k ) - cp ) );
385 }
386 }
387 }
388 }
389 }
390
391 difference_type csize = list.size( );
392 while( csize > 0 )
393 {
394 difference_type cindex = 0;
395 for( ; cindex < csize ; cindex++ )
396 {
397 if( out[ list[ cindex ] ] != 0 )
398 {
399 break;
400 }
401 }
402
403 if( cindex == csize )
404 {
405 cindex = 0;
406 }
407
408 difference_type index = list[ cindex ];
409 list[ cindex ] = list[ csize - 1 ];
410 csize--;
411
412 if( out[ index ] == 0 )
413 {
414 // 他の領域から塗られていないので，新しいラベルとする
415 label_count++;
416 if( label_count > label_max )
417 {
418 // 最大のラベル数を超えた場合には，最大ラベルとする
419 label_count = label_max + 1;
420 }
421 current_label = static_cast< value_type >( label_count );
422 }
423 else
424 {
425 // 既にラベルが割り当てられているので，そのラベルで塗りつぶす
426 current_label = out[ index ];
427 }
428
429 for( size_type i = 0 ; i < sphere.size( ) ; i++ )
430 {
431 const position_type &pt = sphere[ i ];
432 const difference_type indx = index + pt.index;
433 if( mask[ indx ] < pt.radius )
434 {
435 if( dist[ indx ] + pt.RADIUS < RR )
436 {
437 mask[ indx ] = pt.radius;
438 out[ indx ] = static_cast< value_type >( current_label );
439 }
440 }
441 }
442 }
443 }
444 else
445 {
446 // 距離値１の画素だけ特殊な処理を行う
447 difference_type csize = list.size( );
448 while( csize > 0 )
449 {
450 current_label = 0;
451
452 difference_type cindex = 0;
453 for( ; cindex < csize ; cindex++ )
454 {
455 difference_type indx = list[ cindex ];
456 difference_type z = indx / ( out.width( ) * out.height( ) );
457 difference_type y = ( indx - z * out.width( ) * out.height( ) ) / out.width( );
458 difference_type x = indx - ( y + z * out.height( ) ) * out.width( );
459
460 // ８・２６近傍から既に塗られているものを探してくる
461 double min = 1.0e10;
462 for( difference_type k = -rz ; k <= rz ; k++ )
463 {
464 size_type pz = k + z;
465 if( pz < 0 || pz >= in.depth( ) ) continue;
466 double kk = static_cast< double >( k * k );
467
468 for( difference_type j = -ry ; j <= ry ; j++ )
469 {
470 size_type py = j + y;
471 if( py < 0 || py >= in.height( ) ) continue;
472 double jj = static_cast< double >( j * j );
473
474 for( difference_type i = -rx ; i <= rx ; i++ )
475 {
476 size_type px = i + x;
477 if( px < 0 || px >= in.width( ) ) continue;
478 double rR = static_cast< double >( i * i + jj + kk );
479
480 size_type cl = static_cast< size_type >( out( px, py, pz ) );
481 if( cl != 0 && min > rR )
482 {
483 current_label = cl;
484 min = rR;
485 }
486 }
487 }
488 }
489 if( current_label != 0 )
490 {
491 break;
492 }
493 }
494
495 if( cindex == csize )
496 {
497 cindex = 0;
498 }
499
500 difference_type index = list[ cindex ];
501 list[ cindex ] = list[ csize - 1 ];
502 csize--;
503
504 if( current_label == 0 )
505 {
506 // 他の領域から塗られていないので，新しいラベルとする
507 label_count++;
508 if( label_count > label_max )
509 {
510 // 最大のラベル数を超えた場合には，最大ラベルとする
511 label_count = label_max + 1;
512 }
513 out[ index ] = static_cast< value_type >( label_count );
514 }
515 else
516 {
517 // 既にラベルが割り当てられているので，そのラベルで塗りつぶす
518 out[ index ] = static_cast< value_type >( current_label );
519 }
520 }
521 }
522
523 f( ++loop_count, label_count, r, in, out );
524 }
525
526 // MAXラベル数を超えたものを除去
527 for( size_type j = 0 ; j < out.size( ) ; j++ )
528 {
529 out[ j ] = static_cast< size_type >( out[ j ] ) > label_max ? 0 : out[ j ];
530 }
531
532 return( label_count );
533 }
534
545 template < class Array1, class Array2 >
546 typename Array1::size_type figure_decomposition( const Array1 &in, Array2 &out, double max_distance = -1 )
547 {
548 return( figure_decomposition( in, out, max_distance, __figure_dedomposition__::__mist_dmy_fd_callback__( ) ) );
549 }
550
551
553 // 図形分割グループの終わり
554
555
556 // mist名前空間の終わり
557 _MIST_END
558
559
560 #endif // __INCLUDE_MIST_FIGURE_DECOMPOSITION__
561