# Thread: Point in polygon test - spherical coords

1. ## Point in polygon test - spherical coords

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_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.

2. 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.
Here, I presume:
Code:
```			if(xi > X_MAX)
xi-=X_RANGE;
else if(xi < X_MIN)
xi+=X_RANGE;```
That works reasonably well, unless xi is, say, over X_RANGE outside of the valid range. A fix for this would be to use while loops instead of if statements. However, the modulus operator will take care of the positive side for you. This
Code:
`while(xi > X_MAX) xi -= X_RANGE;`
(should that be ">="?) is the same as
Code:
`xi &#37;= X_RANGE;`
Of course, modulus doesn't work for floating point numbers -- for those, you need to use fmod(): http://www.cprogramming.com/fod/fmod.html

You should also put that code into a function to avoid repeating it.