OsgEarth下量取球面面积

#ifndef AREA_H
#define AREA_H

#include "handleadapter.h"

class Area : public HandleAdapter
{
public:
	Area(GraphicsView* view);
	~Area();

protected:
	virtual void slotPicked(osg::Vec3d pos);
	virtual void slotMoveing(osg::Vec3d pos);
	virtual void slotRightHandle();

private:
	void drawPolygon(osg::Vec3d pos = osg::Vec3d());

private:
	std::vector<osg::Vec3d> m_vecPoints;
	osgEarth::Symbology::Style m_polygonStyle;
	osgEarth::Symbology::Style m_textStyle;
	osgEarth::Annotation::FeatureNode* m_pFeatureNode;
	osgEarth::Annotation::FeatureEditor* m_pPolygonEdit;
	osgEarth::Annotation::PlaceNode* m_pPlaceNode;

	/添加橡皮筋功能
	osgEarth::Symbology::Style m_stippleLineStyle;
	osgEarth::Annotation::FeatureNode* m_pStippleFeatureNode;
};

#endif // AREA_H

#include "area.h"
#include "calcmath.h"

Area::Area(GraphicsView* view)
: HandleAdapter(view)
{
	m_vecPoints.clear();
	m_pFeatureNode = NULL;
	//m_pPolygonEdit = NULL;
	m_pPlaceNode = NULL;
	m_pStippleFeatureNode = NULL;


	m_polygonStyle.getOrCreate<osgEarth::Symbology::PolygonSymbol>()->fill()->color() = osgEarth::Symbology::Color::Yellow;
	m_polygonStyle.getOrCreate<osgEarth::Symbology::AltitudeSymbol>()->clamping() = osgEarth::Symbology::AltitudeSymbol::CLAMP_TO_TERRAIN;
	m_polygonStyle.getOrCreate<osgEarth::Symbology::AltitudeSymbol>()->technique() = osgEarth::Symbology::AltitudeSymbol::TECHNIQUE_DRAPE;
	m_polygonStyle.getOrCreate<osgEarth::Symbology::AltitudeSymbol>()->verticalOffset() = 0.1;
	/*********橡皮筋**********/
	m_stippleLineStyle.getOrCreate<osgEarth::Symbology::LineSymbol>()->stroke()->color() = osgEarth::Symbology::Color::Red;
	m_stippleLineStyle.getOrCreate<osgEarth::Symbology::LineSymbol>()->stroke()->width() = 2.0;
	m_stippleLineStyle.getOrCreate<osgEarth::Symbology::LineSymbol>()->tessellation() = 20.0;
	m_stippleLineStyle.getOrCreate<osgEarth::Symbology::AltitudeSymbol>()->clamping() = osgEarth::Symbology::AltitudeSymbol::CLAMP_TO_TERRAIN;
	m_stippleLineStyle.getOrCreate<osgEarth::Symbology::AltitudeSymbol>()->technique() = osgEarth::Symbology::AltitudeSymbol::TECHNIQUE_DRAPE;
	m_stippleLineStyle.getOrCreate<osgEarth::Symbology::AltitudeSymbol>()->verticalOffset() = 0.1;
	m_stippleLineStyle.getOrCreate<osgEarth::Symbology::LineSymbol>()->stroke()->stipple() = 255;

	m_textStyle.getOrCreate<osgEarth::Symbology::TextSymbol>()->alignment()= osgEarth::Symbology::TextSymbol::ALIGN_LEFT_TOP;
	m_textStyle.getOrCreate<osgEarth::Symbology::TextSymbol>()->encoding()= osgEarth::Symbology::TextSymbol::ENCODING_UTF8;
	m_textStyle.getOrCreate<osgEarth::Symbology::TextSymbol>()->declutter()	= true;
	m_textStyle.getOrCreate<osgEarth::Symbology::TextSymbol>()->halo()= osgEarth::Symbology::Color::White;
	m_textStyle.getOrCreate<osgEarth::Symbology::TextSymbol>()->haloOffset() = 0.05;
	m_textStyle.getOrCreate<osgEarth::Symbology::TextSymbol>()->fill()->color()	= osg::Vec4(1, 0.5, 0, 1);
	m_textStyle.getOrCreate<osgEarth::Symbology::RenderSymbol>()->lighting() = false;
	m_textStyle.getOrCreate<osgEarth::Symbology::TextSymbol>()->font() = "simhei.ttf";
	m_textStyle.getOrCreate<osgEarth::Symbology::TextSymbol>()->size() = 25.0;
}

