[Gdal-dev] Center of polygon

Chapman, Martin MChapman at sanz.com
Mon Apr 5 16:01:40 EDT 2004


Yo,

I know how to do this in c++ code:

// example method call:
vector< float> pCenter = pPolygon->GetCenter(pPolygon->GetXPoints(),
pPolygon->GetYPoints(), pEx);


// the supporting routines
double CSimpleGeometry::CalculateProjectedArea(vector< float>* pXPoints,
vector< float>* pYPoints)
{
	try
	{
		int j = 0;
		double nArea = 0.0;
		long nSize = pXPoints->size();

		if (nSize < 2)
			return 0.0;

		for (int i = 0; i < nSize; ++i) 
		{
			j = (i + 1) % nSize;
			nArea += pXPoints->at(i) * pYPoints->at(j);
			nArea -= pYPoints->at(i) * pXPoints->at(j);
		}

		nArea = nArea / 2.0;
		return (nArea < 0.0 ? -nArea : nArea);
	}
	catch(...)
	{
		throw CString("An unspecified error occurred");
	}
}

vector< float> CSimpleGeometry::GetCenter(vector< float>* pXPoints,
vector< float>* pYPoints, vector< float>* pExtents)
{
	try
	{
		float cx = 0.0f;
		float cy = 0.0f;
		float nArea = (float) CalculateProjectedArea(pXPoints,
pYPoints);
		int j = 0;
		float nFactor = 0.0f;
		vector< float> pCenter;

		long nSize = pXPoints->size();

		if (nSize < 2)
		{
			pCenter.push_back(0.0f);
			pCenter.push_back(0.0f);
			return pCenter;
		}

		if (nSize == 2)
		{
			pCenter.push_back(pXPoints->at(0));
			pCenter.push_back(pXPoints->at(1));
			return pCenter;
		}

		for (int i = 0; i < nSize; i++) 
		{
			j = (i + 1) % nSize;
	      
			nFactor = pXPoints->at(i) * pYPoints->at(j) -
pXPoints->at(j) * pYPoints->at(i);
			cx += (pXPoints->at(i) + pXPoints->at(j)) *
nFactor;
			cy += (pYPoints->at(i) + pYPoints->at(j)) *
nFactor;
		}
	    
		nArea = nArea * 6.0f;
		nFactor = 1 / nArea;
		cx = abs(cx * nFactor);
		cy = abs(cy * nFactor);

		float pXExtents[5] = {pExtents->at(1), pExtents->at(1),
pExtents->at(3), pExtents->at(3), pExtents->at(1)};
		float pYExtents[5] = {pExtents->at(2), pExtents->at(0),
pExtents->at(0), pExtents->at(2), pExtents->at(2)};

		if ( IsPointInPolygon(5, pXExtents, pYExtents, cx, cy))
		{
			pCenter.push_back(cy);
			pCenter.push_back(cx);
		}
		else if ( IsPointInPolygon(5, pXExtents, pYExtents, cx *
-1.0f, cy))
		{
			pCenter.push_back(cy);
			pCenter.push_back(cx * -1);
		}
		else if ( IsPointInPolygon(5, pXExtents, pYExtents, cx,
cy * -1.0f))
		{
			pCenter.push_back(cy * -1);
			pCenter.push_back(cx);
		}
		else if ( IsPointInPolygon(5, pXExtents, pYExtents, cx *
-1.0f, cy * -1.0f))
		{
			pCenter.push_back(cy * -1);
			pCenter.push_back(cx * -1);
		}
		else
		{
			pCenter.push_back(cy);
			pCenter.push_back(cx);
		}

		return pCenter;
	}
	catch(...)
	{
		throw CString("An unspecified error occurred");
	}
}

