基于osgEarth的卫星瞬时覆盖绘制方法

---

一、问题提出

卫星是空间态势中的重要对象，在空间态势中不仅要表现卫星的位置、轨道信息，还需要表现卫星的覆盖能力。

表现卫星覆盖包括2个问题：卫星探测在空间的形状，在场景中一般采用具有一定透明度的几何体表示；卫星在地球表面形成的覆盖范围，用离散的线集表示。

如在STK中，卫星可附着sensor对象，其覆盖形状可为简单圆锥、四棱锥、SAR等。sensor的设置非常灵活，可以与卫星存在各种位置关系，而不是一定在卫星坐标系原点。同时当这些sensor的形状超出地球表面范围时，可以仍然只显示相交区域。

OSG本身是既支持预设的多种几何形状，也可以直接构造primitiveset，如以三角形扇表示圆锥。我们希望获得如下图所示的绘制效果。

但如果不进行任何处理直接以三角形扇绘制，则随着视点位置的改变，表示覆盖的透明区域会出现下图所示的情况。

圆锥中出现了部分重复绘制的部分，即透明绘制是不稳定的。究其原因，是透明覆盖本身是存在先后关系的，圆锥是由多个三角形表示的，在不进行任何处理的情况下其绘制顺序固定。如果靠近视点一侧的三角形都先于远离视点一侧三角形，则远离一侧三角形在绘制时由于Z缓冲值已经小于当前值，所以不会绘制，形成理想的效果；否则，远离一侧的部分三角形先绘制完成后，靠近一侧三角形还会被绘制，形成上图效果。

需要指出的是，上述效果的出现前提是：绘制透明覆盖时不进行背面剔除（CULL_FACE）。如果进行背面剔除，那么通过控制三角形方向，对于简单圆锥形状的透明覆盖，是不会出现上述现象的。但如果是复杂的覆盖形状，仍会出现上述现象；同时，当视点进入覆盖内部的时候，应该可见的面不会绘制。

对于地表形成的覆盖区域边界，可以用LINE_LOOP表示，但如果直接表示，由于离散三角形与地表的交点构成的线段与表示地表的网格并不共面，因此会出现下图所示的消隐错误情况（显示不连贯、闪烁等），而且这是不能通过polygon_offset解决的（线与地球表面三角形并不共面）。当然，在计算时可以不是与地球表面求交，而是与地球上方如8km高程的椭球求交，那么可以避免消隐不正确，但是显然当视点靠近地表时其效果明显不好。

---

二、控制三角形扇的绘制顺序实现透明覆盖的稳定绘制

对于透明覆盖在不同视点时表现不稳定问题，应该可以通过shader解决，暂时还未尝试，现在介绍另外一种根据视点调整三角形扇绘制顺序的方法。

1.三角形扇数据生成

如图所示，根据卫星位置（严格来说是传感器位置，其原点可位于卫星坐标系中某点）、指向以及覆盖角度，采用离散的三角形扇表示覆盖锥形（如固定120个三角形或者间隔若干度形成一个三角形）。此时，需要计算离散后的圆锥与地球表面的交点，这些交点既用于表示三角形扇，也用于表示地表覆盖区域。虽然osgEarth支持地形数据，但如果与地形实时求交，其计算量过大。因此与地球椭球计算交点。交点计算通过矢量运算完成，代码如下

	CVector3D pos = _frame->getOrigin();
	osg::Vec3Array* va0 = dynamic_cast<osg::Vec3Array*>(_geometryPlane->getVertexArray());
	(*va)[0] = osg::Vec3(pos.x, pos.y, pos.z);
	osg::Vec3Array* va1 = dynamic_cast<osg::Vec3Array*>(_geometryLoop->getVertexArray());
	CVector3D xAxis = _frame->getAxisX();
	CVector3D yAxis = _frame->getAxisY();
	CVector3D zAxis = _frame->getAxisZ();
	double dis = 1000.0;
	double tmp = 1000.0*tan(_detectAngle* osg::PI / 180.0);
	for (int i = 0; i < discreteNum; i++){
		CVector3D p1 = pos + xAxis * tmp * _sincos[2 * i] + yAxis * tmp * _sincos[2 * i + 1] + zAxis * dis;
		CVector3D p;
		segInterEarth(pos, p1, p);
		(*va)[i + 1] = osg::Vec3(p.x, p.y, p.z);
		(*va1)[i] = osg::Vec3(p.x, p.y, p.z);}
	(*va)[discreteNum + 1] = (*va)[1];
	(*va1)[discreteNum] = (*va1)[0];

为优化效率，节点所需数据空间预先分配，所需三角函数计算预先计算并存储，数据更新只有当仿真时间改变、卫星位置改变时才发生。

2.根据视点生成图元，控制三角形绘制顺序

如图所示，原始的数据是由a点开始的若干个点（按顺时针顺序），加上卫星位置点，构成三角形扇表示空间覆盖。如果按固定的图元绘制，设由a开始顺时针构成三角形扇绘制顺序。

当视点处于图中位置A时，绘制m到n之间的三角形时，由a到m之间的三角形已经绘制完成，因此就会重复绘制；而n之后的三角形，由于z缓冲区中值大于其对应深度，因此不会绘制。因此，就形成了前述不稳定绘制的现象。视点位于图中B时，也可作类似分析。但当视点位于图中位置C时，由于可见的三角形先于不可见三角形绘制，因此绘制效果为合理的效果。

此处给出的方法是根据视点位置构造绘制图元，控制三角形的绘制顺序，确保可见的三角形先于不可见三角形绘制。

如上图中视点位于A时，先构造由m到最后一个点的三角形扇，再构造由a到m的三角形扇，也可以由m开始构造一个三角形扇（由于数据都已经预先生成，所以构造以索引方式描述的几何图元）。

