地球楕円面上の多角形面積の計算(C++コード)

昨日突然テストした时に以前の制品の中で书いた地球の楕球面上の面积の计算のコードが少し问题があることを発见して、そこで今日彻底的に修正して、QGISの中からコードをほじくり出してC++で书き直して、新しいコードは比较的に正确に楕球面上の多角形の面积を计算することができて、この基础関数は空间量算机能の中の面积の量测に対してとても重要で、ここで共有して参考にしたり、直接持って行ったりします.
ヘッダファイルは次のとおりです.

/**
* @file              DistanceArea.h
* @brief                             
* @details           
* @author            zxg
* @date              2015 5 15 
* @version           1.0.0.1
* @par               Copyright (c):   
* @par               History:
*/
#ifndef __DISTANCEAREA_H_C0C6C4C8_B1C4_49C3_87E4_EEF1ABC862F0__
#define __DISTANCEAREA_H_C0C6C4C8_B1C4_49C3_87E4_EEF1ABC862F0__




class DistanceArea
{
public:
	DistanceArea();

    ~DistanceArea();

	/**
	* @brief SetEllipsePara          
	* @param[in] double a    
	* @param[in] double b    
	* @return void
	* @author zxg
	* @date 2015 5 15 
	* @note 
	*/
	void SetEllipsePara(double a,double b);

	/**
	* @brief ComputePolygonArea               
	* @param[in]  const double *padX X    
	* @param[in] const double* padY Y    
	* @param[in] int nCount     
	* @return double       
	* @author zxg
	* @date 2015 5 15 
	* @note 
	*/
	double ComputePolygonArea( const double *padX,const double* padY,int nCount );

private:

    double mSemiMajor, mSemiMinor, mInvFlattening;

    double GetQ( double x );
    double GetQbar( double x );

	void ComputeAreaInit();

    //         

    double m_QA, m_QB, m_QC;
    double m_QbarA, m_QbarB, m_QbarC, m_QbarD;
    double m_AE;  /* a^2(1-e^2) */
    double m_Qp;  
    double m_E;  
    double m_TwoPI;

};


#endif // end of __DISTANCEAREA_H_C0C6C4C8_B1C4_49C3_87E4_EEF1ABC862F0__

実装コードは次のとおりです.

#include "DistanceArea.h"

#define DEG2RAD(x)    ((x)*M_PI/180)


DistanceArea::DistanceArea()
{
	//
}


DistanceArea::~DistanceArea()
{
}

void DistanceArea::SetEllipsePara(double a,double b)
{
	mSemiMajor = a;
	mSemiMinor = b;
	//mInvFlattening = mSemiMajor

	ComputeAreaInit();
}

double DistanceArea::GetQ( double x )
{
	double sinx, sinx2;

	sinx = sin( x );
	sinx2 = sinx * sinx;

	return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
}


double DistanceArea::GetQbar( double x )
{
	double cosx, cosx2;

	cosx = cos( x );
	cosx2 = cosx * cosx;

	return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
}


void DistanceArea::ComputeAreaInit()
{
	double a2 = ( mSemiMajor * mSemiMajor );
	double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
	double e4, e6;

	m_TwoPI = M_PI + M_PI;

	e4 = e2 * e2;
	e6 = e4 * e2;

	m_AE = a2 * ( 1 - e2 );

	m_QA = ( 2.0 / 3.0 ) * e2;
	m_QB = ( 3.0 / 5.0 ) * e4;
	m_QC = ( 4.0 / 7.0 ) * e6;

	m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4  - ( 4.0 / 7.0 ) * e6;
	m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4  + ( 4.0 / 7.0 ) * e6;
	m_QbarC =                     - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
	m_QbarD = ( 4.0 / 49.0 ) * e6;

	m_Qp = GetQ( M_PI / 2 );
	m_E  = 4 * M_PI * m_Qp * m_AE;
	if ( m_E < 0.0 )
	m_E = -m_E;
}


double DistanceArea::ComputePolygonArea( const double *padX,const double* padY,int nCount )
{
	double x1, y1, dx, dy;
	double Qbar1, Qbar2;

	if (NULL == padX || NULL == padY)
	{
		return 0;
	}
	

	if (nCount < 3)
	{
		return 0;
	}
	
	

	double x2 = DEG2RAD( padX[nCount-1] );
	double y2 = DEG2RAD( padY[nCount-1] );
	Qbar2 = GetQbar( y2 );

	double area = 0.0;

	for ( int i = 0; i < nCount; i++ )
	{
		x1 = x2;
		y1 = y2;
		Qbar1 = Qbar2;

		x2 = DEG2RAD( padX[i] );
		y2 = DEG2RAD( padY[i] );
		Qbar2 = GetQbar( y2 );

		if ( x1 > x2 )
		  while ( x1 - x2 > M_PI )
			x2 += m_TwoPI;
		else if ( x2 > x1 )
		  while ( x2 - x1 > M_PI )
			x1 += m_TwoPI;

		dx = x2 - x1;
		area += dx * ( m_Qp - GetQ( y2 ) );

		if (( dy = y2 - y1 ) != 0.0 )
		  area += dx * GetQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
	}
	if (( area *= m_AE ) < 0.0 )
		area = -area;

	if ( area > m_E )
		area = m_E;
	if ( area > m_E / 2 )
		area = m_E - area;

	return area;
}

主な関数はComputePolygonAreaであり,計算された面積単位は平方メートルである.
テスト例は次のとおりです.

std::vector<double> vecX;
	std::vector<double> vecY;
	vecX.push_back(116.0120);
	vecX.push_back(116.0121);
	vecX.push_back(116.01205);

	vecY.push_back(40.004);
	vecY.push_back(40.004);
	vecY.push_back(40.0041);

	DistanceArea area;
	area.SetEllipsePara(6378140.0,6356755.0);
	double aaa = area.ComputePolygonArea(&vecX[0],&vecY[0],vecY.size());

テストを経て、要求を満たすことができます.