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;
}