当视点位于B时，也可以先构造j到i的三角形扇，然后构造i到j的三角形扇，也可以由j开始构造一个三角形扇。

具体算法为：
1）判断第一个三角形是否可见；
2）循环判断所有三角形的可见性，记录可见性发生改变的位置，即为上图中i、j、m、n之类，可见性发生改变的点或为0个，或为2个；
3）根据第一个三角形的可见性情况，按顺序构造三角形扇。
其中，三角形可见性判断：由于是一个简单圆锥，因此在构造数据时，令三角形向内一侧法向为正，构造卫星位置到视点的矢量，通过该矢量与三角形法向的数量积来判断三角形是否可见。

	CVector3D pos = _frame->getOrigin();
	osg::ref_ptr<osg::Vec3Array> va = dynamic_cast<osg::Vec3Array*>(_geometryPlane->getVertexArray());
	osgViewer::Viewer* viewer = g_pDetectPerformanceEvaluation->getSpaceWindow3D()->getOsgViewer();
	osg::Vec3d eye, center, up;
	viewer->getCamera()->getViewMatrixAsLookAt(eye, center, up, 1000);
	osg::Vec3d v0 = va->at(0);
	osg::Vec3d v1 = va->at(1);
	osg::Vec3d v2 = va->at(2);
	up = eye - v0;
	up.normalize();
	v1 -= v0;
	v2 -= v0;
	osg::Vec3d vv = v1 ^ v2;
	vv.normalize();
	double tmp = up * vv;
	if (tmp < 0)	{
		int i = 2;
		for (; i < discreteNum + 1; i++){
			v1 = va->at(i);
			v2 = va->at(i + 1);
			v1 -= v0;
			v2 -= v0;
			vv = v1 ^ v2;
			vv.normalize();
			double tmp1 = up * vv;
			if (tmp1 > 0)
				break;	}
		if (i == discreteNum + 1)	//无可见性改变点{
			_geometryPlane->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::TRIANGLE_FAN, 0, discreteNum + 2));
			return;	}
		int pos =  i;
		i++;
		for (; i < discreteNum + 1; i++){
			v1 = va->at(i);
			v2 = va->at(i + 1);
			v1 -= v0;
			v2 -= v0;
			vv = v1 ^ v2;
			vv.normalize();
			double tmp1 = up * vv;
			if (tmp1 < 0)
				break;	}
		osg::ref_ptr<osg::DrawElementsUShort> primitive1 = new osg::DrawElementsUShort(GL_TRIANGLE_FAN);
		primitive1->push_back(0);
		int j;
		for (j = i; j <= discreteNum + 1; j++)
			primitive1->push_back(j);
		for (j = 2; j <= pos; j++)
			primitive1->push_back(j);
		_geometryPlane->addPrimitiveSet(primitive1);
		osg::ref_ptr<osg::DrawElementsUShort> primitive2 = new osg::DrawElementsUShort(GL_TRIANGLE_FAN);
		primitive2->push_back(0);
		for (; j <= i; j++)
			primitive2->push_back(j);
		_geometryPlane->addPrimitiveSet(primitive2);	}
	else//第一个三角形不可见，略

---

三、强制计算相交线的可见性实现地表区域线的稳定绘制

离散后的圆锥相当于由卫星位置发出的射线与地球表面求交，交点连线与地球表面的位置关系不共面，因此消隐不正确。为此，根据视点位置，实时计算出各交点的可见性，可见的线不受Z缓冲消隐影响，强制绘制。

交点可见性判断方法：构造由视点到交点的射线（参数方程表示，则视点对应参数0，交点参数为1），该射线与地球表面有2个交点，其参数为t0和t1，且t0<t1；如果t0<1，则证明另一个交点更靠近视点，因此不可见；否则可见。

如图所示，a为构成卫星地表覆盖的离散交点，当视点位于A、B时，与a共线，此时直线AB与地球表面的另一个交点为b。a所对应的参数方程参数为1。当视点位于B时，两个参数中小者为1，可见；当视点位于视点A时，两个参数中小者小于1，不可见。

osg::ref_ptr<osg::Vec3Array> va = dynamic_cast<osg::Vec3Array*>(_geometryLoop->getVertexArray());
	osgViewer::Viewer* viewer = Application::_viewer3D;
	osg::Vec3d eye, center, up;
	viewer->getCamera()->getViewMatrixAsLookAt(eye, center, up, 1000);
	osg::Vec3d v0 = va->at(0);
	CPoint3D eyep(eye.x(), eye.y(), eye.z());
	CPoint3D p0(v0.x(), v0.y(), v0.z());
	bool flag = isFacedEye(eyep, p0);
	if (flag)	//可见	{
		int i = 1;
		for (; i < discreteNum; i++)	{
			osg::Vec3d v1 = va->at(i);
			CPoint3D p1(v1.x(), v1.y(), v1.z());
			bool flag1 = isFacedEye(eyep, p1);
			if (!flag1)break;}
		if (i == discreteNum)
			_geometryLoop->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_LOOP, 0, discreteNum));
		else{
			if (i - 1 > 0)
				_geometryLoop->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP, 0, i));
			for (; i < discreteNum; i++)	{
				osg::Vec3d v1 = va->at(i);
				CPoint3D p1(v1.x(), v1.y(), v1.z());
				bool flag1 = isFacedEye(eyep, p1);
				if (flag1)break;}
			//此时，i又开始可见，反向构造图元
			if (discreteNum - i > 0)
			_geometryLoop->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::LINE_STRIP, i - 1, discreteNum + 2 - i));}}
	else//不可见，略

对于表示覆盖边线的图元，需强制深度测试通过。