[geos-devel] Computing Radials inside a Polygon

[Jo]

First of all, I apologize to post this question when my code is not
using GEOS but I think this is a generic geometry problem so it really
does not matter which API, I am using (I'm using OGR).
I need to compute a certain number of equally spaced radials inside a
polygon, given that each segment should start on the centroide and end
somewhere in its boundary (or external ring).
The algorithm I designed, follows this steps:
- given a certain number of radials(n), calculate the angle(fangle)
that equally divides the circle in n/2 sectors;
- fetch a point from the perimeter of the minimum enclosing circle (to
guarantee the point is outside the perimeter);
- loop n times, rotating this point around using fangle for this
rotation; at each step, the rotated vector will be a vector that links
the centroide to the rotated point; to get a vector inside the
polygon, this vector is intersected with the polygon and we use the
resulting geometry.

I am getting a strange (wrong) result, so I guess I am not seeing
things right (or there is a bug in my code :-)) because I get equally
spaced radials covering approximately a quarter if the polygon, and
the angle from the first point to the first rotated point does not
seem correct... image here:

http://ladybug.no-ip.org/files/radials.png


Here is my code! I would be very grateful for any suggestions that
could point me in the right direction... also, if you have any ideas
of computing this using GEOS I would not mind to go for them;
                               thanks in advance,
                                                           cheers, Jo

	OGRPoint centroide;					
	if (Centroid(&centroide)!=OGRERR_NONE) return false; // get centroid of polygon
	const OGRLinearRing * pRing=minCircle.getExteriorRing(); //get
external ring of the minimum enclosing circle

	if (pRing->getNumPoints()<1) {printf("Ring not valid!\n");return false;}
	OGRPoint startPt;
	pRing->getPoint(0, &startPt);//get one point in the minimum enclosing
circle external ring

	OGRLineString aLine;
	aLine.addPoint(&centroide);
	aLine.addPoint(&startPt);
	aLine.addPoint(&centroide);

	if (!aLine.IsValid()) {printf("Line is not valid!\n");return false;}

	OGRGeometry* ints=Intersection(&aLine);
	radials.addGeometry(ints);

	double ns=n/0.5; //n is the number of radials
	double fangle=PI/(ns*0.5);// assumed that the whole circumference is
2 PI radians

	for (size_t i=0; i < n; ++i){

		OGRLineString dLine;
		dLine.addPoint(&centroide);

		double rx=startPt.getX()*cos(fangle)-startPt.getY()*sin(fangle);
//rotated x- is this correct?
		double ry=startPt.getX()*sin(fangle)+startPt.getY()*cos(fangle);
//rotated y- is this correct?

		startPt.setX(rx);
		startPt.setY(ry);
		dLine.addPoint(&startPt);
		dLine.addPoint(&centroide);

		OGRGeometry* ints=Intersection(&dLine);
		radials.addGeometry(ints);

	}