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

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

facet.h
1 //
2 // Copyright (c) 2003-2011, MIST Project, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28
33
34 #ifndef __INCLUDE_MIST_FACET__
35 #define __INCLUDE_MIST_FACET__
36
37
38
39 #ifndef __INCLUDE_MIST_CONF_H__
40 #include "config/mist_conf.h"
41 #endif
42
43 #ifndef __INCLUDE_MIST_TYPE_TRAIT_H__
44 #include "config/type_trait.h"
45 #endif
46
47 #ifndef __INCLUDE_MIST_MATRIX__
48 #include "matrix.h"
49 #endif
50
51 #ifndef __INCLUDE_MIST_VECTOR__
52 #include "vector.h"
53 #endif
54
55
56 #include <vector>
57 #include <map>
58 #include <set>
59 #include <deque>
60 #include <algorithm>
61
62 // mist名前空間の始まり
63 _MIST_BEGIN
64
65
66 #define __SHOW_FACET_DEBUG_INFORMATION__ 0
67
75
80 template < class T >
81 class facet
82 {
83 public:
84 typedef T value_type;
85 typedef size_t size_type;
86 typedef ptrdiff_t difference_type;
87 typedef vector3< T > vector_type;
88 typedef typename float_type< T >::value_type float_type;
89
90 public:
91 vector_type normal;
92 vector_type p1;
93 vector_type p2;
94 vector_type p3;
95
97 facet( ) : normal( 0, 0, 1 ){ }
98
100 facet( const vector_type &N, const vector_type &P1, const vector_type &P2, const vector_type &P3 ) : normal( N ), p1( P1 ), p2( P2 ), p3( P3 ){ }
101
103 facet( const vector_type &P1, const vector_type &P2, const vector_type &P3 ) : normal( ( ( P2 - P1 ).outer( P3 - P1 ) ).unit( ) ), p1( P1 ), p2( P2 ), p3( P3 ){ }
104
106 template < class TT >
107 facet( const facet< TT > &f ) : normal( f.normal ), p1( f.p1 ), p2( f.p2 ), p3( f.p3 ){ }
108
109
111 template < class TT >
112 const facet &operator =( const facet< TT > &f )
113 {
114 if( &f != this )
115 {
116 normal = f.normal;
117 p1 = f.p1;
118 p2 = f.p2;
119 p3 = f.p3;
120 }
121 return ( *this );
122 }
123 };
124
125
137 template < class T > inline std::ostream &operator <<( std::ostream &out, const facet< T > &f )
138 {
139 out << "(" << f.normal << ") - ";
140 out << "(" << f.p1 << "), ";
141 out << "(" << f.p2 << "), ";
142 out << "(" << f.p3 << ")";
143 return( out );
144 }
145
146
151 template < class T >
152 class facet_list : public std::vector< facet< T > >
153 {
154 private:
155 typedef std::vector< facet< T > > base;
156
157 public:
158 typedef facet< T > facet_type;
159 typedef size_t size_type;
160 typedef ptrdiff_t difference_type;
161 typedef typename facet_type::vector_type vector_type;
162 typedef typename float_type< T >::value_type float_type;
163
164 public:
165 std::string name;
166
168 facet_list( ){ }
169
171 facet_list( const std::string &name_ ) : name( name_ ){ }
172
174 template < class TT >
175 facet_list( const facet_list< TT > &f ) : base( f ), name( f.name ){ }
176
177
179 template < class TT >
180 const facet_list &operator =( const facet_list< TT > &f )
181 {
182 if( &f != this )
183 {
184 base::operator =( f );
185 name = f.name;
186 }
187 return ( *this );
188 }
189 };
190
191
203 template < class T, class T1, class T2 >
204 inline bool convert_to_vertex_face_list( const facet_list< T > &facets, std::vector< vector3< T1 > > &vertices, std::vector< vector3< T2 > > &faces, const double eps = 1.0e-6 )
205 {
206 if( facets.empty( ) )
207 {
208 return( false );
209 }
210
211 typedef typename facet_list< T >::size_type size_type;
212 typedef typename facet_list< T >::difference_type difference_type;
213 typedef vector3< T1 > vector_type;
214 typedef vector3< T2 > ivector_type;
215
216 typedef std::pair< size_type, vector_type > value_type;
217 typedef std::multimap< difference_type, value_type > map_type;
218 typedef std::pair< difference_type, value_type > key_type;
219 map_type table;
220
221 faces.clear( );
222 faces.reserve( facets.size( ) );
223
224 double SCALE = 0.1 / eps;
225
226 // データをテーブルに登録する
227 for( size_type i = 0 ; i < facets.size( ) ; i++ )
228 {
229 const typename facet_list< T >::facet_type &f = facets[ i ];
230 size_type indx[ 3 ] = { 0, 0, 0 };
231 vector_type v[ 3 ] = { f.p1, f.p2, f.p3 };
232
233 for( size_type j = 0 ; j < 3 ; j++ )
234 {
235 const vector_type &vec = v[ j ];
236 difference_type k1 = static_cast< difference_type > ( vec.x * SCALE );
237 difference_type k2 = static_cast< difference_type > ( vec.y * SCALE );
238 difference_type k3 = static_cast< difference_type > ( vec.z * SCALE );
239 difference_type key = k1 ^ k2 ^ k3;
240
241 typename map_type::iterator ite = table.find( key );
242 if( ite == table.end( ) )
243 {
244 indx[ j ] = table.size( );
245 table.insert( key_type( key, value_type( indx[ j ], vec ) ) );
246 }
247 else
248 {
249 typename map_type::iterator upper = table.upper_bound( key );
250 for( ; ite != upper ; ++ite )
251 {
252 if( ( ite->second.second - vec ).length( ) < eps )
253 {
254 break;
255 }
256 }
257
258 if( ite == upper )
259 {
260 indx[ j ] = table.size( );
261 table.insert( key_type( key, value_type( indx[ j ], vec ) ) );
262 }
263 else
264 {
265 indx[ j ] = ite->second.first;
266 }
267 }
268 }
269
270 // 頂点を反時計回りに並べる
271 if( f.normal.inner( ( v[ 1 ] - v[ 0 ] ).outer( v[ 2 ] - v[ 0 ] ) ) < 0 )
272 {
273 size_type tmp = indx[ 1 ];
274 indx[ 1 ] = indx[ 2 ];
275 indx[ 2 ] = tmp;
276 }
277
278 if( indx[ 0 ] != indx[ 1 ] && indx[ 1 ] != indx[ 2 ] && indx[ 2 ] != indx[ 0 ] )
279 {
280 faces.push_back( ivector_type( indx[ 0 ], indx[ 1 ], indx[ 2 ] ) );
281 }
282 }
283
284 vertices.clear( );
285 vertices.resize( table.size( ) );
286 for( typename map_type::iterator ite = table.begin( ) ; ite != table.end( ) ; ++ite )
287 {
288 vertices[ ite->second.first ] = ite->second.second;
289 }
290
291
292 return( true );
293 }
294
295
296 namespace __mc__
297 {
298 struct __face__
299 {
300 typedef size_t size_type;
301 typedef ptrdiff_t difference_type;
302
303 difference_type eid1;
304 difference_type eid2;
305 difference_type eid3;
306 bool flag;
307
308 __face__( ) : eid1( 0 ), eid2( 0 ), eid3( 0 ), flag( true ){ }
309 };
310
311 template < class VERTEX >
312 struct __edge__
313 {
314 typedef size_t size_type;
315 typedef ptrdiff_t difference_type;
316 typedef VERTEX vector_type;
317
318 difference_type v1;
319 difference_type v2;
320 difference_type fid1;
321 difference_type fid2;
322 vector_type v;
323 double err;
324
325 __edge__( ) : fid1( 0 ), fid2( 0 ){ }
326 __edge__( difference_type vv1, difference_type vv2 ) : v1( vv1 ), v2( vv2 ), fid1( 0 ), fid2( 0 ), err( 0.0 ) { }
327 __edge__( difference_type vv1, difference_type vv2, difference_type id1 ) : v1( vv1 ), v2( vv2 ), fid1( id1 ), fid2( 0 ), err( 0.0 ) { }
328
329 bool operator <( const __edge__ &v ) const
330 {
331 return( err < v.err );
332 }
333 };
334
335 template < class EDGE >
336 struct __edge_comp__
337 {
338 typedef EDGE edge_type;
339 typedef typename edge_type::difference_type difference_type;
340
341 const std::vector< edge_type > &edges;
342
343 __edge_comp__( const std::vector< edge_type > &edge_list ) : edges( edge_list )
344 {
345 }
346
347 bool operator ()( difference_type v1, difference_type v2 ) const
348 {
349 return( edges[ v1 ].err < edges[ v2 ].err );
350 }
351 };
352
353 template< class T1, class T2 >
354 inline T1 create_key( const T1 v1, const T1 v2, const T2 stride )
355 {
356 if( v1 < v2 )
357 {
358 return( v1 * stride + v2 );
359 }
360 else
361 {
362 return( v2 * stride + v1 );
363 }
364 }
365
366
367 template< class T >
368 inline T find_edge( const T V1, const T V2, const T V3, const T v1, const T v2 )
369 {
370 if( V1 == v1 )
371 {
372 if( V2 == v2 )
373 {
374 return( 0 );
375 }
376 else if( V3 == v2 )
377 {
378 return( 2 );
379 }
380 }
381 else if( V2 == v1 )
382 {
383 if( V3 == v2 )
384 {
385 return( 1 );
386 }
387 else if( V1 == v2 )
388 {
389 return( 0 );
390 }
391 }
392 else if( V3 == v1 )
393 {
394 if( V1 == v2 )
395 {
396 return( 2 );
397 }
398 else if( V2 == v2 )
399 {
400 return( 1 );
401 }
402 }
403
404 return( -1 );
405 }
406
407 template < class T1, class T2, class T3 >
408 inline bool convert_to_vertex_face_list( const std::vector< vector3< T1 > > &vertices, const std::vector< vector3< T2 > > &faces, std::vector< __face__ > &face_lists, std::vector< __edge__< T3 > > &edge_lists )
409 {
410 if( vertices.empty( ) || faces.empty( ) )
411 {
412 return( false );
413 }
414
415 typedef typename vector3< T1 >::size_type size_type;
416 typedef typename vector3< T1 >::difference_type difference_type;
417 typedef vector3< T1 > vector_type;
418 typedef vector3< T2 > ivector_type;
419 typedef __face__ face_type;
420 typedef __edge__< T3 > edge_type;
421
422 typedef std::map< difference_type, difference_type > map_type;
423 typedef typename map_type::iterator iterator;
424 map_type table;
425
426 edge_lists.clear( );
427 face_lists.clear( );
428 face_lists.reserve( faces.size( ) );
429
430 // ダミーを1つだけ挿入する
431 face_lists.push_back( face_type( ) );
432 edge_lists.push_back( edge_type( ) );
433
434 // 辺をテーブルに登録する
435 for( size_type i = 1 ; i <= faces.size( ) ; i++ )
436 {
437 const ivector_type &f = faces[ i - 1 ];
438
439 difference_type key12 = create_key( f.x, f.y, vertices.size( ) );
440 difference_type key23 = create_key( f.y, f.z, vertices.size( ) );
441 difference_type key31 = create_key( f.z, f.x, vertices.size( ) );
442
443 // 辺1-2を調べる
444 iterator ite = table.find( key12 );
445 if( ite == table.end( ) )
446 {
447 // 新規に辺と面の対応情報を追加する
448 table[ key12 ] = edge_lists.size( );
449
450 // エッジを登録する
451 if( f.x < f.y )
452 {
453 edge_lists.push_back( edge_type( f.x, f.y, i ) );
454 }
455 else if( f.y < f.x )
456 {
457 edge_lists.push_back( edge_type( f.y, f.x, i ) );
458 }
459 else
460 {
461 std::cerr << "Incorrect edge is found." << std::endl;
462 }
463 }
464 else if( edge_lists[ ite->second ].fid2 == 0 )
465 {
466 edge_lists[ ite->second ].fid2 = i;
467 }
468 else
469 {
470 std::cerr << "Edge may be shared among more than three faces." << std::endl;
471 }
472
473 // 辺2-3を調べる
474 ite = table.find( key23 );
475 if( ite == table.end( ) )
476 {
477 // 新規に辺情報を追加する
478 table[ key23 ] = edge_lists.size( );
479
480 // エッジを登録する
481 if( f.y < f.z )
482 {
483 edge_lists.push_back( edge_type( f.y, f.z, i ) );
484 }
485 else if( f.z < f.y )
486 {
487 edge_lists.push_back( edge_type( f.z, f.y, i ) );
488 }
489 else
490 {
491 std::cerr << "Incorrect edge is found." << std::endl;
492 }
493 }
494 else if( edge_lists[ ite->second ].fid2 == 0 )
495 {
496 edge_lists[ ite->second ].fid2 = i;
497 }
498 else
499 {
500 std::cerr << "Edge may be shared among more than three faces." << std::endl;
501 }
502
503 // 辺2-3を調べる
504 ite = table.find( key31 );
505 if( ite == table.end( ) )
506 {
507 // 新規に辺情報を追加する
508 table[ key31 ] = edge_lists.size( );
509
510 // エッジを登録する
511 if( f.z < f.x )
512 {
513 edge_lists.push_back( edge_type( f.z, f.x, i ) );
514 }
515 else if( f.x < f.z )
516 {
517 edge_lists.push_back( edge_type( f.x, f.z, i ) );
518 }
519 else
520 {
521 std::cerr << "Incorrect edge is found." << std::endl;
522 }
523 }
524 else if( edge_lists[ ite->second ].fid2 == 0 )
525 {
526 edge_lists[ ite->second ].fid2 = i;
527 }
528 else
529 {
530 std::cerr << "Edge may be shared among more than three faces." << std::endl;
531 }
532 }
533
534 // 面と辺を関連付ける
535 for( size_type i = 1 ; i <= faces.size( ) ; i++ )
536 {
537 const ivector_type &F = faces[ i - 1 ];
538 face_type f;
539
540 difference_type key12 = create_key( F.x, F.y, vertices.size( ) );
541 difference_type key23 = create_key( F.y, F.z, vertices.size( ) );
542 difference_type key31 = create_key( F.z, F.x, vertices.size( ) );
543
544 // 辺1-2を調べる
545 iterator ite = table.find( key12 );
546 if( ite == table.end( ) )
547 {
548 std::cerr << "Can't find edge 1 from edge table." << std::endl;
549 }
550 else
551 {
552 edge_type &e = edge_lists[ ite->second ];
553 f.eid1 = F.x == e.v1 ? ite->second : -ite->second;
554 }
555
556 // 辺2-3を調べる
557 ite = table.find( key23 );
558 if( ite == table.end( ) )
559 {
560 std::cerr << "Can't find edge 2 from edge table." << std::endl;
561 }
562 else
563 {
564 edge_type &e = edge_lists[ ite->second ];
565 f.eid2 = F.y == e.v1 ? ite->second : -ite->second;
566 }
567
568 // 辺2-3を調べる
569 ite = table.find( key31 );
570 if( ite == table.end( ) )
571 {
572 std::cerr << "Can't find edge 3 from edge table." << std::endl;
573 }
574 else
575 {
576 edge_type &e = edge_lists[ ite->second ];
577 f.eid3 = F.z == e.v1 ? ite->second : -ite->second;
578 }
579
580 face_lists.push_back( f );
581 }
582
583 return( true );
584 }
585
586 template < class T >
587 inline T ABS( T val )
588 {
589 return( val < 0 ? -val : val );
590 }
591
592 template < class T >
593 void contract_edges( std::vector< __edge__< T > > &edges, std::vector< __face__ > &faces, __face__::difference_type eid0, __face__::difference_type eid1, __face__::difference_type EID, __face__::difference_type FID1, __face__::difference_type FID2 )
594 {
595 typedef __edge__< T > edge_type;
596 typedef __face__ face_type;
597 typedef __face__::difference_type difference_type;
598
599 difference_type fid2 = 0;
600 const edge_type &E = edges[ EID ];
601 edge_type &e0 = edges[ eid0 ];
602 edge_type &e1 = edges[ eid1 ];
603
604 if( e0.fid1 == FID1 )
605 {
606 if( e1.fid1 == FID2 )
607 {
608 fid2 = e1.fid2;
609 e0.fid1 = e1.fid2;
610 }
611 else if( e1.fid2 == FID2 )
612 {
613 fid2 = e1.fid1;
614 e0.fid1 = e1.fid1;
615 }
616 else
617 {
618 std::cerr << "Face-edge correspondence is incorrect!!" << std::endl;
619 }
620
621 }
622 else if( e0.fid2 == FID1 )
623 {
624 if( e1.fid1 == FID2 )
625 {
626 fid2 = e1.fid2;
627 e0.fid2 = e1.fid2;
628 }
629 else if( e1.fid2 == FID2 )
630 {
631 fid2 = e1.fid1;
632 e0.fid2 = e1.fid1;
633 }
634 else
635 {
636 std::cerr << "Face-edge correspondence is incorrect!!" << std::endl;
637 }
638 }
639 else
640 {
641 std::cerr << "Unknown error occured!!" << std::endl;
642 }
643
644 if( fid2 != 0 )
645 {
646 face_type &f2 = faces[ fid2 ];
647
648 if( ABS( f2.eid1 ) == eid1 )
649 {
650 if( e0.v1 == e1.v1 || e0.v2 == e1.v2 )
651 {
652 f2.eid1 = f2.eid1 < 0 ? -eid0 : eid0;
653 }
654 else if( e0.v1 == e1.v2 || e0.v2 == e1.v1 )
655 {
656 f2.eid1 = f2.eid1 < 0 ? eid0 : -eid0;
657 }
658 else
659 {
660 std::cerr << "Vertex-edge map is incorrect!! " << "(" << e1.v1 << ", " << e1.v2 << ") <-> (" << E.v1 << ", " << E.v2 << ")" << std::endl;
661 }
662 }
663 else if( ABS( f2.eid2 ) == eid1 )
664 {
665 if( e0.v1 == e1.v1 || e0.v2 == e1.v2 )
666 {
667 f2.eid2 = f2.eid2 < 0 ? -eid0 : eid0;
668 }
669 else if( e0.v1 == e1.v2 || e0.v2 == e1.v1 )
670 {
671 f2.eid2 = f2.eid2 < 0 ? eid0 : -eid0;
672 }
673 else
674 {
675 std::cerr << "Vertex-edge map is incorrect!! " << "(" << e1.v1 << ", " << e1.v2 << ") <-> (" << E.v1 << ", " << E.v2 << ")" << std::endl;
676 }
677 }
678 else if( ABS( f2.eid3 ) == eid1 )
679 {
680 if( e0.v1 == e1.v1 || e0.v2 == e1.v2 )
681 {
682 f2.eid3 = f2.eid3 < 0 ? -eid0 : eid0;
683 }
684 else if( e0.v1 == e1.v2 || e0.v2 == e1.v1 )
685 {
686 f2.eid3 = f2.eid3 < 0 ? eid0 : -eid0;
687 }
688 else
689 {
690 std::cerr << "Vertex-edge map is incorrect!! " << "(" << e1.v1 << ", " << e1.v2 << ") <-> (" << E.v1 << ", " << E.v2 << ")" << std::endl;
691 }
692 }
693 else
694 {
695 std::cerr << "Face-edge correspondence is incorrect!!" << std::endl;
696 }
697 }
698 }
699
700 template < class T, class Allocator >
701 inline bool inv3x3( matrix< T, Allocator > &a, const double eps = 1.0e-6 )
702 {
703 typedef typename matrix< T, Allocator >::value_type value_type;
704
705 double a11 = a( 0, 0 ), a12 = a( 0, 1 ), a13 = a( 0, 2 );
706 double a21 = a( 1, 0 ), a22 = a( 1, 1 ), a23 = a( 1, 2 );
707 double a31 = a( 2, 0 ), a32 = a( 2, 1 ), a33 = a( 2, 2 );
708
709 // 共通で利用する係数の計算
710 double _22x33_23x32_ = a22 * a33 - a23 * a32;
711 double _21x32_22x31_ = a21 * a32 - a22 * a31;
712 double _23x31_21x33_ = a23 * a31 - a21 * a33;
713
714 // 行列式を計算
715 double detA = a11 * _22x33_23x32_ + a13 * _21x32_22x31_ + a12 * _23x31_21x33_;
716
717 // 逆行列が存在する場合のも逆行列を計算する
718 if( std::abs( detA ) > eps )
719 {
720 // 各要素の値を計算
721 double A11 = _22x33_23x32_;
722 double A12 = a13 * a32 - a12 * a33;
723 double A13 = a12 * a23 - a13 * a22;
724 double A21 = _23x31_21x33_;
725 double A22 = a11 * a33 - a13 * a31;
726 double A23 = a13 * a21 - a11 * a23;
727 double A31 = _21x32_22x31_;
728 double A32 = a12 * a31 - a11 * a32;
729 double A33 = a11 * a22 - a12 * a21;
730
731 double _1_detA = 1.0 / detA;
732 a( 0, 0 ) = static_cast< value_type >( A11 * _1_detA );
733 a( 0, 1 ) = static_cast< value_type >( A12 * _1_detA );
734 a( 0, 2 ) = static_cast< value_type >( A13 * _1_detA );
735 a( 1, 0 ) = static_cast< value_type >( A21 * _1_detA );
736 a( 1, 1 ) = static_cast< value_type >( A22 * _1_detA );
737 a( 1, 2 ) = static_cast< value_type >( A23 * _1_detA );
738 a( 2, 0 ) = static_cast< value_type >( A31 * _1_detA );
739 a( 2, 1 ) = static_cast< value_type >( A32 * _1_detA );
740 a( 2, 2 ) = static_cast< value_type >( A33 * _1_detA );
741
742 return( true );
743 }
744 else
745 {
746 return( false );
747 }
748 }
749
750 template < class Matrix, class Vector >
751 inline double compute_vertex_error( const Matrix &Q, const Vector &v )
752 {
753 double e0 = Q( 0, 0 ) * v.x + Q( 0, 1 ) * v.y + Q( 0, 2 ) * v.z + Q( 0, 3 );
754 double e1 = Q( 1, 0 ) * v.x + Q( 1, 1 ) * v.y + Q( 1, 2 ) * v.z + Q( 1, 3 );
755 double e2 = Q( 2, 0 ) * v.x + Q( 2, 1 ) * v.y + Q( 2, 2 ) * v.z + Q( 2, 3 );
756 double e3 = Q( 3, 0 ) * v.x + Q( 3, 1 ) * v.y + Q( 3, 2 ) * v.z + Q( 3, 3 );
757 return( e0 * v.x + e1 * v.y + e2 * v.z + e3 );
758 }
759
760 template < class T1, class T2, class Matrix >
761 inline void update_edge( const std::vector< vector3< T1 > > &vertices, const std::vector< Matrix > &Q, __edge__< T2 > &edge, bool use_optimal_vertex_placement )
762 {
763 typedef Matrix matrix_type;
764 typedef vector3< T1 > vector_type;
765
766 // 各頂点の誤差評価行列を求める
767 matrix_type QQ = Q[ edge.v1 ] + Q[ edge.v2 ];
768 matrix_type QQQ = matrix_type::_33( QQ( 0, 0 ), QQ( 0, 1 ), QQ( 0, 2 ),
769 QQ( 1, 0 ), QQ( 1, 1 ), QQ( 1, 2 ),
770 QQ( 2, 0 ), QQ( 2, 1 ), QQ( 2, 2 ) );
771
772
773 if( use_optimal_vertex_placement && inv3x3( QQQ ) )
774 {
775 matrix_type V = QQQ * matrix_type::_31( QQ( 0, 3 ), QQ( 1, 3 ), QQ( 2, 3 ) );
776 edge.v = vector_type( -V[ 0 ], -V[ 1 ], -V[ 2 ] );
777 edge.err = compute_vertex_error( QQ, edge.v );
778 }
779 else
780 {
781 vector_type vs = vertices[ edge.v1 ];
782 vector_type ve = vertices[ edge.v2 ];
783 edge.err = 1.0e300;
784
785 for( double s = 0.0 ; s <= 1.0 ; s += 0.0625 )
786 {
787 vector_type v = vs + ( ve - vs ) * s;
788 double err = compute_vertex_error( QQ, v );
789
790 if( err < edge.err )
791 {
792 edge.v = v;
793 edge.err = err;
794 }
795 }
796 }
797 }
798
799 template < class T1, class T2 >
800 inline double compute_penalty( const std::vector< vector3< T1 > > &vertices, const std::vector< __edge__< T2 > > &edges, const __face__ &f, typename __edge__< T2 >::difference_type vid, const typename __edge__< T2 >::vector_type &v )
801 {
802 typedef __face__ face_type;
803 typedef __edge__< T2 > edge_type;
804 typedef typename edge_type::vector_type vector_type;
805 typedef typename edge_type::difference_type difference_type;
806
807 const vector_type &v1 = v;
808 vector_type v2, v3;
809
810 const edge_type &e1 = edges[ f.eid1 < 0 ? -f.eid1 : f.eid1 ];
811 const edge_type &e2 = edges[ f.eid2 < 0 ? -f.eid2 : f.eid2 ];
812 const edge_type &e3 = edges[ f.eid3 < 0 ? -f.eid3 : f.eid3 ];
813 difference_type vid1 = f.eid1 < 0 ? e1.v2 : e1.v1;
814 difference_type vid2 = f.eid2 < 0 ? e2.v2 : e2.v1;
815 difference_type vid3 = f.eid3 < 0 ? e3.v2 : e3.v1;
816
817 if( vid1 == vid )
818 {
819 v2 = vertices[ vid2 ];
820 v3 = vertices[ vid3 ];
821 }
822 else if( vid2 == vid )
823 {
824 v2 = vertices[ vid1 ];
825 v3 = vertices[ vid3 ];
826 }
827 else if( vid3 == vid )
828 {
829 v2 = vertices[ vid1 ];
830 v3 = vertices[ vid2 ];
831 }
832 else
833 {
834 std::cout << "Can't find vertices from the specified face." << std::endl;
835 }
836
837 const double coeff = 6.9282032302755091741097853660235;
838 double w = 0.5 * ( ( v2 - v1 ).outer( v3 - v1 ) ).length( );
839 double l1 = ( v2 - v1 ).inner( v2 - v1 );
840 double l2 = ( v3 - v2 ).inner( v3 - v2 );
841 double l3 = ( v1 - v3 ).inner( v1 - v3 );
842
843 return( coeff * w / ( l1 + l2 + l3 ) );
844 }
845
846 template < class T1, class T2 >
847 inline bool is_collapsed_after_contraction( const std::vector< vector3< T1 > > &vertices, const std::vector< __edge__< T2 > > &edges, const __face__ &f, typename __edge__< T2 >::difference_type vid, const typename __edge__< T2 >::vector_type &v )
848 {
849 typedef __face__ face_type;
850 typedef __edge__< T2 > edge_type;
851 typedef typename edge_type::vector_type vector_type;
852 typedef typename edge_type::difference_type difference_type;
853
854 const vector_type &v0 = vertices[ vid ];
855 const vector_type &v1 = v;
856 vector_type v2, v3;
857
858 const edge_type &e1 = edges[ f.eid1 < 0 ? -f.eid1 : f.eid1 ];
859 const edge_type &e2 = edges[ f.eid2 < 0 ? -f.eid2 : f.eid2 ];
860 const edge_type &e3 = edges[ f.eid3 < 0 ? -f.eid3 : f.eid3 ];
861 difference_type vid1 = f.eid1 < 0 ? e1.v2 : e1.v1;
862 difference_type vid2 = f.eid2 < 0 ? e2.v2 : e2.v1;
863 difference_type vid3 = f.eid3 < 0 ? e3.v2 : e3.v1;
864
865 if( vid1 == vid )
866 {
867 v2 = vertices[ vid2 ];
868 v3 = vertices[ vid3 ];
869 }
870 else if( vid2 == vid )
871 {
872 v2 = vertices[ vid1 ];
873 v3 = vertices[ vid3 ];
874 }
875 else if( vid3 == vid )
876 {
877 v2 = vertices[ vid1 ];
878 v3 = vertices[ vid2 ];
879 }
880 else
881 {
882 std::cout << "Can't find vertices from the specified face." << std::endl;
883 }
884
885 vector_type n1 = ( v2 - v0 ).outer( v3 - v0 );
886 vector_type n2 = ( v2 - v1 ).outer( v3 - v1 );
887
888 return( n1.inner( n2 ) < 0.0 );
889 }
890
891 template < class MULTIMAP, class T1, class T2 >
892 inline void apply_penalties( const MULTIMAP &vertex_edge_map, const std::vector< vector3< T1 > > &vertices, const std::vector< __face__ > &faces, std::vector< __edge__< T2 > > &edges, typename __edge__< T2 >::difference_type EID, double threshold_for_triangle_compactness )
893 {
894 typedef MULTIMAP vertex_edge_map_type;
895 typedef typename vertex_edge_map_type::const_iterator const_iterator;
896
897 typedef __face__ face_type;
898 typedef __edge__< T2 > edge_type;
899 typedef typename edge_type::vector_type vector_type;
900 typedef typename edge_type::difference_type difference_type;
901
902 edge_type &edge = edges[ EID ];
903
904 double penalty = 1.0e100;
905
906 std::set< difference_type > face_map1;
907 std::set< difference_type > face_map2;
908
909 {
910 const_iterator ite = vertex_edge_map.find( edge.v1 );
911 if( ite != vertex_edge_map.end( ) )
912 {
913 const_iterator upper = vertex_edge_map.upper_bound( edge.v1 );
914 for( ; ite != upper ; ++ite )
915 {
916 face_map1.insert( edges[ ite->second ].fid1 );
917 face_map1.insert( edges[ ite->second ].fid2 );
918 }
919 }
920
921 for( typename std::set< difference_type >::iterator ite = face_map1.begin( ) ; ite != face_map1.end( ) ; ++ite )
922 {
923 if( *ite != 0 )
924 {
925 double ppp = compute_penalty( vertices, edges, faces[ *ite ], edge.v1, edge.v );
926 if( ppp < penalty )
927 {
928 penalty = ppp;
929 }
930 }
931 }
932 }
933
934 {
935 std::set< difference_type > face_map;
936 const_iterator ite = vertex_edge_map.find( edge.v2 );
937 if( ite != vertex_edge_map.end( ) )
938 {
939 const_iterator upper = vertex_edge_map.upper_bound( edge.v2 );
940 for( ; ite != upper ; ++ite )
941 {
942 face_map2.insert( edges[ ite->second ].fid1 );
943 face_map2.insert( edges[ ite->second ].fid2 );
944 }
945 }
946
947 for( typename std::set< difference_type >::iterator ite = face_map2.begin( ) ; ite != face_map2.end( ) ; ++ite )
948 {
949 if( *ite != 0 )
950 {
951 double ppp = compute_penalty( vertices, edges, faces[ *ite ], edge.v2, edge.v );
952 if( ppp < penalty )
953 {
954 penalty = ppp;
955 }
956 }
957 }
958 }
959
960 if( penalty < threshold_for_triangle_compactness )
961 {
962 edge.err += 1.0 - penalty;
963 }
964
965
966 // 法線が入れ替わるかどうかをチェックする
967 for( typename std::set< difference_type >::iterator ite = face_map1.begin( ) ; ite != face_map1.end( ) ; ++ite )
968 {
969 if( *ite != 0 && face_map2.find( *ite ) == face_map2.end( ) )
970 {
971 if( is_collapsed_after_contraction( vertices, edges, faces[ *ite ], edge.v1, edge.v ) )
972 {
973 edge.err += 1.0;
974 }
975 }
976 }
977
978 for( typename std::set< difference_type >::iterator ite = face_map2.begin( ) ; ite != face_map2.end( ) ; ++ite )
979 {
980 if( *ite != 0 && face_map1.find( *ite ) == face_map1.end( ) )
981 {
982 if( is_collapsed_after_contraction( vertices, edges, faces[ *ite ], edge.v2, edge.v ) )
983 {
984 edge.err += 1.0;
985 }
986 }
987 }
988 }
989
990 template < class T >
991 inline bool is_sharp_edge( const __edge__< T > &edge )
992 {
993 return( edge.fid1 * edge.fid2 == 0 );
994 }
995
996 template < class MULTIMAP, class T >
997 inline size_t number_of_sharp_edges( const MULTIMAP &vertex_edge_map, const std::vector< __edge__< T > > &edges, typename MULTIMAP::difference_type vid )
998 {
999 typedef MULTIMAP vertex_edge_map_type;
1000 typedef typename vertex_edge_map_type::const_iterator const_iterator;
1001 typedef typename vertex_edge_map_type::difference_type difference_type;
1002
1003 difference_type count = 0;
1004
1005 const_iterator ite = vertex_edge_map.find( vid );
1006 if( ite != vertex_edge_map.end( ) )
1007 {
1008 const_iterator upper = vertex_edge_map.upper_bound( vid );
1009 for( ; ite != upper ; ++ite )
1010 {
1011 if( is_sharp_edge( edges[ ite->second ] ) )
1012 {
1013 count++;
1014 }
1015 }
1016 }
1017
1018 return( count );
1019 }
1020
1021 template < class T >
1022 inline void compute_edge_connections( const std::vector< __edge__< T > > &edges, __face__::difference_type EID, __face__::difference_type eid,
1023 __face__::difference_type vs, __face__::difference_type vt, __face__::difference_type &vv,
1024 __face__::difference_type &eid1, __face__::difference_type &eid2 )
1025 {
1026 typedef __edge__< T > edge_type;
1027
1028 if( eid != EID )
1029 {
1030 const edge_type &e = edges[ eid ];
1031 if( e.v1 == vs && eid1 == 0 )
1032 {
1033 // vs へ接続する辺
1034 eid1 = eid;
1035 vv = e.v2;
1036 }
1037 else if( e.v2 == vs && eid1 == 0 )
1038 {
1039 // vs へ接続する辺
1040 eid1 = eid;
1041 vv = e.v1;
1042 }
1043 else if( e.v1 == vt )
1044 {
1045 // vt へ接続する辺
1046 eid2 = eid;
1047 vv = e.v2;
1048 }
1049 else if( e.v2 == vt )
1050 {
1051 // vt へ接続する辺
1052 eid2 = eid;
1053 vv = e.v1;
1054 }
1055 else
1056 {
1057 std::cerr << "Incorrect edge pairs." << std::endl;
1058 }
1059 }
1060 }
1061
1062 template < class T >
1063 inline bool compute_edge_connections( const std::vector< __face__ > &faces, const std::vector< __edge__< T > > &edges, __face__::difference_type EID,
1064 __face__::difference_type vs, __face__::difference_type vt, __face__::difference_type &vl, __face__::difference_type &vr, __face__::difference_type eid[ 4 ] )
1065 {
1066 typedef __edge__< T > edge_type;
1067 typedef __face__ face_type;
1068 typedef __face__::difference_type difference_type;
1069
1070 const edge_type &EDGE = edges[ EID ];
1071
1072 if( is_sharp_edge( EDGE ) )
1073 {
1074 // 縁に接している辺は削除対象外
1075 return( false );
1076 }
1077
1078 const face_type &f1 = faces[ EDGE.fid1 ];
1079 const face_type &f2 = faces[ EDGE.fid2 ];
1080
1081 eid[ 0 ] = eid[ 1 ] = eid[ 2 ] = eid[ 3 ] = 0;
1082
1083 const edge_type &e11 = edges[ f1.eid1 < 0 ? -f1.eid1 : f1.eid1 ];
1084 const edge_type &e12 = edges[ f1.eid2 < 0 ? -f1.eid2 : f1.eid2 ];
1085 const edge_type &e13 = edges[ f1.eid3 < 0 ? -f1.eid3 : f1.eid3 ];
1086 const edge_type &e21 = edges[ f2.eid1 < 0 ? -f2.eid1 : f2.eid1 ];
1087 const edge_type &e22 = edges[ f2.eid2 < 0 ? -f2.eid2 : f2.eid2 ];
1088 const edge_type &e23 = edges[ f2.eid3 < 0 ? -f2.eid3 : f2.eid3 ];
1089
1090 // 面1をチェック
1091 compute_edge_connections( edges, EID, f1.eid1 < 0 ? -f1.eid1 : f1.eid1, vs, vt, vl, eid[ 0 ], eid[ 1 ] );
1092 compute_edge_connections( edges, EID, f1.eid2 < 0 ? -f1.eid2 : f1.eid2, vs, vt, vl, eid[ 0 ], eid[ 1 ] );
1093 compute_edge_connections( edges, EID, f1.eid3 < 0 ? -f1.eid3 : f1.eid3, vs, vt, vl, eid[ 0 ], eid[ 1 ] );
1094
1095 // 面2をチェック
1096 compute_edge_connections( edges, EID, f2.eid1 < 0 ? -f2.eid1 : f2.eid1, vs, vt, vr, eid[ 3 ], eid[ 2 ] );
1097 compute_edge_connections( edges, EID, f2.eid2 < 0 ? -f2.eid2 : f2.eid2, vs, vt, vr, eid[ 3 ], eid[ 2 ] );
1098 compute_edge_connections( edges, EID, f2.eid3 < 0 ? -f2.eid3 : f2.eid3, vs, vt, vr, eid[ 3 ], eid[ 2 ] );
1099
1100 return( eid[ 0 ] * eid[ 1 ] * eid[ 2 ] * eid[ 3 ] != 0 );
1101 }
1102
1103 template < class MULTIMAP >
1104 inline void remove_edge_from_map( MULTIMAP &vertex_edge_map, typename MULTIMAP::key_type key, typename MULTIMAP::mapped_type val )
1105 {
1106 typename MULTIMAP::iterator ite = vertex_edge_map.find( key );
1107 if( ite != vertex_edge_map.end( ) )
1108 {
1109 typename MULTIMAP::iterator upper = vertex_edge_map.upper_bound( key );
1110 for( ; ite != upper ; )
1111 {
1112 if( ite->second == val )
1113 {
1114 // 削除する
1115 ite = vertex_edge_map.erase( ite );
1116 }
1117 else
1118 {
1119 ++ite;
1120 }
1121 }
1122 }
1123 }
1124
1125 template < class SET >
1126 inline void remove_edge_from_set( SET &edge_set, typename SET::value_type EID )
1127 {
1128 typename SET::iterator ite = edge_set.find( EID );
1129 if( ite != edge_set.end( ) )
1130 {
1131 edge_set.erase( ite );
1132 }
1133 }
1134
1135 template < class MULTIMAP, class T >
1136 inline bool check_topology_change( const MULTIMAP &vertex_edge_map, const std::vector< __face__ > &faces, const std::vector< __edge__< T > > &edges, __face__::difference_type EID )
1137 {
1138 typedef MULTIMAP vertex_edge_map_type;
1139 typedef typename vertex_edge_map_type::const_iterator const_iterator;
1140
1141 typedef __edge__< T > edge_type;
1142 typedef __face__ face_type;
1143 typedef __face__::difference_type difference_type;
1144
1145 const edge_type &EDGE = edges[ EID ];
1146
1147 if( is_sharp_edge( EDGE ) )
1148 {
1149 // 縁に接している辺は削除対象外
1150 return( true );
1151 }
1152
1153 // 処理対象の辺の登録と，テーブル内からの削除等を行う
1154 difference_type vs = EDGE.v1;
1155 difference_type vt = EDGE.v2;
1156 difference_type vl, vr;
1157 difference_type eid[ 4 ];
1158 if( !__mc__::compute_edge_connections( faces, edges, EID, vs, vt, vl, vr, eid ) )
1159 {
1160 return( true );
1161 }
1162
1163 if( vl == vr )
1164 {
1165 return( true );
1166 }
1167
1168 // 指定した辺を削除した場合にトポロジーが変化するかどうかをチェックする
1169 if( is_sharp_edge( edges[ eid[ 0 ] ] ) && is_sharp_edge( edges[ eid[ 1 ] ] ) )
1170 {
1171 return( true );
1172 }
1173 else if( is_sharp_edge( edges[ eid[ 3 ] ] ) && is_sharp_edge( edges[ eid[ 2 ] ] ) )
1174 {
1175 return( true );
1176 }
1177
1178 difference_type nvs = number_of_sharp_edges( vertex_edge_map, edges, vs );
1179 difference_type nvt = number_of_sharp_edges( vertex_edge_map, edges, vt );
1180
1181 if( nvs >= 1 && nvt >= 1 )
1182 {
1183 return( true );
1184 }
1185
1186 std::set< difference_type > vertex_list;
1187
1188 {
1189 const_iterator ite = vertex_edge_map.find( EDGE.v1 );
1190 if( ite != vertex_edge_map.end( ) )
1191 {
1192 const_iterator upper = vertex_edge_map.upper_bound( EDGE.v1 );
1193 for( ; ite != upper ; ++ite )
1194 {
1195 if( ite->second != EID )
1196 {
1197 const edge_type &e = edges[ ite->second ];
1198 if( e.fid1 != EDGE.fid1 && e.fid2 != EDGE.fid2 && e.fid1 != EDGE.fid2 && e.fid2 != EDGE.fid1 )
1199 {
1200 if( e.v1 != EDGE.v1 )
1201 {
1202 vertex_list.insert( e.v1 );
1203 }
1204 else
1205 {
1206 vertex_list.insert( e.v2 );
1207 }
1208 }
1209 }
1210 }
1211 }
1212 }
1213
1214 {
1215 const_iterator ite = vertex_edge_map.find( EDGE.v2 );
1216 if( ite != vertex_edge_map.end( ) )
1217 {
1218 const_iterator upper = vertex_edge_map.upper_bound( EDGE.v2 );
1219 for( ; ite != upper ; ++ite )
1220 {
1221 if( ite->second != EID )
1222 {
1223 const edge_type &e = edges[ ite->second ];
1224 if( e.fid1 != EDGE.fid1 && e.fid2 != EDGE.fid2 && e.fid1 != EDGE.fid2 && e.fid2 != EDGE.fid1 )
1225 {
1226 if( e.v1 != EDGE.v2 )
1227 {
1228 if( vertex_list.find( e.v1 ) != vertex_list.end( ) )
1229 {
1230 // 面の存在しないループを発見
1231 return( true );
1232 }
1233 }
1234 else
1235 {
1236 if( vertex_list.find( e.v2 ) != vertex_list.end( ) )
1237 {
1238 // 面の存在しないループを発見
1239 return( true );
1240 }
1241 }
1242 }
1243 }
1244 }
1245 }
1246 }
1247
1248 return( false );
1249 }
1250 }
1251
1268 template < class T >
1269 inline bool surface_simplification( facet_list< T > &facets, size_t number_of_facets, bool use_optimal_vertex_placement = true, double threshold_for_triangle_compactness = 0.0, const double eps = 1.0e-3 )
1270 {
1271 typedef typename facet_list< T >::facet_type facet_type;
1272 typedef typename facet_type::size_type size_type;
1273 typedef typename facet_type::difference_type difference_type;
1274 typedef typename facet_type::vector_type vector_type;
1275 typedef vector3< difference_type > ivector_type;
1276 typedef __mc__::__edge__< vector_type > edge_type;
1277 typedef __mc__::__face__ face_type;
1278 typedef matrix< double > matrix_type;
1279 typedef std::vector< edge_type > edge_list_type;
1280
1281
1282 std::vector< vector_type > vertices;
1283 std::vector< ivector_type > face_ids;
1284 std::vector< face_type > faces;
1285 edge_list_type edges;
1286
1287 // 頂点と面のリストに変換する
1288 if( !convert_to_vertex_face_list( facets, vertices, face_ids, eps ) )
1289 {
1290 return( false );
1291 }
1292
1293 // 面情報にエッジ情報を付与する
1294 if( !convert_to_vertex_face_list( vertices, face_ids, faces, edges ) )
1295 {
1296 return( false );
1297 }
1298 else
1299 {
1300 face_ids.clear( );
1301 }
1302
1303 // Contraction Error を計算するための行列集合を設定する
1304 std::vector< matrix_type > Q( vertices.size( ) );
1305 for( size_type i = 0 ; i < Q.size( ) ; i++ )
1306 {
1307 Q[ i ].resize( 4, 4 );
1308 }
1309
1310 // 接続する頂点の候補集合を求める
1311 for( size_type i = 1 ; i < faces.size( ) ; i++ )
1312 {
1313 const face_type &f = faces[ i ];
1314 const difference_type v1 = f.eid1 < 0 ? edges[ -f.eid1 ].v2 : edges[ f.eid1 ].v1;
1315 const difference_type v2 = f.eid2 < 0 ? edges[ -f.eid2 ].v2 : edges[ f.eid2 ].v1;
1316 const difference_type v3 = f.eid3 < 0 ? edges[ -f.eid3 ].v2 : edges[ f.eid3 ].v1;
1317 const vector_type &p1 = vertices[ v1 ];
1318 const vector_type &p2 = vertices[ v2 ];
1319 const vector_type &p3 = vertices[ v3 ];
1320 vector_type n = ( ( p2 - p1 ).outer( p3 - p1 ) ).unit( );
1321
1322 double a = n.x;
1323 double b = n.y;
1324 double c = n.z;
1325
1326 vector_type p0 = ( p1 + p2 + p3 ) / 3.0;
1327 vector_type d12 = ( p2 - p1 ).unit( );
1328 vector_type d23 = ( p3 - p2 ).unit( );
1329 vector_type d31 = ( p1 - p3 ).unit( );
1330 vector_type p12 = p1 + d12 * ( p0 - p1 ).inner( d12 );
1331 vector_type p23 = p2 + d23 * ( p0 - p2 ).inner( d23 );
1332 vector_type p31 = p3 + d31 * ( p0 - p3 ).inner( d31 );
1333
1334 {
1335 double area1 = ( ( p0 - p1 ).outer( p12 - p1 ) ).length( ) * 0.5;
1336 double area2 = ( ( p0 - p1 ).outer( p31 - p1 ) ).length( ) * 0.5;
1337 double d = -( a * p1.x + b * p1.y + c * p1.z );
1338 Q[ v1 ] += matrix_type::_44( a * a, a * b, a * c, a * d,
1339 a * b, b * b, b * c, b * d,
1340 a * c, b * c, c * c, c * d,
1341 a * d, b * d, c * d, d * d ) * ( area1 + area2 );
1342 }
1343
1344 {
1345 double area1 = ( ( p0 - p2 ).outer( p12 - p1 ) ).length( ) * 0.5;
1346 double area2 = ( ( p0 - p2 ).outer( p23 - p1 ) ).length( ) * 0.5;
1347 double d = -( a * p2.x + b * p2.y + c * p2.z );
1348 Q[ v2 ] += matrix_type::_44( a * a, a * b, a * c, a * d,
1349 a * b, b * b, b * c, b * d,
1350 a * c, b * c, c * c, c * d,
1351 a * d, b * d, c * d, d * d ) * ( area1 + area2 );
1352 }
1353
1354 {
1355 double area1 = ( ( p0 - p3 ).outer( p23 - p1 ) ).length( ) * 0.5;
1356 double area2 = ( ( p0 - p3 ).outer( p31 - p1 ) ).length( ) * 0.5;
1357 double d = -( a * p3.x + b * p3.y + c * p3.z );
1358 Q[ v3 ] += matrix_type::_44( a * a, a * b, a * c, a * d,
1359 a * b, b * b, b * c, b * d,
1360 a * c, b * c, c * c, c * d,
1361 a * d, b * d, c * d, d * d ) * ( area1 + area2 );
1362 }
1363
1364 }
1365
1366 // 頂点とエッジのテーブルを作成する
1367 typedef __mc__::__edge_comp__< edge_type > edge_compare_type;
1368 typedef std::set< difference_type, edge_compare_type > edge_map_type;
1369 typedef std::multimap< size_type, difference_type > vertex_edge_map_type;
1370 typedef std::pair< size_type, difference_type > vertex_edge_map_pair_type;
1371
1372 vertex_edge_map_type vertex_edge_map;
1373 edge_compare_type ecomp( edges );
1374 edge_map_type edge_map( ecomp );
1375 for( size_type i = 1 ; i < edges.size( ) ; i++ )
1376 {
1377 edge_type &e = edges[ i ];
1378
1379 // 削除対象にできるかどうかを判定
1380 if( !__mc__::is_sharp_edge( edges[ i ] ) )
1381 {
1382 // 各頂点の誤差評価を行う
1383 __mc__::update_edge( vertices, Q, edges[ i ], use_optimal_vertex_placement );
1384 edge_map.insert( i );
1385 }
1386
1387 vertex_edge_map.insert( vertex_edge_map_pair_type( e.v1, i ) );
1388 vertex_edge_map.insert( vertex_edge_map_pair_type( e.v2, i ) );
1389 }
1390
1391 size_t num_facets = faces.size( ) - 1;
1392
1393 // 頂点をコスト順に削減する
1394 for( ; num_facets - 2 >= number_of_facets && !edge_map.empty( ) ; num_facets -= 2 )
1395 {
1396 typename edge_map_type::iterator mite = edge_map.end( );
1397
1398 for( typename edge_map_type::iterator ite = edge_map.begin( ) ; ite != edge_map.end( ) ; )
1399 {
1400 const edge_type &e = edges[ *ite ];
1401
1402 if( __mc__::is_sharp_edge( e ) )
1403 {
1404 ite = edge_map.erase( ite );
1405 }
1406 else if( !__mc__::check_topology_change( vertex_edge_map, faces, edges, *ite ) )
1407 {
1408 // 削除可能な辺を発見
1409 mite = ite;
1410 break;
1411 }
1412 else
1413 {
1414 ++ite;
1415 }
1416 }
1417
1418 if( mite == edge_map.end( ) )
1419 {
1420 // 削除できる辺が見つからなかったので終了する
1421 break;
1422 }
1423
1424 difference_type EID = *mite < 0 ? -( *mite ) : *mite;
1425 edge_type &EDGE = edges[ EID ];
1426
1427 #if defined( __SHOW_FACET_DEBUG_INFORMATION__ ) && __SHOW_FACET_DEBUG_INFORMATION__ >= 1
1428 std::cout << "Contraction [" << num_facets << "] : #" << EID << " <" << EDGE.v2 << " -> " << EDGE.v1 << ">" << std::endl;
1429 #endif
1430
1431 // 辺を削除する
1432 edge_map.erase( mite );
1433
1434 // 処理対象の辺の登録と，テーブル内からの削除等を行う
1435 difference_type vs = EDGE.v1;
1436 difference_type vt = EDGE.v2;
1437 difference_type vl, vr;
1438 difference_type eid[ 4 ];
1439 if( !__mc__::compute_edge_connections( faces, edges, EID, vs, vt, vl, vr, eid ) )
1440 {
1441 std::cerr << "Unknown error occured!!" << std::endl;
1442 break;
1443 }
1444
1445 // 面を削除する
1446 face_type &f1 = faces[ EDGE.fid1 ];
1447 face_type &f2 = faces[ EDGE.fid2 ];
1448 f1.flag = false;
1449 f2.flag = false;
1450
1451 // 頂点を移動させる
1452 vertices[ EDGE.v1 ] = EDGE.v;
1453 vertices[ EDGE.v2 ] = EDGE.v;
1454
1455 // 辺の付け替えを行う
1456 __mc__::contract_edges( edges, faces, eid[ 0 ], eid[ 1 ], EID, EDGE.fid1, EDGE.fid1 );
1457 __mc__::contract_edges( edges, faces, eid[ 3 ], eid[ 2 ], EID, EDGE.fid2, EDGE.fid2 );
1458
1459 // 誤差計算行列を更新する
1460 matrix_type QQQ = Q[ EDGE.v1 ] + Q[ EDGE.v2 ];
1461 Q[ EDGE.v1 ] = QQQ;
1462 Q[ EDGE.v2 ] = QQQ;
1463
1464 std::set< difference_type > emap;
1465 {
1466 typename vertex_edge_map_type::iterator ite = vertex_edge_map.find( vs );
1467 if( ite != vertex_edge_map.end( ) )
1468 {
1469 typename vertex_edge_map_type::iterator upper = vertex_edge_map.upper_bound( vs );
1470 for( ; ite != upper ; )
1471 {
1472 if( ite->second == EID )
1473 {
1474 // 削除する
1475 ite = vertex_edge_map.erase( ite );
1476 }
1477 else
1478 {
1479 emap.insert( ite->second );
1480 ++ite;
1481 }
1482 }
1483 }
1484 }
1485 {
1486 typename vertex_edge_map_type::iterator ite = vertex_edge_map.find( vt );
1487 if( ite != vertex_edge_map.end( ) )
1488 {
1489 std::vector< difference_type > combine_edge;
1490 typename vertex_edge_map_type::iterator upper = vertex_edge_map.upper_bound( vt );
1491 for( ; ite != upper ; )
1492 {
1493 if( ite->second == EID )
1494 {
1495 // 削除する
1496 ite = vertex_edge_map.erase( ite );
1497 }
1498 else
1499 {
1500 if( ite->second == eid[ 1 ] || ite->second == eid[ 2 ] )
1501 {
1502 }
1503 else
1504 {
1505 emap.insert( ite->second );
1506 combine_edge.push_back( ite->second );
1507
1508 if( edges[ ite->second ].v1 == EDGE.v2 )
1509 {
1510 edges[ ite->second ].v1 = EDGE.v1;
1511 }
1512 else if( edges[ ite->second ].v2 == EDGE.v2 )
1513 {
1514 edges[ ite->second ].v2 = EDGE.v1;
1515 }
1516 else
1517 {
1518 std::cerr << "Incorrect edge is found!!" << std::endl;
1519 }
1520 }
1521
1522 ite = vertex_edge_map.erase( ite );
1523 }
1524 }
1525
1526 // 不要な辺をすべて削除する
1527 __mc__::remove_edge_from_map( vertex_edge_map, vl, eid[ 1 ] );
1528 __mc__::remove_edge_from_map( vertex_edge_map, vr, eid[ 2 ] );
1529
1530 // 処理対象の枝からも除外する
1531 __mc__::remove_edge_from_set( edge_map, eid[ 1 ] );
1532 __mc__::remove_edge_from_set( edge_map, eid[ 2 ] );
1533
1534 // 変更された辺を再登録する
1535 for( size_type ii = 0 ; ii < combine_edge.size( ) ; ii++ )
1536 {
1537 vertex_edge_map.insert( vertex_edge_map_pair_type( vs, combine_edge[ ii ] ) );
1538 }
1539 }
1540 }
1541
1542 // 各頂点を共有するエッジを更新する
1543 // ただし，統合して領域の端に接した場合は以降の対象から除く
1544 for( typename std::set< difference_type >::iterator ite = emap.begin( ) ; ite != emap.end( ) ; ++ite )
1545 {
1546 __mc__::remove_edge_from_set( edge_map, *ite );
1547
1548 // 縁に接しているかどうかを判定する
1549 if( !__mc__::is_sharp_edge( edges[ *ite ] ) )
1550 {
1551 // 各辺の評価値を更新する
1552 __mc__::update_edge( vertices, Q, edges[ *ite ], use_optimal_vertex_placement );
1553 __mc__::apply_penalties( vertex_edge_map, vertices, faces, edges, *ite, threshold_for_triangle_compactness );
1554 edge_map.insert( *ite );
1555 }
1556 }
1557
1558 #if defined( __SHOW_FACET_DEBUG_INFORMATION__ ) && __SHOW_FACET_DEBUG_INFORMATION__ >= 2
1559 for( vertex_edge_map_type::iterator ite = vertex_edge_map.begin( ) ; ite != vertex_edge_map.end( ) ; ++ite )
1560 {
1561 edge_type &e = edges[ ite->second ];
1562 if( e.fid1 > 0 && !faces[ e.fid1 ].flag )
1563 {
1564 std::cerr << "Edge " << ite->second << " connects to disappeared face." << std::endl;
1565 }
1566 if( e.fid2 > 0 && !faces[ e.fid2 ].flag )
1567 {
1568 std::cerr << "Edge " << ite->second << " connects to disappeared face." << std::endl;
1569 }
1570
1571 difference_type vs = e.v1;
1572 difference_type vt = e.v2;
1573 difference_type vl, vr;
1574 difference_type eid[ 4 ];
1575 if( !__mc__::compute_edge_connections( faces, edges, ite->second, vs, vt, vl, vr, eid ) )
1576 {
1577 std::cerr << "Edge " << ite->second << " is invalid." << std::endl;
1578 }
1579
1580 if( vl == vr )
1581 {
1582 std::cerr << "Edge " << ite->second << " share duplicated faces." << std::endl;
1583 }
1584 }
1585
1586 // 面の接続関係が正しく保たれているかどうかをチェックする
1587 for( size_type i = 1 ; i < faces.size( ) ; i++ )
1588 {
1589 const face_type &f = faces[ i ];
1590 if( f.flag )
1591 {
1592 const difference_type eid1 = f.eid1 < 0 ? -f.eid1 : f.eid1;
1593 const difference_type eid2 = f.eid2 < 0 ? -f.eid2 : f.eid2;
1594 const difference_type eid3 = f.eid3 < 0 ? -f.eid3 : f.eid3;
1595
1596 if( eid1 != 0 && edges[ eid1 ].fid1 != i && edges[ eid1 ].fid2 != i )
1597 {
1598 std::cerr << "Face-edge relationships is broken." << std::endl;
1599 }
1600 if( eid2 != 0 && edges[ eid2 ].fid1 != i && edges[ eid2 ].fid2 != i )
1601 {
1602 std::cerr << "Face-edge relationships is broken." << std::endl;
1603 }
1604 if( eid3 != 0 && edges[ eid3 ].fid1 != i && edges[ eid3 ].fid2 != i )
1605 {
1606 std::cerr << "Face-edge relationships is broken." << std::endl;
1607 }
1608
1609 if( eid1 * eid2 * eid3 != 0 )
1610 {
1611 const edge_type &e1 = edges[ eid1 ];
1612 const edge_type &e2 = edges[ eid2 ];
1613 const edge_type &e3 = edges[ eid3 ];
1614 difference_type v11 = f.eid1 < 0 ? e1.v2 : e1.v1;
1615 difference_type v12 = f.eid1 < 0 ? e1.v1 : e1.v2;
1616 difference_type v21 = f.eid2 < 0 ? e2.v2 : e2.v1;
1617 difference_type v22 = f.eid2 < 0 ? e2.v1 : e2.v2;
1618 difference_type v31 = f.eid3 < 0 ? e3.v2 : e3.v1;
1619 difference_type v32 = f.eid3 < 0 ? e3.v1 : e3.v2;
1620
1621 if( v12 == v21 && v22 == v31 && v32 == v11 )
1622 {
1623 }
1624 else
1625 {
1626 std::cerr << "Face-edge connection is broken." << std::endl;
1627 std::cout << "(" << v11 << ", " << v12 << ") -> (" << v21 << ", " << v22 << ") -> (" << v31 << ", " << v32 << ")" << std::endl;
1628 }
1629 }
1630 }
1631 }
1632 #endif
1633
1634 #if defined( __SHOW_FACET_DEBUG_INFORMATION__ ) && __SHOW_FACET_DEBUG_INFORMATION__ == 3
1635 for( size_type i = 1 ; i < faces.size( ) ; i++ )
1636 {
1637 const face_type &f = faces[ i ];
1638 if( f.flag )
1639 {
1640 std::cout << i << ": ";
1641 if( f.eid1 < 0 )
1642 {
1643 std::cout << edges[ -f.eid1 ].v2 << ", " << edges[ -f.eid1 ].v1 << "[" << -f.eid1 << "] -> ";
1644 }
1645 else
1646 {
1647 std::cout << edges[ f.eid1 ].v1 << ", " << edges[ f.eid1 ].v2 << "[" << f.eid1 << "] -> ";
1648 }
1649
1650 if( f.eid2 < 0 )
1651 {
1652 std::cout << edges[ -f.eid2 ].v2 << ", " << edges[ -f.eid2 ].v1 << "[" << -f.eid2 << "] -> ";
1653 }
1654 else
1655 {
1656 std::cout << edges[ f.eid2 ].v1 << ", " << edges[ f.eid2 ].v2 << "[" << f.eid2 << "] -> ";
1657 }
1658
1659 if( f.eid3 < 0 )
1660 {
1661 std::cout << edges[ -f.eid3 ].v2 << ", " << edges[ -f.eid3 ].v1 << "[" << -f.eid3 << "]" << std::endl;
1662 }
1663 else
1664 {
1665 std::cout << edges[ f.eid3 ].v1 << ", " << edges[ f.eid3 ].v2 << "[" << f.eid3 << "]" << std::endl;
1666 }
1667 }
1668 }
1669 std::cout << std::endl;
1670 #elif defined( __SHOW_FACET_DEBUG_INFORMATION__ ) && __SHOW_FACET_DEBUG_INFORMATION__ >= 4
1671 for( size_type i = 1 ; i < faces.size( ) ; i++ )
1672 {
1673 const face_type &f = faces[ i ];
1674 if( f.flag )
1675 {
1676 printf( "%2d: ", i );
1677
1678 difference_type eid1 = f.eid1 < 0 ? -f.eid1 : f.eid1;
1679 difference_type eid2 = f.eid2 < 0 ? -f.eid2 : f.eid2;
1680 difference_type eid3 = f.eid3 < 0 ? -f.eid3 : f.eid3;
1681 edge_type &e1 = edges[ eid1 ];
1682 edge_type &e2 = edges[ eid2 ];
1683 edge_type &e3 = edges[ eid3 ];
1684
1685 if( f.eid1 < 0 )
1686 {
1687 printf( "%2d, %2d [%2d] (%2d, %2d) -> ", e1.v2, e1.v1, eid1, e1.fid1, e1.fid2 );
1688 }
1689 else
1690 {
1691 printf( "%2d, %2d [%2d] (%2d, %2d) -> ", e1.v1, e1.v2, eid1, e1.fid1, e1.fid2 );
1692 }
1693
1694 if( f.eid2 < 0 )
1695 {
1696 printf( "%2d, %2d [%2d] (%2d, %2d) -> ", e2.v2, e2.v1, eid2, e2.fid1, e2.fid2 );
1697 }
1698 else
1699 {
1700 printf( "%2d, %2d [%2d] (%2d, %2d) -> ", e2.v1, e2.v2, eid2, e2.fid1, e2.fid2 );
1701 }
1702
1703 if( f.eid3 < 0 )
1704 {
1705 printf( "%2d, %2d [%2d] (%2d, %2d)\n", e3.v2, e3.v1, eid3, e3.fid1, e3.fid2 );
1706 }
1707 else
1708 {
1709 printf( "%2d, %2d [%2d] (%2d, %2d)\n", e3.v1, e3.v2, eid3, e3.fid1, e3.fid2 );
1710 }
1711 }
1712 }
1713 std::cout << std::endl;
1714 #endif
1715 }
1716
1717 facets.clear( );
1718
1719 for( size_type i = 1 ; i < faces.size( ) ; i++ )
1720 {
1721 const face_type &f = faces[ i ];
1722 if( f.flag )
1723 {
1724 const difference_type v1 = f.eid1 < 0 ? edges[ -f.eid1 ].v2 : edges[ f.eid1 ].v1;
1725 const difference_type v2 = f.eid2 < 0 ? edges[ -f.eid2 ].v2 : edges[ f.eid2 ].v1;
1726 const difference_type v3 = f.eid3 < 0 ? edges[ -f.eid3 ].v2 : edges[ f.eid3 ].v1;
1727 facets.push_back( facet_type( vertices[ v1 ], vertices[ v2 ], vertices[ v3 ] ) );
1728 }
1729 }
1730
1731 return( true );
1732 }
1733
1741 template < class T >
1742 inline bool maximum_connected_region( facet_list< T > &facets, const double eps = 1.0e-3 )
1743 {
1744 typedef typename facet_list< T >::facet_type facet_type;
1745 typedef typename facet_type::size_type size_type;
1746 typedef typename facet_type::difference_type difference_type;
1747 typedef typename facet_type::vector_type vector_type;
1748 typedef vector3< difference_type > ivector_type;
1749 typedef __mc__::__edge__< vector_type > edge_type;
1750 typedef __mc__::__face__ face_type;
1751 typedef matrix< double > matrix_type;
1752 typedef std::vector< edge_type > edge_list_type;
1753
1754
1755 std::vector< vector_type > vertices;
1756 std::vector< ivector_type > faces;
1757
1758 // 頂点と面のリストに変換する
1759 if( !convert_to_vertex_face_list( facets, vertices, faces, eps ) )
1760 {
1761 return( false );
1762 }
1763
1764
1765 // 頂点と面のテーブルを作成する
1766 typedef std::multimap< size_type, difference_type > vertex_face_map_type;
1767 typedef typename std::multimap< size_type, difference_type >::iterator iterator;
1768 typedef std::pair< size_type, difference_type > vertex_face_map_pair_type;
1769
1770 typedef std::multimap< size_type, size_type > group_map_type;
1771 typedef std::pair< size_type, size_type > group_map_type_pair_type;
1772
1773 vertex_face_map_type vertex_face_map;
1774 group_map_type gmap;
1775 std::vector< size_type > groups( faces.size( ) );
1776 for( size_type i = 0 ; i < faces.size( ) ; i++ )
1777 {
1778 ivector_type &f = faces[ i ];
1779
1780 vertex_face_map.insert( vertex_face_map_pair_type( f.x, i ) );
1781 vertex_face_map.insert( vertex_face_map_pair_type( f.y, i ) );
1782 vertex_face_map.insert( vertex_face_map_pair_type( f.z, i ) );
1783
1784 groups[ i ] = i;
1785 gmap.insert( group_map_type_pair_type( i, i ) );
1786 }
1787
1788 // 各面をグルーピングしていく
1789 for( size_type i = 0 ; i < faces.size( ) ; i++ )
1790 {
1791 ivector_type &f = faces[ i ];
1792 size_type &g1 = groups[ i ];
1793
1794 {
1795 iterator ite = vertex_face_map.find( f.x );
1796 if( ite != vertex_face_map.end( ) )
1797 {
1798 iterator upper = vertex_face_map.upper_bound( f.x );
1799 for( ; ite != upper ; ++ite )
1800 {
1801 if( ite->second != i )
1802 {
1803 size_type &g2 = groups[ ite->second ];
1804
1805 if( g1 < g2 )
1806 {
1807 typename group_map_type::iterator lower = gmap.find( g2 );
1808 typename group_map_type::iterator upper = gmap.upper_bound( g2 );
1809
1810 std::vector< size_type > gtmp;
1811 for( typename group_map_type::iterator gite = lower ; gite != upper ; ++gite )
1812 {
1813 gtmp.push_back( gite->second );
1814 }
1815
1816 gmap.erase( lower, upper );
1817
1818 for( size_type i = 0 ; i < gtmp.size( ) ; i++ )
1819 {
1820 groups[ gtmp[ i ] ] = g1;
1821 gmap.insert( group_map_type_pair_type( g1, gtmp[ i ] ) );
1822 }
1823 }
1824 else if( g1 > g2 )
1825 {
1826 typename group_map_type::iterator lower = gmap.find( g1 );
1827 typename group_map_type::iterator upper = gmap.upper_bound( g1 );
1828
1829 std::vector< size_type > gtmp;
1830 for( typename group_map_type::iterator gite = lower ; gite != upper ; ++gite )
1831 {
1832 gtmp.push_back( gite->second );
1833 }
1834
1835 gmap.erase( lower, upper );
1836
1837 for( size_type i = 0 ; i < gtmp.size( ) ; i++ )
1838 {
1839 groups[ gtmp[ i ] ] = g2;
1840 gmap.insert( group_map_type_pair_type( g2, gtmp[ i ] ) );
1841 }
1842 }
1843 }
1844 }
1845 }
1846 }
1847
1848 {
1849 iterator ite = vertex_face_map.find( f.y );
1850 if( ite != vertex_face_map.end( ) )
1851 {
1852 iterator upper = vertex_face_map.upper_bound( f.y );
1853 for( ; ite != upper ; ++ite )
1854 {
1855 if( ite->second != i )
1856 {
1857 size_type &g2 = groups[ ite->second ];
1858
1859 if( g1 < g2 )
1860 {
1861 typename group_map_type::iterator lower = gmap.find( g2 );
1862 typename group_map_type::iterator upper = gmap.upper_bound( g2 );
1863
1864 std::vector< size_type > gtmp;
1865 for( typename group_map_type::iterator gite = lower ; gite != upper ; ++gite )
1866 {
1867 gtmp.push_back( gite->second );
1868 }
1869
1870 gmap.erase( lower, upper );
1871
1872 for( size_type i = 0 ; i < gtmp.size( ) ; i++ )
1873 {
1874 groups[ gtmp[ i ] ] = g1;
1875 gmap.insert( group_map_type_pair_type( g1, gtmp[ i ] ) );
1876 }
1877 }
1878 else if( g1 > g2 )
1879 {
1880 typename group_map_type::iterator lower = gmap.find( g1 );
1881 typename group_map_type::iterator upper = gmap.upper_bound( g1 );
1882
1883 std::vector< size_type > gtmp;
1884 for( typename group_map_type::iterator gite = lower ; gite != upper ; ++gite )
1885 {
1886 gtmp.push_back( gite->second );
1887 }
1888
1889 gmap.erase( lower, upper );
1890
1891 for( size_type i = 0 ; i < gtmp.size( ) ; i++ )
1892 {
1893 groups[ gtmp[ i ] ] = g2;
1894 gmap.insert( group_map_type_pair_type( g2, gtmp[ i ] ) );
1895 }
1896 }
1897 }
1898 }
1899 }
1900 }
1901
1902
1903 {
1904 iterator ite = vertex_face_map.find( f.z );
1905 if( ite != vertex_face_map.end( ) )
1906 {
1907 iterator upper = vertex_face_map.upper_bound( f.z );
1908 for( ; ite != upper ; ++ite )
1909 {
1910 if( ite->second != i )
1911 {
1912 size_type &g2 = groups[ ite->second ];
1913
1914 if( g1 < g2 )
1915 {
1916 typename group_map_type::iterator lower = gmap.find( g2 );
1917 typename group_map_type::iterator upper = gmap.upper_bound( g2 );
1918
1919 std::vector< size_type > gtmp;
1920 for( typename group_map_type::iterator gite = lower ; gite != upper ; ++gite )
1921 {
1922 gtmp.push_back( gite->second );
1923 }
1924
1925 gmap.erase( lower, upper );
1926
1927 for( size_type i = 0 ; i < gtmp.size( ) ; i++ )
1928 {
1929 groups[ gtmp[ i ] ] = g1;
1930 gmap.insert( group_map_type_pair_type( g1, gtmp[ i ] ) );
1931 }
1932 }
1933 else if( g1 > g2 )
1934 {
1935 typename group_map_type::iterator lower = gmap.find( g1 );
1936 typename group_map_type::iterator upper = gmap.upper_bound( g1 );
1937
1938 std::vector< size_type > gtmp;
1939 for( typename group_map_type::iterator gite = lower ; gite != upper ; ++gite )
1940 {
1941 gtmp.push_back( gite->second );
1942 }
1943
1944 gmap.erase( lower, upper );
1945
1946 for( size_type i = 0 ; i < gtmp.size( ) ; i++ )
1947 {
1948 groups[ gtmp[ i ] ] = g2;
1949 gmap.insert( group_map_type_pair_type( g2, gtmp[ i ] ) );
1950 }
1951 }
1952 }
1953 }
1954 }
1955 }
1956 }
1957
1958 size_type mgroup = 0;
1959 size_type mcount = 0;
1960 for( typename group_map_type::iterator ite = gmap.begin( ) ; ite != gmap.end( ) ; )
1961 {
1962 size_type count = gmap.count( ite->first );
1963 if( count > mcount )
1964 {
1965 mcount = count;
1966 mgroup = ite->first;
1967 }
1968
1969 ite = gmap.upper_bound( ite->first );
1970 if( ite == gmap.end( ) )
1971 {
1972 break;
1973 }
1974 else
1975 {
1976 ++ite;
1977 }
1978 }
1979
1980 facets.clear( );
1981
1982 typename group_map_type::iterator lower = gmap.find( mgroup );
1983 if( lower != gmap.end( ) )
1984 {
1985 typename group_map_type::iterator upper = gmap.upper_bound( mgroup );
1986 for( ; lower != upper ; ++lower )
1987 {
1988 ivector_type &f = faces[ lower->second ];
1989 facets.push_back( facet_type( vertices[ f.x ], vertices[ f.y ], vertices[ f.z ] ) );
1990 }
1991 }
1992
1993 return( true );
1994 }
1995
1997 // Marching Cubesグループの終わり
1998
1999
2000 // mist名前空間の終わり
2001 _MIST_END
2002
2003
2004 #endif // __INCLUDE_MIST_FACET__
2005