Hi Jo,<br><br>I think that each rotate angles should be calculated more simply.<br>(I think rotation matrix is valid when centroid is equal to origine(0, 0).)<br><br>The following code is pseudo code. (Sorry I haven't OGR environment...)<br>
<br>===================================<br>// get minCircle radius. (If radius is not calculated yet)<br>double radius=Distance(&centroid, &startPt); // I don't know that Distance function is implemented in OGR...<br>
// get start angle.<br>double startAngle=atan2(&startPt.getY()-&centroid.getY(), &startPt.getX()-&centroid.getX());<br><br>for (size_t i=0; i <n; ++i){<br> OGRLineString dLine;<br> dLine.addPoint(&centroid);<br>
<br> double rx=centroid.getX() + radius*cos(startAngle+fangle*i);<br> double ry=centroid.getY() + radius*sin(startAngle+fangle*i);<br> :<br>}<br>===================================<br><br>Regards,<br>
<br>Sanak.<br><br><div class="gmail_quote">2009/8/24 Jo <span dir="ltr"><<a href="mailto:doublebyte@gmail.com">doublebyte@gmail.com</a>></span><br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
First of all, I apologize to post this question when my code is not<br>
using GEOS but I think this is a generic geometry problem so it really<br>
does not matter which API, I am using (I'm using OGR).<br>
I need to compute a certain number of equally spaced radials inside a<br>
polygon, given that each segment should start on the centroide and end<br>
somewhere in its boundary (or external ring).<br>
The algorithm I designed, follows this steps:<br>
- given a certain number of radials(n), calculate the angle(fangle)<br>
that equally divides the circle in n/2 sectors;<br>
- fetch a point from the perimeter of the minimum enclosing circle (to<br>
guarantee the point is outside the perimeter);<br>
- loop n times, rotating this point around using fangle for this<br>
rotation; at each step, the rotated vector will be a vector that links<br>
the centroide to the rotated point; to get a vector inside the<br>
polygon, this vector is intersected with the polygon and we use the<br>
resulting geometry.<br>
<br>
I am getting a strange (wrong) result, so I guess I am not seeing<br>
things right (or there is a bug in my code :-)) because I get equally<br>
spaced radials covering approximately a quarter if the polygon, and<br>
the angle from the first point to the first rotated point does not<br>
seem correct... image here:<br>
<br>
<a href="http://ladybug.no-ip.org/files/radials.png" target="_blank">http://ladybug.no-ip.org/files/radials.png</a><br>
<br>
<br>
Here is my code! I would be very grateful for any suggestions that<br>
could point me in the right direction... also, if you have any ideas<br>
of computing this using GEOS I would not mind to go for them;<br>
thanks in advance,<br>
cheers, Jo<br>
<br>
OGRPoint centroide;<br>
if (Centroid(&centroide)!=OGRERR_NONE) return false; // get centroid of polygon<br>
const OGRLinearRing * pRing=minCircle.getExteriorRing(); //get<br>
external ring of the minimum enclosing circle<br>
<br>
if (pRing->getNumPoints()<1) {printf("Ring not valid!\n");return false;}<br>
OGRPoint startPt;<br>
pRing->getPoint(0, &startPt);//get one point in the minimum enclosing<br>
circle external ring<br>
<br>
OGRLineString aLine;<br>
aLine.addPoint(&centroide);<br>
aLine.addPoint(&startPt);<br>
aLine.addPoint(&centroide);<br>
<br>
if (!aLine.IsValid()) {printf("Line is not valid!\n");return false;}<br>
<br>
OGRGeometry* ints=Intersection(&aLine);<br>
radials.addGeometry(ints);<br>
<br>
double ns=n/0.5; //n is the number of radials<br>
double fangle=PI/(ns*0.5);// assumed that the whole circumference is<br>
2 PI radians<br>
<br>
for (size_t i=0; i < n; ++i){<br>
<br>
OGRLineString dLine;<br>
dLine.addPoint(&centroide);<br>
<br>
double rx=startPt.getX()*cos(fangle)-startPt.getY()*sin(fangle);<br>
//rotated x- is this correct?<br>
double ry=startPt.getX()*sin(fangle)+startPt.getY()*cos(fangle);<br>
//rotated y- is this correct?<br>
<br>
startPt.setX(rx);<br>
startPt.setY(ry);<br>
dLine.addPoint(&startPt);<br>
dLine.addPoint(&centroide);<br>
<br>
OGRGeometry* ints=Intersection(&dLine);<br>
radials.addGeometry(ints);<br>
<br>
}<br>
_______________________________________________<br>
geos-devel mailing list<br>
<a href="mailto:geos-devel@lists.osgeo.org">geos-devel@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/geos-devel" target="_blank">http://lists.osgeo.org/mailman/listinfo/geos-devel</a><br>
</blockquote></div><br>