[geos-devel] Computing Radials inside a Polygon

[Sanak Goe]

Hi Jo,

I think that each rotate angles should be calculated more simply.
(I think rotation matrix is valid when centroid is equal to origine(0, 0).)

The following code is pseudo code. (Sorry I haven't OGR environment...)

===================================
// get minCircle radius. (If radius is not calculated yet)
double radius=Distance(&centroid, &startPt); // I don't know that Distance
function is implemented in OGR...
// get start angle.
double startAngle=atan2(&startPt.getY()-&centroid.getY(),
&startPt.getX()-&centroid.getX());

for (size_t i=0; i <n; ++i){
        OGRLineString dLine;
        dLine.addPoint(&centroid);

        double rx=centroid.getX() + radius*cos(startAngle+fangle*i);
        double ry=centroid.getY() + radius*sin(startAngle+fangle*i);
         :
}
===================================

Regards,

Sanak.

2009/8/24 Jo <doublebyte at gmail.com>

> 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);
>
>        }
> _______________________________________________
> geos-devel mailing list
> geos-devel at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/geos-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/geos-devel/attachments/20090824/554629f1/attachment.html