DOPの描画(4) - 古いコード -
DOPのrealignmentにはDOPの頂点が必要なのですが、描画の時に使った方法(全ての3平面の交点を列挙し、その中から頂点を求める)ではどうにも効率が悪すぎます。そこで、おとなしくqhull(qhalf)を使うことにしました。
ちなみに、DOPのように平面で構成されたものは、日本語では「半空間の積」と呼ばれるようです。計算幾何学の本をパラパラ読んでみたら乗ってました。qhullのように半空間の積を求めるには、「双対点(dual point)」という言葉がキーワードになりそうです。そのうち勉強しよう。
qhullで全て求めることにしたので、前に作ったコードは必要ありません。ちょっともったいないのでコードを残しておきます。
//================================================================= // Func: CreateHSpaceIntersection // 平面からそれらの交差を元に多面体を作る // Params: // planes - 平面(正規化されている必要がある) // numPlanes - 平面の数 // verts - 作った多面体の頂点 // idxs - 作った多面体の面構成 //================================================================= bool CreateHSpaceIntersection( const Plane *planes, int numPlanes, std::vector<Vector3d> &verts, std::vector<int> &idxs ) { verts.clear(); idxs.clear(); // 頂点となる点を求める std::vector<Vector3d> points; for ( int i=0; i<numPlanes; ++i ) { for ( int j=i+1; j<numPlanes; ++j ) { for ( int k=j+1; k<numPlanes; ++k ) { // 三つの面の交点が頂点の候補となる Vector3d p; bool res = IntersectPlanes( planes[i], planes[j], planes[k], p ); // 三つの面の交点pのうち、全ての平面内にあるものだけが頂点となる for ( int l=0; l<numPlanes && res; ++l ) { if ( DistPointPlane(p, planes[l]) > EPSILON ) res=false; } if ( res ) { // 既に登録されている頂点と同じ頂点は登録しない for ( int l=0; l<points.size(); ++l ) { if ( IsEqual( p, points[l] ) ) break; } if ( l==points.size() ) points.push_back( Vector3d(p) ); } } } } // 同じ頂点と思われるものは削除(頂点間の距離が平均の0.01%よりも短いものを削除の対象とする) double sumDist = 0; for ( i=0; i<points.size(); ++i ) { for ( int j=i+1; j<points.size(); ++j ) sumDist += Distance( points[i], points[j] ); } double avgDist = sumDist / (points.size()*(points.size()-1))/2.0; i=0; while ( i<points.size() ) { int delIdx = -1; for ( i=0; i<points.size(); ++i ) { int j; for ( j=i+1; j<points.size(); ++j ) { if ( Distance(points[i],points[j]) < avgDist*0.01 ) { delIdx = i; break; } } if ( j<points.size() ) break; } if ( i<points.size() ) points.erase( &points[delIdx] ); } // 多面体であるためには4頂点以上必要 if ( points.size() < 4 ) return true; // 全ての点が同一平面上にあるかどうかチェック(同一平面にある場合は凸包は行わない) Plane pl = ComputePlane( points[0], points[1], points[2] ); for ( i=3; i<points.size(); ++i ) { if ( !IsZero( DistPointPlane( points[i], pl ) ) ) break; } if ( i == points.size() ) return true; // 多面体にならないので凸包を行わない // 凸包の実行 if ( qh_new_qhull( 3, points.size(), points[0].vertex, False, "qhull ", NULL, stderr) ) return false; // 凸包に失敗 // 頂点を取得 pointT *point, *pointtemp; FORALLpoints { verts.push_back( Vector3d(point[0],point[1],point[2]) ); } // 面を取得 facetT *facet; vertexT *vertex, **vertexp; FORALLfacets { setT *vertices = qh_facet3vertex(facet); // 逆時計回り順にする idxs.push_back(qh_setsize(vertices)); FOREACHvertex_(vertices) { idxs.push_back( qh_pointid(vertex->point) ); } qh_settempfree(&vertices); } qh_freeqhull(!qh_ALL); // free long memory int curlong, totlong; // memory remaining after qh_memfreeshort qh_memfreeshort (&curlong, &totlong); // free short memory and memory allocator if (curlong || totlong) { fprintf (stderr, "qhull internal warning (user_eg, #1): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); return false; } return true; }