Hello there,
I have an algorithm that tests a point for being inside a polygon using the rays algorithm.
It works fine and i was able to make a small optimisation but i find it hard to adapt the code to work with spherical coordinates: x,y are degrees and x axis is actually a circle and x belongs to [-180*10^5; 180*10^5) range.
What i did to the algorithm is to substract 360*10^5 degrees if the value is more than 180*10^5, and add 360*10^5 if value is less than -180*10^5.
I tested the code with several polygons but some points fail.
The integer coordinates are with 5 grades of accuracy
and x in [-18000000; 18000000)
y in [-9000000; 9000000]
Code:
/*polygon_point_test.h*/
#ifndef _POLYGON_POINT_TEST_H_
#define _POLYGON_POINT_TEST_H_
typedef int IntCoord_t;
#define X_MAX 18000000 // the opposite of greenwich, west
#define X_MIN -18000000 // the opposite of greenwich, east
#define Y_MAX 9000000 // north pole
#define Y_MIN -9000000 // south pole
#define X_RANGE 36000000
typedef struct PointInt2d
{
IntCoord_t x, y;
}PointInt2d_t;
#endif
this is the modified function: is_inside_polygon_rays_int3
Code:
/*polygon_point_test.c*/
#include "polygon_point_test.h"
int is_inside_polygon_rays_int3(PointInt2d_t *poly, size_t n, PointInt2d_t testPt)
{
IntCoord_t dx, dy, a, xi, xj, xmin, xmax;
int oddNodes;
size_t i, j=n-1;
oddNodes=0;
for(i=0;i<n;i++)
{
if((poly[i].y < testPt.y && poly[j].y >= testPt.y) ||
(poly[j].y < testPt.y && poly[i].y >= testPt.y))
{
dx = poly[j].x - poly[i].x;
xi = testPt.x - poly[i].x;
xj = testPt.x - poly[j].x;
xmax = testPt.x-X_MAX;
xmin = testPt.x-X_MIN;
//xmin = X_MIN
if(xi > X_MAX)
xi-=X_RANGE;
else if(xi < X_MIN)
xi+=X_RANGE;
if(xj > X_MAX)
xj-=X_RANGE;
else if(xj < X_MIN)
xj+=X_RANGE;
if(dx > X_MAX)
dx-=X_RANGE;
else if(dx < X_MIN)
dx+=X_RANGE;
if(xmin < X_MIN)
xmin+=X_RANGE;
else if(xmin > X_MAX)
xmin-=X_RANGE;
if(xmax < X_MIN)
xmax+=X_RANGE;
else if(xmax > X_MAX)
xmax-=X_RANGE;
printf("xmin=%d, xmax=%d, xi=%d, xj=%d\n",xmin,xmax,xi,xj);
if(0 < xi && 0 < xj )//&& xi < xmax && xj < xmax)
{
// the ray is crossing the side
oddNodes=1-oddNodes;
}
else if(!(0 > xi && 0 > xj ))//&& xi > xmin && xj > xmin))
{
dy = poly[j].y - poly[i].y;
a = testPt.y - poly[i].y;
if(dy > 0) {
if (a*dx < xi*dy) // the ray is crossing the side
oddNodes=1-oddNodes;
}
else if (a*dx > xi*dy)
oddNodes=1-oddNodes;
}
}
j=i;
}
return oddNodes;
}
this is the working version in normal flat coordinates:
Code:
/*...*/
/*polygon_point_test.c*/
int is_inside_polygon_rays_int21(PointInt2d_t *poly, size_t n, PointInt2d_t testPt)
{
int oddNodes;
size_t i, j;
IntCoord_t dy, dxi, dxj;
oddNodes=0;
j=n-1;
for(i=0;i<n;i++)
{
if((poly[i].y < testPt.y && poly[j].y >= testPt.y) ||
(poly[j].y < testPt.y && poly[i].y >= testPt.y))
{
dxi = testPt.x - poly[i].x;
dxj = testPt.x - poly[j].x;
if(0 < dxi && 0 < dxj)
oddNodes=1-oddNodes;
else if(!(0 > dxi && 0 > dxj))
{
dy = poly[j].y - poly[i].y;
if(dy > 0)
{
if ((testPt.y-poly[i].y)*(poly[j].x-poly[i].x) < dxi*dy)
oddNodes=1-oddNodes;
}
else if((testPt.y-poly[i].y)*(poly[j].x-poly[i].x)> dxi*dy)
oddNodes=1-oddNodes;
}
}
j=i;
}
return oddNodes;
}
for convinience here is the testing function and main.c
Code:
/**
This tests the extreeme polygons that is positioned on the opposite side of the greenwich meridian
*/
int test_rays2_extremities()
{
// int calc;
//odd polygon
PointInt2d_t poly[] = {
{ 16000000, -1000000},
{ 16000000, 1000000},
{-16000000, 1000000},
{-16000000, -1000000},
};
PointInt2d_t testPt ={ 17000000, 0}; //inside the polygon
PointInt2d_t testPt6={-18000000, 0}; //inside
PointInt2d_t testPt5={-17000000, 0}; //inside the polygon
PointInt2d_t testPt2={-15000000, 0}; //left of the polygon, => rays2_int shouldnt calculate
PointInt2d_t testPt3={ 15000000, 0}; //right of the polygon, => rays2__int should calculate twice
PointInt2d_t testPt4={ 0, 0}; //opposite outside
/*
//normal polygon
PointInt2d_t poly[] = {
{ 2000000, -1000000},
{ 2000000, 1000000},
{ -2000000, 1000000},
{ -2000000, -1000000},
};
PointInt2d_t testPt ={ 1, 0}; //inside the polygon
PointInt2d_t testPt6={ 0, 0}; //inside
PointInt2d_t testPt5={ -1, 0}; //inside the polygon
PointInt2d_t testPt2={ -2500000, 0}; //left of the polygon, => rays2_int shouldnt calculate
PointInt2d_t testPt3={ 2500000, 0}; //right of the polygon, => rays2__int should calculate twice
PointInt2d_t testPt4={ 17900000, 0}; //opposite outside
*/
int result, ret;
ret=0;
//expect in -> expect result=1
printf("\ntesting center...\n");
result = is_inside_polygon_rays_int3(poly, ARRAY_COUNT(poly), testPt);
printf("point(%d,%d) is ",testPt.x,testPt.y);
if(result == 1)
printf("inside! -> SUCCESS!\n");
else if(result == 0)
printf("outside! -> FAILED!\n");
//
printf("\ntesting center...\n");
result = is_inside_polygon_rays_int3(poly, ARRAY_COUNT(poly), testPt5);
printf("point(%d,%d) is ",testPt5.x,testPt5.y);
if(result == 1)
printf("inside! -> SUCCESS!\n");
else if(result == 0)
printf("outside! -> FAILED!\n");
//
printf("\ntesting center...\n");
result = is_inside_polygon_rays_int3(poly, ARRAY_COUNT(poly), testPt6);
printf("point(%d,%d) is ",testPt6.x,testPt6.y);
if(result == 1)
printf("inside! -> SUCCESS!\n");
else if(result == 0)
printf("outside! -> FAILED!\n");
//left
printf("\ntesting left...\n");
result = is_inside_polygon_rays_int3(poly, ARRAY_COUNT(poly), testPt2);
printf("point(%d,%d) is ",testPt2.x,testPt2.y);
if(result == 1)
printf("inside! -> FAILED!\n");
else if(result == 0)
printf("outside! -> SUCCESS!\n");
printf("\ntesting right...\n");
result = is_inside_polygon_rays_int3(poly, ARRAY_COUNT(poly), testPt3);
printf("point(%d,%d) is ",testPt3.x,testPt3.y);
if(result == 1)
printf("inside! -> FAILED!\n");
else if(result == 0)
printf("outside! -> SUCCESS!\n");
printf("\ntesting opposite center...\n");
result = is_inside_polygon_rays_int3(poly, ARRAY_COUNT(poly), testPt4);
printf("point(%d,%d) is ",testPt4.x,testPt4.y);
if(result == 1)
printf("inside! -> FAILED!\n");
else if(result == 0)
printf("outside! -> SUCCESS!\n");
return ret;
}
int main(int argc, char *argv[])
{
//test_differencies();
test_rays2_extremities();
//test_rays_from_header_data();
//test_rays_from_header_data2();
//test_rays2_int();
//test_rays2_int_02();
//test_speed();
//test_accuracy();
//printf("%g\n", 10.0/2.0*2.0);
return 0;
}
Any kind of help will be appreciated.