Area::~Area()
{

}

void Area::slotPicked(osg::Vec3d pos)
{
	/*m_vecPoints.push_back(pos);
	if (m_vecPoints.size() <= 2)
	{
		return;
	}

	drawPolygon();*/
	
	m_vecPoints.push_back(pos);
	if (m_vecPoints.size() <= 1)
	{
		return;
	}

	/*if (m_pPolygonEdit != NULL)
	{
		m_pPolygonEdit->removeChildren(0, m_pPolygonEdit->getNumChildren());
		m_pPolygonEdit = NULL;
	}*/
	if (m_pFeatureNode == NULL)
	{
		osgEarth::Features::Feature* pFeature = new osgEarth::Features::Feature(new osgEarth::Symbology::Polygon, m_pMap3D->getSRS(), m_polygonStyle);
		m_pFeatureNode = new osgEarth::Annotation::FeatureNode(m_pMap3D->getMapNode(), pFeature);

		m_pPlaceNode = new osgEarth::Annotation::PlaceNode(m_pMap3D->getMapNode(), osgEarth::GeoPoint::GeoPoint(), "", m_textStyle);
		m_pPlaceNode->setDynamic(true);

		m_pLayerGroup->addChild(m_pFeatureNode);
		m_pLayerGroup->addChild(m_pPlaceNode);
		
	}
	osgEarth::Symbology::Geometry* pGeometry = m_pFeatureNode->getFeature()->getGeometry();
	pGeometry->clear();
	m_pFeatureNode->setStyle(m_polygonStyle);
	for (auto& n : m_vecPoints)
	{
		pGeometry->push_back(n);
	
	}	
	CalcMath math;
	QString str = QStringLiteral("面积：%1 ㎡").arg(QString::number(math.calcArea(m_vecPoints), 'g', 9));
	m_pPlaceNode->setText(str.toStdString());
	m_pPlaceNode->setPosition(osgEarth::GeoPoint::GeoPoint(m_pMap3D->getSRS(), pos));
	///
	m_pFeatureNode->init();

	if (m_pStippleFeatureNode != NULL)
	{
		m_pStippleFeatureNode->getFeature()->getGeometry()->clear();
		m_pStippleFeatureNode->init();
	}

	/*if (m_pPolygonEdit == NULL)
	{
		m_pPolygonEdit = new osgEarth::Annotation::FeatureEditor(m_pFeatureNode);
		m_pLayerGroup->addChild(m_pPolygonEdit);
	}*/
}

