Point in polygon test - spherical coords

This is a discussion on Point in polygon test - spherical coords within the C Programming forums, part of the General Programming Boards category; Hello there, I have an algorithm that tests a point for being inside a polygon using the rays algorithm. It ...

  1. #1
    Registered User
    Join Date
    Sep 2007
    Posts
    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_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.

  2. #2
    Frequently Quite Prolix dwks's Avatar
    Join Date
    Apr 2005
    Location
    Canada
    Posts
    8,045
    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.
    dwk

    Seek and ye shall find. quaere et invenies.

    "Simplicity does not precede complexity, but follows it." -- Alan Perlis
    "Testing can only prove the presence of bugs, not their absence." -- Edsger Dijkstra
    "The only real mistake is the one from which we learn nothing." -- John Powell


    Other boards: DaniWeb, TPS
    Unofficial Wiki FAQ: cpwiki.sf.net

    My website: http://dwks.theprogrammingsite.com/
    Projects: codeform, xuni, atlantis, nort, etc.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Why is my program freezing?
    By ShadowMetis in forum Windows Programming
    Replies: 8
    Last Post: 08-20-2004, 03:20 PM
  2. Question About Linker Errors In Dev-C++
    By ShadowMetis in forum C++ Programming
    Replies: 9
    Last Post: 08-18-2004, 08:42 PM
  3. Replies: 5
    Last Post: 10-30-2002, 09:23 PM
  4. point lies in circle?
    By cozman in forum Game Programming
    Replies: 3
    Last Post: 12-20-2001, 03:39 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21