bool CSimpleGeometry::IsPointInPolygon(int nPolySize, float* pXPoints,
float* pYPoints, float nX, float nY)
{
	try
	{
		long i, j;
		bool bIsInside = false;

		for (i = 0, j = nPolySize - 1; i < nPolySize; j = i++) 
		{
			if ((((pYPoints[i] <= nY) && (nY < pYPoints[j]))
||
			    ((pYPoints[j] <= nY) && (nY < pYPoints[i])))
&&
				(nX < (pXPoints[j] - pXPoints[i]) * 
				(nY - pYPoints[i]) / (pYPoints[j] -
pYPoints[i]) + pXPoints[i]))
			{
				bIsInside = !bIsInside;
			}
		}
		return bIsInside;
	}
	catch(...)
	{
		throw CString("An unspecified error occurred");
	}
}

vector< float> CSimpleGeometry::GetMidPoint(vector< float>* pPoints)
{
	try
	{
		float nY1 = 0.0f, nX1 = 0.0f, nY2 = 0.0f, nX2 = 0.0f,
nDistance = 0.0f;
		int j = 0;
		vector< float> pMidPoint;
		nDistance = (float) GetDistance(pPoints);
		float nMidPoint = nDistance / 2;
		nDistance = 0.0f;
		vector< float>::iterator it;

		for (it = pPoints->begin(); it != pPoints->end(); j++)
		{
			if (j == 0)
			{
				nY1 = *it;
				it++;
				nX1 = *it;
				it++;
			}

			nY2 = *it;
			it++;
			nX2 = *it;
			it++;

			nDistance += sqrt(((nX2 - nX1) * (nX2 - nX1)) +
((nY2 - nY1) * (nY2 - nY1)));

			if (nDistance > nMidPoint)
			{
				nY2 = nY1 + ((nY2 - nY1) / 2);
				nX2 = nX1 + ((nX2 - nX1) / 2);
				break;
			}

			nY1 = nY2;
			nX1 = nX2;
		}

		pMidPoint.push_back(nY2);
		pMidPoint.push_back(nX2);

		return pMidPoint;
	}
	catch(...)
	{
		throw CString("An unspecified error occurred");
	}
}

double CSimpleGeometry::GetDistance(vector< float>* pPoints)
{
	try
	{
		vector< float>::iterator it;
		float nY1 = 0.0f, nX1 = 0.0f, nY2 = 0.0f, nX2 = 0.0f,
nDistance = 0.0f;
		int j = 0;
		long nSize = pPoints->size();

		for (it = pPoints->begin(); it != pPoints->end(); j++)
		{
			if (j == 0)
			{
				nY1 = *it;
				it++;
				nX1 = *it;
				it++;
			}

			nY2 = *it;
			it++;
			nX2 = *it;
			it++;

			nDistance += sqrt(((nX2 - nX1) * (nX2 - nX1)) +
((nY2 - nY1) * (nY2 - nY1)));
			
			nY1 = nY2;
			nX1 = nX2;
		}

		return nDistance;
	}
	catch(...)
	{
		throw CString("An unspecified error occurred");
	}
}

There you go, you can use these functions to find the center point of an
n-gon and much more...

Martin Chapman

-----Original Message-----
From: Frank Warmerdam [mailto:warmerdam at pobox.com] 
Sent: Monday, April 05, 2004 1:35 PM
To: gdal-dev at remotesensing.org
Subject: Re: [Gdal-dev] Center of polygon


Clay, Bruce wrote:
> Is there a method within OGR to find the center of a polygon (where it

> is not part of the data set) or polyline for the purposes of centering

> a label?  If not does anyone have a method they are willing to share 
> to walk through the data of a polygon / ring to determine the center?

Bruce,

There is no support this currently in GDAL/OGR.  You might find
something like this in the general GEOS geometry engine.  There is also
some centroid code in Shapelibs contrib directory if you would like to
adapt that.

Best regards,
-- 
---------------------------------------+--------------------------------
---------------------------------------+------
I set the clouds in motion - turn up   | Frank Warmerdam,
warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent

_______________________________________________
Gdal-dev mailing list
Gdal-dev at remotesensing.org
http://remotesensing.org/mailman/listinfo/gdal-dev



More information about the Gdal-dev mailing list