# [geos-devel] Computing Radials inside a Polygon

Jo doublebyte at gmail.com
Sun Aug 23 13:18:37 EDT 2009

```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:

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;
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;

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

OGRGeometry* ints=Intersection(&aLine);

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

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

OGRLineString dLine;

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);