[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:
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(¢roide)!=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(¢roide);
aLine.addPoint(&startPt);
aLine.addPoint(¢roide);
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(¢roide);
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(¢roide);
OGRGeometry* ints=Intersection(&dLine);
radials.addGeometry(ints);
}
More information about the geos-devel
mailing list