void Area::slotMoveing(osg::Vec3d pos)
{
	/*if (m_vecPoints.size() < 2)
	{
		return;
	}
	drawPolygon(pos);*/
	if (m_vecPoints.size() < 2)
	{
		return;
	}
	if (m_pStippleFeatureNode == NULL)
	{
		osgEarth::Features::Feature* pFeature = new osgEarth::Features::Feature(new osgEarth::Annotation::LineString, m_pMap3D->getSRS(), m_stippleLineStyle);
		m_pStippleFeatureNode = new osgEarth::Annotation::FeatureNode(m_pMap3D->getMapNode(), pFeature);
		m_pLayerGroup->addChild(m_pStippleFeatureNode);
	
	}
	std::vector<osg::Vec3d> vecPoints = m_vecPoints;//??
	osgEarth::Symbology::Geometry* pGeometry = m_pStippleFeatureNode->getFeature()->getGeometry();
	pGeometry->clear();
	m_pStippleFeatureNode->setStyle(m_stippleLineStyle);
	pGeometry->push_back(m_vecPoints[0]);
	pGeometry->push_back(pos);
	pGeometry->push_back(m_vecPoints[m_vecPoints.size() - 1]);

	//
	if (pos != osg::Vec3d())//?????
	{
		//pGeometry->push_back(pos);
		vecPoints.push_back(pos);
	}

	CalcMath math;
	QString str = QStringLiteral("面积：%1 ㎡").arg(QString::number(math.calcArea(vecPoints), 'g', 9));
	m_pPlaceNode->setText(str.toStdString());
	m_pPlaceNode->setPosition(osgEarth::GeoPoint::GeoPoint(m_pMap3D->getSRS(), pos));

	m_pStippleFeatureNode->init();
}

void Area::slotRightHandle()
{
	
	if (m_pStippleFeatureNode != NULL)
	{
		m_pStippleFeatureNode->getFeature()->getGeometry()->clear();
		m_pStippleFeatureNode->init();
		CalcMath math;
		QString str = QStringLiteral("面积：%1 ㎡").arg(QString::number(math.calcArea(m_vecPoints), 'g', 9));
		m_pPlaceNode->setText(str.toStdString());
		m_pPlaceNode->setPosition(osgEarth::GeoPoint::GeoPoint(m_pMap3D->getSRS(), m_vecPoints[0]));
	}
	if (m_pFeatureNode!= nullptr)
	{
		m_pFeatureNode = nullptr;
		m_vecPoints.clear();
	}	
	
}

#ifndef MATH_H
#define MATH_H

#include "link.h"
class CalcMath
{
public:
	CalcMath();
	~CalcMath();

	// 球面面积计算公式
	double calcArea(std::vector<osg::Vec3d> points);

private:
	
};

#endif // MATH_H

#include "calcmath.h"

CalcMath::CalcMath()
{

}

CalcMath::~CalcMath()
{

}

