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

Quote:

---

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 %= 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.