[geos-devel] Re: Computing Radials inside a Polygon

Jo doublebyte at gmail.com
Sat Aug 29 17:53:41 EDT 2009


Thinking better, I could see what I was doing wrong.
I was rotating a point around (0,0), rather than rotating it around my
origin (the centroid). Therefore, I had to add the difference between
the point and centroid coordinates to the sin and cosine terms:
	
		double rx=centroid->getX()+cos(fangle)*(x-centroid->getX())-sin(fangle)*(y-centroid->getY());
		double ry=centroid->getY()+sin(fangle)*(x-centroid->getX())+cos(fangle)*(y-centroid->getY());

And calling consecutively this function, gives me the equally spaced
radials from a centroid of a polygon to its external ring; in this
example, using a 45 degrees angle:

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

                                   cheers,
                                              Jo


2009/8/29  <geos-devel-request at lists.osgeo.org>:
> Send geos-devel mailing list submissions to
>        geos-devel at lists.osgeo.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
>        http://lists.osgeo.org/mailman/listinfo/geos-devel
> or, via email, send a message with subject or body 'help' to
>        geos-devel-request at lists.osgeo.org
>
> You can reach the person managing the list at
>        geos-devel-owner at lists.osgeo.org
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of geos-devel digest..."
>
>
> Today's Topics:
>
>   1. Re: Computing Radials inside a Polygon (Jo)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Sat, 29 Aug 2009 16:16:23 +0100
> From: Jo <doublebyte at gmail.com>
> Subject: [geos-devel] Re: Computing Radials inside a Polygon
> To: geos-devel at lists.osgeo.org
> Message-ID:
>        <23ab5f0a0908290816j66928c44y58de481b26525657 at mail.gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1
>
> 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.
>
> Hi Sanak,
> First of all, thanks for your reply!
> I think I didn't explain myself clearly though...
>
> The idea is to compute radials from the centroid of the polygon, to
> its external ring, at equal intervals. Something like this:
>
> http://ladybug.no-ip.org/files/radials_final.png
>
> I had the idea of starting with a vector that starts in the centroid
> and ends up in the external ring (at any point), and therefore its
> length equals the radius; to make things simple I choose a point that
> is in the same x-alignment as the centroid, and so the start angle in
> your code would be null. My start vector is something like this:
>
> http://ladybug.no-ip.org/files/start.png
>
> Than my idea, was to rotate this vector around the centroid, at a
> constant angle, and get the radials.
> Lets assume for simplicity that the angle is 90 degrees (PI/2 radians)
> and therefore I want to end up with four radials.
> As I compute the new locations of the rotated endpoint of the vector,
> I store this point on a stl container, so that in the end I can loop
> through all these locations and generate my vectors.
> Here is a GEOS version of my code:
>
>        //////////GEOS VERSION /////
>        std::vector<geos::geom::Point*> pts;
>
>        geos::geom::GeometryFactory factory;
>        geos::geom::Coordinate coordC(ogrCentroid.getX(),ogrCentroid.getY());
>
>        //create centroid from given point
>        geos::geom::Point* centroid = factory.createPoint(coordC);
>
>        //create start pt as any point on the perimeter (for instance, a
> point aligned in x with the centroid)
>        geos::geom::Coordinate coordS(ogrCentroid.getX()+radius,ogrCentroid.getY());
>        geos::geom::Point* startPt = factory.createPoint(coordS);
>
>        pts.push_back(startPt);
>
>        //nr radials=4; angle=90 deg= PI/2
>        double fangle=PI/2.0;
>
>        for (size_t i=0; i <4; ++i){
>
>                double x=pts.back()->getX();
>                double y=pts.back()->getY();
>
>                double rx=x*cos(fangle)-y*sin(fangle);
>                double ry=x*sin(fangle)+y*cos(fangle);
>
>                geos::geom::Coordinate coords(x+rx,y+ry);
>                geos::geom::Point* aPoint = factory.createPoint(coords);
>
>                pts.push_back(aPoint);
>        }
>
> Let me mention that this code is *not* doing what I want, since:
> - I am getting the wrong angle between the start point and the first
> iteration point;
> - the angles do not appear to be 90 degrees;
> - I think also the angles are not constant;
>
> Here is a screenshot of my output:
>
> http://ladybug.no-ip.org/files/output.png
>
> So I guess there must be some assumption (or more than one, actually)
> that is wrong in my code :-)
> Again, I would be very grateful for any observations or suggestions...
> thanks in advance for your help!
>
>                                       cheers,
>                                                     Jo
>
>
> ------------------------------
>
> _______________________________________________
> geos-devel mailing list
> geos-devel at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/geos-devel
>
> End of geos-devel Digest, Vol 82, Issue 25
> ******************************************
>



-- 
"#define QUESTION ((bb) || !(bb))"  (Shakespeare)


More information about the geos-devel mailing list