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

[Jo]

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