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

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

mesh.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_MESH__
34 #define __INCLUDE_MIST_MESH__
35
36
37 #ifndef __INCLUDE_MIST_CONF_H__
38 #include "../config/mist_conf.h"
39 #endif
40
41 #ifndef __INCLUDE_MIST_MATRIX__
42 #include "../matrix.h"
43 #endif
44
45 #ifndef __INCLUDE_MIST_VECTOR__
46 #include "../vector.h"
47 #endif
48
49 #ifndef __INCLUDE_CONVERTER__
50 #include "../converter.h"
51 #endif
52
53 #ifndef __INCLUDE_MIST_LABELING__
54 #include "../filter/labeling.h"
55 #endif
56
57 #ifndef __INCLUDE_MIST_DRAWING__
58 #include "../drawing.h"
59 #endif
60
61
62 #include <deque>
63 #include <set>
64 #include <algorithm>
65
66 // mist名前空間の始まり
67 _MIST_BEGIN
68
69
77
78 namespace __mesh_utility__
79 {
80 struct __mesh_information__
81 {
82 typedef size_t size_type;
83 typedef ptrdiff_t difference_type;
84 typedef vector2< double > vector_type;
85
86 difference_type label;
87 difference_type count;
88 double round;
89 vector_type pos;
90 bool enabled;
91 difference_type mark;
92 difference_type row;
93 difference_type col;
94 double length;
95
96 __mesh_information__( ) : label( -1 ), count( 0 ), round( 0.0 ), pos( 0.0, 0.0 ), enabled( true ), mark( -1 ), row( -1 ), col( -1 ), length( 1.0e10 )
97 {
98 }
99
100 // 周囲長の大きさで並べ替える
101 bool operator <( const __mesh_information__ &m ) const
102 {
103 return( this->round > m.round );
104 }
105
106 static bool sort_by_round( const __mesh_information__ &m1, const __mesh_information__ &m2 )
107 {
108 return( m1.round > m2.round );
109 }
110
111 static bool sort_by_length( const __mesh_information__ &m1, const __mesh_information__ &m2 )
112 {
113 return( m1.length < m2.length );
114 }
115 };
116
117 inline double __minimum__( double v1, double v2 )
118 {
119 return( v1 < v2 ? v1 : v2 );
120 }
121
122 inline double __compute_border_distance__( const vector2< double > &pos, double image_width, double image_height )
123 {
124 double xmin = __minimum__( std::abs( pos.x ), std::abs( image_width - pos.x ) );
125 double ymin = __minimum__( std::abs( pos.y ), std::abs( image_height - pos.y ) );
126 return( __minimum__( xmin, ymin ) );
127 }
128 }
129
147 template < class T, class Allocator >
148 bool extract_mesh( array2< T, Allocator > &chart, matrix< vector2< double > > &grid, typename array2< T, Allocator >::difference_type row, typename array2< T, Allocator >::difference_type col, double threshold_for_circular_ratio = 0.4 )
149 {
150 typedef __mesh_utility__::__mesh_information__ __mesh_information__;
151 typedef array2< T, Allocator >::size_type size_type;
152 typedef array2< T, Allocator >::difference_type difference_type;
153 typedef vector2< double > vector_type;
154 typedef std::deque< __mesh_information__ > mesh_list_type;
155 typedef typename mesh_list_type::iterator mesh_iterator;
156 typedef typename mesh_list_type::reverse_iterator mesh_reverse_iterator;
157 array2< size_type > binary;
158
159
160 // ２値画像に変換
161 convert( chart, binary );
162
163 size_type i, j;
164 difference_type r, c;
165 difference_type rows = grid.rows( );
166 difference_type cols = grid.cols( );
167
168 size_type labelnum = labeling4( binary, binary );
169 //std::cout << "LabelNum: " << labelnum << std::endl;
170
171 mesh_list_type mesh_list;
172 array< __mesh_information__ > meshes( labelnum + 1 );
173
174 for( i = 0 ; i <= labelnum ; i++ )
175 {
176 meshes[ i ].label = i;
177 }
178
179 // 全てのグリッドを初期化する
180 for( size_type i = 0 ; i < grid.size( ) ; i++ )
181 {
182 grid[ i ].x = -1;
183 grid[ i ].y = -1;
184 }
185
186 // 背景は無効にする
187 meshes[ 0 ].enabled = false;
188
189 // 各ラベルの画素数，重心位置，周囲長を求める
190 const double _2 = std::sqrt( 2.0 );
191 for( j = 0 ; j < chart.height( ) - 1 ; j++ )
192 {
193 for( i = 0 ; i < chart.width( ) - 1 ; i++ )
194 {
195 __mesh_information__ &mesh = meshes[ binary( i, j ) ];
196 mesh.count++;
197 mesh.pos.x += i;
198 mesh.pos.y += j;
199
200 // 画像の周囲長を計算する
201 __mesh_information__ &m0 = mesh;
202 __mesh_information__ &m1 = meshes[ binary( i + 1, j ) ];
203 __mesh_information__ &m2 = meshes[ binary( i , j + 1 ) ];
204 __mesh_information__ &m3 = meshes[ binary( i + 1, j + 1 ) ];
205 size_type l0 = m0.label;
206 size_type l1 = m1.label;
207 size_type l2 = m2.label;
208 size_type l3 = m3.label;
209
210 if( l0 > 0 )
211 {
212 size_type num_eq = 0;
213 if( l1 == l0 ) num_eq++;
214 if( l2 == l0 ) num_eq++;
215 if( l3 == l0 ) num_eq++;
216
217 if( num_eq == 2 )
218 {
219 m0.round += _2;
220 }
221 else if( num_eq == 1 && l0 != l3 )
222 {
223 m0.round += 1.0;
224 }
225 }
226 else if( l1 > 0 )
227 {
228 if( l1 == l3 )
229 {
230 if( l1 == l2 )
231 {
232 m1.round += _2;
233 }
234 else
235 {
236 m1.round += 1.0;
237 }
238 }
239 }
240 else if( l2 > 0 )
241 {
242 if( l2 == l3 )
243 {
244 m2.round += 1.0;
245 }
246 }
247 }
248 }
249
250
251 // 各領域の重心位置を計算する
252 for( i = 1 ; i <= labelnum ; i++ )
253 {
254 __mesh_information__ &mesh = meshes[ i ];
255 mesh.pos.x /= mesh.count;
256 mesh.pos.y /= mesh.count;
257 }
258
259 // 画面の境界に接している領域を取り除く
260 for( i = 0 ; i < chart.width( ) ; i++ )
261 {
262 meshes[ binary( i, 0 ) ].enabled = false;
263 meshes[ binary( i, chart.height( ) - 1 ) ].enabled = false;
264 }
265
266 // 画面の境界に接している領域を取り除く
267 for( j = 0 ; j < chart.height( ) ; j++ )
268 {
269 meshes[ binary( 0, j ) ].enabled = false;
270 meshes[ binary( chart.width( ) - 1, j ) ].enabled = false;
271 }
272
273 // 円形度を評価して，不要な領域を除去する
274 difference_type minimum_count = 15; // この個数に満たない領域は削除する
275 for( i = 1 ; i <= labelnum ; i++ )
276 {
277 __mesh_information__ &mesh = meshes[ i ];
278 if( mesh.count < minimum_count )
279 {
280 mesh.enabled = false;
281 }
282 else if( mesh.enabled )
283 {
284 const double pai = 3.1415926535897932384626433832795;
285 const double S = mesh.count;
286 const double L = mesh.round;
287 double e = 4.0 * pai * S / ( L * L );
288 if( e < threshold_for_circular_ratio )
289 {
290 mesh.enabled = false;
291 }
292 else
293 {
294 mesh_list.push_back( mesh );
295 }
296 //std::cout << "ラベル: " << i << ", 円形度: " << e << std::endl;
297 }
298 }
299
300 if( mesh_list.size( ) < 10 )
301 {
302 return( false );
303 }
304
305 // 周囲長の大きさを基準として並べ替える
306 std::sort( mesh_list.begin( ), mesh_list.end( ), __mesh_information__::sort_by_round );
307
308 // 周囲長の大き差を基準として，基準点列を抽出する
309 {
310 mesh_iterator ite = mesh_list.begin( );
311 difference_type base_mesh_num = 1;
312 double oround = ite->round;
313 double maxround = ite->round;
314
315 ite->mark = 1;
316 ++ite;
317
318 for( ; ite != mesh_list.end( ) ; ++ite )
319 {
320 double round = ite->round;
321
322 //std::cout << "周囲長: " << round << ", 比: " << round / oround << std::endl;
323
324 // 一つ前の周囲長の90％以上で，最大の円の周囲長の80%以上の長さがある場合，基準点としてマークする
325 if( round > oround * 0.9 && round > maxround * 0.5 )
326 {
327 ite->mark = 1;
328 base_mesh_num++;
329 oround = round;
330 }
331 else
332 {
333 break;
334 }
335 }
336 for( ; ite != mesh_list.end( ) ; ++ite )
337 {
338 ite->mark = 100;
339 }
340
342 // 基準点の数が少ない場合には終了
343 if( base_mesh_num < 2 )
344 {
345 for( ite = mesh_list.begin( ) ; ite != mesh_list.end( ) ; ++ite )
346 {
347 meshes[ ite->label ].mark = ite->mark;
348 }
349
350 for( i = 0 ; i < chart.size( ) ; i++ )
351 {
352 if( meshes[ binary[ i ] ].mark >= 0 && meshes[ binary[ i ] ].enabled )
353 {
354 chart[ i ] = static_cast< unsigned char >( meshes[ binary[ i ] ].mark );
355 }
356 else
357 {
358 chart[ i ] = 255;
359 }
360 }
361
362 return( false );
363 }
364
365 // まず，大きい基準点の座標を用いて基準点を結ぶ直線を求める
366 double x1 = mesh_list[ 0 ].pos.x;
367 double y1 = mesh_list[ 0 ].pos.y;
368 double x2 = mesh_list[ 1 ].pos.x;
369 double y2 = mesh_list[ 1 ].pos.y;
370 double A = y2 - y1;
371 double B = x1 - x2;
372 double C = x2 * y1 - x1 * y2;
373
374 mesh_iterator site = mesh_list.begin( );
375 mesh_iterator eite = site + base_mesh_num;
376
377 // 判定ミスの基準点を排除する
378 {
379 double __threshold__ = std::sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) ) * 0.2;
380 double D = 1.0 / std::sqrt( A * A + B * B );
381 for( ite = site ; ite != eite ; ++ite )
382 {
383 if( std::abs( A * ite->pos.x + B * ite->pos.y + C ) * D > __threshold__ )
384 {
385 ite->mark = 100;
386 base_mesh_num--;
387 }
388 }
389 }
390
391 // 基準点を順番に並べ替える評価値を計算する
392 for( ite = site ; ite != eite ; ++ite )
393 {
394 if( ite->mark == 1 )
395 {
396 ite->length = - B * ite->pos.x + A * ite->pos.y;
397 }
398 else
399 {
400 ite->length = 1.0e10;
401 }
402 }
403
404 // 基準線に沿って並び替える
405 std::sort( site, eite, __mesh_information__::sort_by_length );
406 site = mesh_list.begin( );
407 eite = site + base_mesh_num;
408
409 // 大きさが変化する位置を再度検出する
410 {
411 // まず，もっとも周囲長の長いものを見つける
412 mesh_iterator mite = mesh_list.begin( );
413 mesh_iterator oite = mesh_list.begin( );
414 mesh_iterator cite = oite;
415 for( ; cite != eite ; ++cite )
416 {
417 if( mite->round < cite->round )
418 {
419 mite = cite;
420 }
421 }
422
423 for( cite = mite, oite = mite ; cite != eite ; ++cite )
424 {
425 double r1 = oite->round;
426 double r2 = cite->round;
427
428 // 一つ前の周囲長の80％以上であれば，基準点としてマークする
429 //std::cout << "周囲長: " << r1 << ", 比: " << r2 / r1 << std::endl;
430 if( ( r1 > r2 && r2 < r1 * 0.8 ) || ( r1 < r2 && r1 < r2 * 0.8 ) )
431 {
432 break;
433 }
434
435 oite = cite;
436 }
437
438 for( ; cite != eite ; ++cite )
439 {
440 cite->mark = 100;
441 cite->length = 1.0e10;
442 }
443
444 mesh_reverse_iterator rmite( mite + 1 );
445 mesh_reverse_iterator reite = mesh_list.rend( );
446 mesh_reverse_iterator rcite( rmite );
447 for( mesh_reverse_iterator roite = rmite ; rcite != reite ; ++rcite )
448 {
449 double r1 = roite->round;
450 double r2 = rcite->round;
451
452 // 一つ前の周囲長の80％以上であれば，基準点としてマークする
453 //std::cout << "周囲長: " << r1 << ", 比: " << r2 / r1 << std::endl;
454 if( ( r1 > r2 && r2 < r1 * 0.8 ) || ( r1 < r2 && r1 < r2 * 0.8 ) )
455 {
456 break;
457 }
458
459 roite = rcite;
460 }
461
462 for( ; rcite != reite ; ++rcite )
463 {
464 rcite->mark = 100;
465 rcite->length = 1.0e10;
466 }
467
468 for( cite = site, base_mesh_num = 0 ; cite != eite ; ++cite )
469 {
470 if( cite->mark == 1 )
471 {
472 base_mesh_num++;
473 }
474 }
475
476 // 基準線に沿って並び替える
477 std::sort( site, eite, __mesh_information__::sort_by_length );
478 site = mesh_list.begin( );
479 eite = site + base_mesh_num;
480
481 // 基準点が抽出されすぎたときは終了する
482 if( col + 1 < base_mesh_num )
483 {
484 for( ite = mesh_list.begin( ) ; ite != mesh_list.end( ) ; ++ite )
485 {
486 meshes[ ite->label ].mark = ite->mark;
487 }
488
489 for( i = 0 ; i < chart.size( ) ; i++ )
490 {
491 if( meshes[ binary[ i ] ].mark >= 0 && meshes[ binary[ i ] ].enabled )
492 {
493 chart[ i ] = static_cast< unsigned char >( meshes[ binary[ i ] ].mark );
494 }
495 else
496 {
497 chart[ i ] = 255;
498 }
499 }
500
501 return( false );
502 }
503 }
504
505 // 中心をあらわす基準点を決定する
506 // 基準ラインの左右の点のどちらかを基準にする
507 {
508 mesh_iterator ite1 = mesh_list.begin( );
509 mesh_iterator ite2 = ite1 + base_mesh_num - 1;
510 double lth = 2.0 * ( ite1->pos - ite2->pos ).length( ) / static_cast< double >( base_mesh_num - 1 );
511
512 double ratio1 = 1.0e10;
513 double ratio2 = 1.0e10;
514
515 double len0 = __mesh_utility__::__compute_border_distance__( mesh_list[ 0 ].pos, chart.width( ), chart.height( ) );
516 double len1 = __mesh_utility__::__compute_border_distance__( mesh_list[ base_mesh_num - 1 ].pos, chart.width( ), chart.height( ) );
517
518 // 端から一定距離以上離れている場合に判定を行う
519 if( len0 > lth )
520 {
521 vector_type p1 = ite1->pos;
522 vector_type p2 = ( ite1 + 1 )->pos;
523 double length = ( p1 - p2 ).length( );
524 vector_type dir = ( p1 - p2 ).unit( );
525 vector_type p = p1 + dir * length * 0.8;
526
527 double min = 1.0e10;
528 for( mesh_iterator cite = mesh_list.begin( ) ; cite != mesh_list.end( ) ; ++cite )
529 {
530 double l = ( cite->pos - p ).length( );
531 if( l < min && cite != ite1 && l < length )
532 {
533 min = l;
534 ratio1 = cite->round / ite1->round;
535 }
536 }
537 }
538
539 // 端から一定距離以上離れている場合に判定を行う
540 if( len1 > lth )
541 {
542 vector_type p1 = ite2->pos;
543 vector_type p2 = ( ite2 - 1 )->pos;
544 double length = ( p1 - p2 ).length( );
545 vector_type dir = ( p1 - p2 ).unit( );
546 vector_type p = p1 + dir * length * 0.8;
547
548 double min = 1.0e10;
549 for( mesh_iterator cite = mesh_list.begin( ) ; cite != mesh_list.end( ) ; ++cite )
550 {
551 double l = ( cite->pos - p ).length( );
552 if( l < min && cite != ite2 && l < length )
553 {
554 min = l;
555 ratio2 = cite->round / ite2->round;
556 }
557 }
558 }
559
560 if( ratio1 > ratio2 )
561 {
562 // 並び順を逆転させる
563 std::reverse( mesh_list.begin( ), mesh_list.begin( ) + base_mesh_num );
564 site = mesh_list.begin( );
565 eite = site + base_mesh_num;
566 }
567 else if( ratio1 == ratio2 && len0 < len1 )
568 {
569 // 並び順を逆転させる
570 std::reverse( mesh_list.begin( ), mesh_list.begin( ) + base_mesh_num );
571 site = mesh_list.begin( );
572 eite = site + base_mesh_num;
573 }
574 }
575
576 // 順番に値を設定する
577 difference_type mark = 1;
578 difference_type column = col;
579 for( ite = site ; ite != eite ; ++ite )
580 {
581 ite->mark = mark++;
582 ite->row = row;
583 ite->col = column--;
584 //std::cout << "( " << ite->row << ", " << ite->col << " )" << std::endl;
585 grid( ite->row, ite->col ) = ite->pos;
586 }
587 for( ; ite != mesh_list.end( ) ; ++ite )
588 {
589 ite->mark = 100;
590 }
591
592 // 出力データを整形する
593 for( ite = mesh_list.begin( ) ; ite != mesh_list.end( ) ; ++ite )
594 {
595 meshes[ ite->label ].mark = ite->mark;
596 }
597
598 // 基準点の数が少ない場合には終了
599 if( base_mesh_num < 3 )
600 {
601 for( ite = mesh_list.begin( ) ; ite != mesh_list.end( ) ; ++ite )
602 {
603 meshes[ ite->label ].mark = ite->mark;
604 }
605
606 for( i = 0 ; i < chart.size( ) ; i++ )
607 {
608 if( meshes[ binary[ i ] ].mark >= 0 && meshes[ binary[ i ] ].enabled )
609 {
610 chart[ i ] = static_cast< unsigned char >( meshes[ binary[ i ] ].mark );
611 }
612 else
613 {
614 chart[ i ] = 255;
615 }
616 }
617
618 return( false );
619 }
620
621
622 // 求まった大きな３つの円の配置を元にグリッドの座標を更新する
623 // その際，大きいものから３点のみを最初に確定させる
624 __mesh_information__ m1 = mesh_list[ 0 ];
625 __mesh_information__ m2 = mesh_list[ 2 ];
626 site = mesh_list.begin( );
627 eite = site + 3;
628
629 vector_type p1 = m1.pos;
630 vector_type p2 = m2.pos;
631 vector_type dir = ( p1 - p2 ).unit( );
632 vector_type norm( -dir.y, dir.x );
633 double length = ( p1 - p2 ).length( ) / static_cast< double >( 2 );
634
635 // 不用なデータを削除する
636 mesh_list.erase( site, eite );
637
638 // 基準ラインで確定しているものの近傍のみを確定させる
639 difference_type rindex[] = { -1, 1 };
640 for( r = 0 ; r < 2 ; r++ )
641 {
642 for( c = m1.col ; c >= m2.col ; c-- )
643 {
644 vector_type &op = grid( m1.row + rindex[ r ], c );
645 vector_type p = grid( m1.row, c ) + rindex[ r ] * norm * length * 0.5;
646
647 double min = 1.0e10;
648 mesh_iterator cur = mesh_list.begin( );
649 for( mesh_iterator cite = mesh_list.begin( ) ; cite != mesh_list.end( ) ; ++cite )
650 {
651 double l = ( cite->pos - p ).length( );
652 if( l < min )
653 {
654 min = l;
655 cur = cite;
656 }
657 }
658
659 // 一定距離以上はなれた探索点は無視する
660 if( cur != mesh_list.end( ) && min < length && min < cur->length )
661 {
662 // 以前に割り当てられていた場合は切り替える
663 if( cur->row >= 0 || cur->col >= 0 )
664 {
665 grid( cur->row, cur->col ).x = -1;
666 grid( cur->row, cur->col ).y = -1;
667 }
668
669 cur->row = r;
670 cur->col = c;
671 cur->length = min;
672
673 // 対応点が見つかったので，以降の探索から除外候補に追加する
674 op = cur->pos;
675 cur->mark = -2;
676 }
677 }
678 }
679
680 for( mesh_iterator ite = mesh_list.begin( ) ; ite != mesh_list.end( ) ; )
681 {
682 if( ite->mark == -2 )
683 {
684 ite = mesh_list.erase( ite );
685 }
686 else
687 {
688 ++ite;
689 }
690 }
691 }
692
693 {
694 for( i = 0 ; i < chart.size( ) ; i++ )
695 {
696 if( meshes[ binary[ i ] ].mark >= 0 && meshes[ binary[ i ] ].enabled )
697 {
698 chart[ i ] = static_cast< unsigned char >( meshes[ binary[ i ] ].mark );
699 }
700 else
701 {
702 chart[ i ] = 255;
703 }
704 }
705 }
706
707
708 // 反復的に対応点を見つける
709 while( true )
710 {
711 difference_type ncount = 0;
712 for( r = 0 ; r < rows && ncount < rows * cols ; r++ )
713 {
714 for( c = 0 ; c < cols && ncount < rows * cols ; c++ )
715 {
716 vector_type p = grid( r, c );
717
718 if( p.x != -1 || p.y != -1 )
719 {
720 ncount++;
721 continue;
722 }
723
724 double search_length = 0.0;
725
726 // 近傍の状態を使って，グリッド上の点が存在すると思われる位置を予測する
727 if( 0 < c && c < cols - 1 && grid( r, c - 1 ).x >= 0 && grid( r, c + 1 ).x >= 0 )
728 {
729 p = ( grid( r, c - 1 ) + grid( r, c + 1 ) ) / 2.0;
730 search_length = ( grid( r, c - 1 ) - grid( r, c + 1 ) ).length( ) / 2.0;
731 }
732 else if( 0 < r && r < rows - 1 && grid( r - 1, c ).x >= 0 && grid( r + 1, c ).x >= 0 )
733 {
734 p = ( grid( r - 1, c ) + grid( r + 1, c ) ) / 2.0;
735 search_length = ( grid( r - 1, c ) - grid( r + 1, c ) ).length( ) / 2.0;
736 }
737 else if( c < cols - 3 && grid( r, c + 1 ).x >= 0 && grid( r, c + 2 ).x >= 0 && grid( r, c + 3 ).x >= 0 )
738 {
739 double l1 = ( grid( r, c + 2 ) - grid( r, c + 3 ) ).length( );
740 double l2 = ( grid( r, c + 1 ) - grid( r, c + 2 ) ).length( );
741 p = grid( r, c + 1 ) + ( grid( r, c + 1 ) - grid( r, c + 2 ) ).unit( ) * ( 2.0 * l2 - l1 );
742 search_length = ( 2.0 * l2 - l1 ) / 2.0;
743 }
744 else if( c > 3 && grid( r, c - 1 ).x >= 0 && grid( r, c - 2 ).x >= 0 && grid( r, c - 3 ).x >= 0 )
745 {
746 double l1 = ( grid( r, c - 2 ) - grid( r, c - 3 ) ).length( );
747 double l2 = ( grid( r, c - 1 ) - grid( r, c - 2 ) ).length( );
748 p = grid( r, c - 1 ) + ( grid( r, c - 1 ) - grid( r, c - 2 ) ).unit( ) * ( 2.0 * l2 - l1 );
749 search_length = ( 2.0 * l2 - l1 ) / 2.0;
750 }
751 else if( r < rows - 3 && grid( r + 1, c ).x >= 0 && grid( r + 2, c ).x >= 0 && grid( r + 3, c ).x >= 0 )
752 {
753 double l1 = ( grid( r + 2, c ) - grid( r + 3, c ) ).length( );
754 double l2 = ( grid( r + 1, c ) - grid( r + 2, c ) ).length( );
755 p = grid( r + 1, c ) + ( grid( r + 1, c ) - grid( r + 2, c ) ).unit( ) * ( 2.0 * l2 - l1 );
756 search_length = ( 2.0 * l2 - l1 ) / 2.0;
757 }
758 else if( r > 3 && grid( r - 1, c ).x >= 0 && grid( r - 2, c ).x >= 0 && grid( r - 3, c ).x >= 0 )
759 {
760 double l1 = ( grid( r - 2, c ) - grid( r - 3, c ) ).length( );
761 double l2 = ( grid( r - 1, c ) - grid( r - 2, c ) ).length( );
762 p = grid( r - 1, c ) + ( grid( r - 1, c ) - grid( r - 2, c ) ).unit( ) * ( 2.0 * l2 - l1 );
763 search_length = ( 2.0 * l2 - l1 ) / 2.0;
764 }
765 else if( c > 1 && r < rows - 1 && grid( r + 1, c ).x >= 0 && grid( r, c - 1 ).x >= 0 && grid( r + 1, c - 1 ).x >= 0 )
766 {
767 vector_type &p0 = grid( r + 1, c - 1 );
768 vector_type &p1 = grid( r + 1, c );
769 vector_type &p2 = grid( r, c - 1 );
770 vector_type d = p1 + p2 - 2.0 * p0;
771
772 p = d + p0;
773 search_length = 0.5 * d.length( );
774 }
775 else if( c > 1 && r > 1 && grid( r - 1, c ).x >= 0 && grid( r, c - 1 ).x >= 0 && grid( r - 1, c - 1 ).x >= 0 )
776 {
777 vector_type &p0 = grid( r - 1, c - 1 );
778 vector_type &p1 = grid( r - 1, c );
779 vector_type &p2 = grid( r, c - 1 );
780 vector_type d = p1 + p2 - 2.0 * p0;
781
782 p = d + p0;
783 search_length = 0.5 * d.length( );
784 }
785 else if( c < cols - 1 && r > 1 && grid( r, c + 1 ).x >= 0 && grid( r - 1, c ).x >= 0 && grid( r - 1, c + 1 ).x >= 0 )
786 {
787 vector_type &p0 = grid( r - 1, c + 1 );
788 vector_type &p1 = grid( r - 1, c );
789 vector_type &p2 = grid( r, c + 1 );
790 vector_type d = p1 + p2 - 2.0 * p0;
791
792 p = d + p0;
793 search_length = 0.5 * d.length( );
794 }
795 else if( c < cols - 1 && r < rows - 1 && grid( r, c + 1 ).x >= 0 && grid( r + 1, c ).x >= 0 && grid( r + 1, c + 1 ).x >= 0 )
796 {
797 vector_type &p0 = grid( r + 1, c + 1 );
798 vector_type &p1 = grid( r + 1, c );
799 vector_type &p2 = grid( r, c + 1 );
800 vector_type d = p1 + p2 - 2.0 * p0;
801
802 p = d + p0;
803 search_length = 0.5 * d.length( );
804 }
805 else
806 {
807 // 近傍の点が見つからなかった
808 ncount++;
809 continue;
810 }
811
812 // 予測された点を用いて最寄の点を探す
813 double min = 1.0e10;
814 mesh_iterator cur = mesh_list.begin( );
815 for( mesh_iterator cite = mesh_list.begin( ) ; cite != mesh_list.end( ) ; ++cite )
816 {
817 double l = ( cite->pos - p ).length( );
818 if( l < min )
819 {
820 min = l;
821 cur = cite;
822 }
823 }
824
825 if( cur != mesh_list.end( ) && min < search_length && min < cur->length )
826 {
827 // 以前に割り当てられていた場合は切り替える
828 if( cur->row >= 0 || cur->col >= 0 )
829 {
830 grid( cur->row, cur->col ).x = -1;
831 grid( cur->row, cur->col ).y = -1;
832 }
833
834 cur->row = r;
835 cur->col = c;
836 cur->length = min;
837
838 // 対応点が見つかったので，以降の探索から除外候補に追加する
839 grid( r, c ) = cur->pos;
840 cur->mark = -2;
841 }
842 else
843 {
844 // 対応点が見つからなかったので，見つからなかったマークを入れる
845 grid( r, c ).x = -2;
846 grid( r, c ).y = -2;
847 }
848 }
849 }
850
851 for( mesh_iterator ite = mesh_list.begin( ) ; ite != mesh_list.end( ) ; )
852 {
853 if( ite->mark == -2 )
854 {
855 ite = mesh_list.erase( ite );
856 }
857 else
858 {
859 ++ite;
860 }
861 }
862
863 if( ncount == rows * cols )
864 {
865 // 更新すべき点が見つからなかったので終了する
866 break;
867 }
868 }
869
870 // 孤立点を削除する
871 for( r = 0 ; r < rows ; r++ )
872 {
873 for( c = 0 ; c < cols ; c++ )
874 {
875 if( grid( r, c ).x < 0 )
876 {
877 continue;
878 }
879
880 bool b1 = r > 0 ? grid( r - 1, c ).x >= 0 : false;
881 bool b2 = r < rows ? grid( r + 1, c ).x >= 0 : false;
882 bool b3 = c > 0 ? grid( r, c - 1 ).x >= 0 : false;
883 bool b4 = c < cols ? grid( r, c + 1 ).x >= 0 : false;
884
885 size_t count = 0;
886 if( r > 0 && grid( r - 1, c ).x >= 0 )
887 {
888 count++;
889 }
890 if( r < rows && grid( r + 1, c ).x >= 0 )
891 {
892 count++;
893 }
894 if( c > 0 && grid( r, c - 1 ).x >= 0 )
895 {
896 count++;
897 }
898 if( c < cols && grid( r, c + 1 ).x >= 0 )
899 {
900 count++;
901 }
902
903 if( count <= 1 )
904 {
905 grid( r, c ).x = -2;
906 grid( r, c ).y = -2;
907 }
908 }
909 }
910
911 return( true );
912 }
913
914
916 // メッシュ抽出グループの終わり
917
918
919 // mist名前空間の終わり
920 _MIST_END
921
922
923 #endif // __INCLUDE_MIST_TIMER__
924