Computing Radials inside a Polygon
|
|
Jo-4
|
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); } _______________________________________________ geos-devel mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/geos-devel |
||||||||||||||||
sanak
|
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(¢roid, &startPt); // I don't know that Distance function is implemented in OGR... // get start angle. double startAngle=atan2(&startPt.getY()-¢roid.getY(), &startPt.getX()-¢roid.getX()); for (size_t i=0; i <n; ++i){ OGRLineString dLine; dLine.addPoint(¢roid); double rx=centroid.getX() + radius*cos(startAngle+fangle*i); double ry=centroid.getY() + radius*sin(startAngle+fangle*i); : } =================================== Regards, Sanak. 2009/8/24 Jo <[hidden email]> First of all, I apologize to post this question when my code is not _______________________________________________ geos-devel mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/geos-devel |
||||||||||||||||
|
Jo-4
|
In reply to this post
by Jo-4
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(¢roid, &startPt); // I don't know that Distance function is implemented in OGR... // get start angle. double startAngle=atan2(&startPt.getY()-¢roid.getY(), &startPt.getX()-¢roid.getX()); for (size_t i=0; i <n; ++i){ OGRLineString dLine; dLine.addPoint(¢roid); 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 [hidden email] http://lists.osgeo.org/mailman/listinfo/geos-devel |
||||||||||||||||
|
Jo-4
|
In reply to this post
by Jo-4
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 <[hidden email]>: > Send geos-devel mailing list submissions to > [hidden email] > > 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 > [hidden email] > > You can reach the person managing the list at > [hidden email] > > 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 <[hidden email]> > Subject: [geos-devel] Re: Computing Radials inside a Polygon > To: [hidden email] > Message-ID: > <[hidden email]> > 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(¢roid, &startPt); // I don't know that Distance > function is implemented in OGR... > // get start angle. > double startAngle=atan2(&startPt.getY()-¢roid.getY(), > &startPt.getX()-¢roid.getX()); > > for (size_t i=0; i <n; ++i){ > OGRLineString dLine; > dLine.addPoint(¢roid); > > 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 > [hidden email] > http://lists.osgeo.org/mailman/listinfo/geos-devel > > End of geos-devel Digest, Vol 82, Issue 25 > ****************************************** > -- "#define QUESTION ((bb) || !(bb))" (Shakespeare) _______________________________________________ geos-devel mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/geos-devel |
||||||||||||||||
|
sanak
|
Hi Jo,
Thanks for your reply. Ok, your last image seems to be no problem. Regards, Sanak. 2009/8/30 Jo <[hidden email]> Thinking better, I could see what I was doing wrong. _______________________________________________ geos-devel mailing list [hidden email] http://lists.osgeo.org/mailman/listinfo/geos-devel |
||||||||||||||||
| Free Embeddable Forum Powered by Nabble | Help |