double CalcMath::calcArea(std::vector<osg::Vec3d> points)
{
	double dTotalArea = 0.0;
	int iCount = points.size();
	if (iCount > 2)
	{
		double dLowX = 0.0;
		double dLowY = 0.0;
		double dMiddleX = 0.0;
		double dMiddleY = 0.0;
		double dHighX = 0.0;
		double dHighY = 0.0;
		double AM = 0.0;
		double BM = 0.0;
		double CM = 0.0;
		double AL = 0.0;
		double BL = 0.0;
		double CL = 0.0;
		double AH = 0.0;
		double BH = 0.0;
		double CH = 0.0;
		//Coefficient系数
		double dCoefficientL = 0.0;
		double dCoefficientH = 0.0;
		//tangent切线
		double dALtangent = 0.0;
		double dBLtangent = 0.0;
		double dCLtangent = 0.0;
		double dAHtangent = 0.0;
		double dBHtangent = 0.0;
		double dCHtangent = 0.0;
		//NormalLine法线     
		double dANormalLine = 0.0;
		double dBNormalLine = 0.0;
		double dCNormalLine = 0.0;
		//Orientation  Value方向值  
		double dOrientationValue = 0.0;
		//余弦角
		double dAngleCos = 0.0;

		double dSum1 = 0.0;
		double dSum2 = 0.0;
		int iCount2 = 0;
		int iCount1 = 0;

		double dSum = 0.0;

		for (int i = 0; i < iCount; i++)
		{
			if (i == 0)
			{
				//换算成弧度;
				dLowX = osg::DegreesToRadians(points.at(iCount - 1).x());
				dLowY = osg::DegreesToRadians(points.at(iCount - 1).y());
				dMiddleX = osg::DegreesToRadians(points.at(0).x());
				dMiddleY = osg::DegreesToRadians(points.at(0).y());
				dHighX = osg::DegreesToRadians(points.at(1).x());
				dHighY = osg::DegreesToRadians(points.at(1).y());
			}
			else if (i == iCount - 1)
			{
				dLowX = osg::DegreesToRadians(points.at(iCount - 2).x());
				dLowY = osg::DegreesToRadians(points.at(iCount - 2).y());
				dMiddleX = osg::DegreesToRadians(points.at(iCount - 1).x());
				dMiddleY = osg::DegreesToRadians(points.at(iCount - 1).y());
				dHighX = osg::DegreesToRadians(points.at(0).x());
				dHighY = osg::DegreesToRadians(points.at(0).y());
			}
			else
			{
				dLowX = osg::DegreesToRadians(points.at(i - 1).x());
				dLowY = osg::DegreesToRadians(points.at(i - 1).y());
				dMiddleX = osg::DegreesToRadians(points.at(i).x());
				dMiddleY = osg::DegreesToRadians(points.at(i).y());
				dHighX = osg::DegreesToRadians(points.at(i + 1).x());
				dHighY = osg::DegreesToRadians(points.at(i + 1).y());
			}
			AM = cos(dMiddleY) * cos(dMiddleX);
			BM = cos(dMiddleY) * sin(dMiddleX);
			CM = sin(dMiddleY);
			AL = cos(dLowY) * cos(dLowX);
			BL = cos(dLowY) * sin(dLowX);
			CL = sin(dLowY);
			AH = cos(dHighY) * cos(dHighX);
			BH = cos(dHighY) * sin(dHighX);
			CH = sin(dHighY);

			dCoefficientL = (AM*AM + BM*BM + CM*CM) / (AM*AL + BM*BL + CM*CL);
			dCoefficientH = (AM*AM + BM*BM + CM*CM) / (AM*AH + BM*BH + CM*CH);

			dALtangent = dCoefficientL * AL - AM;
			dBLtangent = dCoefficientL * BL - BM;
			dCLtangent = dCoefficientL * CL - CM;
			dAHtangent = dCoefficientH * AH - AM;
			dBHtangent = dCoefficientH * BH - BM;
			dCHtangent = dCoefficientH * CH - CM;

			dAngleCos = (dAHtangent * dALtangent + dBHtangent * dBLtangent + dCHtangent * dCLtangent) /
				(sqrt(dAHtangent * dAHtangent + dBHtangent * dBHtangent + dCHtangent * dCHtangent) *
				sqrt(dALtangent * dALtangent + dBLtangent * dBLtangent + dCLtangent * dCLtangent));

			dAngleCos = acos(dAngleCos);

			dANormalLine = dBHtangent * dCLtangent - dCHtangent * dBLtangent;
			dBNormalLine = 0 - (dAHtangent * dCLtangent - dCHtangent * dALtangent);
			dCNormalLine = dAHtangent * dBLtangent - dBHtangent * dALtangent;

			if (AM != 0)
			{
				dOrientationValue = dANormalLine / AM;
			}
			else if (BM != 0)
			{
				dOrientationValue = dBNormalLine / BM;
			}
			else
			{
				dOrientationValue = dCNormalLine / CM;
			}
			if (dOrientationValue > 0)
			{
				dSum1 += dAngleCos;
				iCount1++;
			}
			else
			{
				dSum2 += dAngleCos;
				iCount2++;
			}
		}
		if (dSum1 > dSum2)
		{
			dSum = dSum1 + (2 * osg::PI*iCount2 - dSum2);
		}
		else
		{
			dSum = (2 * osg::PI*iCount1 - dSum1) + dSum2;
		}

		dTotalArea = (dSum - (iCount - 2)*osg::PI)* osg::WGS_84_RADIUS_EQUATOR *
			osg::WGS_84_RADIUS_EQUATOR;
	}
	return dTotalArea;
}