point in polygon

[Valdas Rapsevicius]

Hello,

in our company we needed a function identify either a point is in the polygon or not. Here is the udf function that does the trick:

//#####################################################
#ifdef STANDARD
#include <stdio.h>
#include <string.h>
#ifdef __WIN__
typedef unsigned __int64 ulonglong; /* Microsofts 64 bit types */
typedef __int64 longlong;
#else
typedef unsigned long long ulonglong;
typedef long long longlong;
#endif /*__WIN__*/
#else
#include <my_global.h>
#include <my_sys.h>
#endif
#include <mysql.h>
#include <m_ctype.h>
#include <m_string.h> // To get strmov()

extern "C" {

my_bool pointinpoly_init(UDF_INIT *initid, UDF_ARGS *args, char *message);
void pointinpoly_deinit(UDF_INIT *initid);
long long pointinpoly(UDF_INIT *initid, UDF_ARGS *args, char *is_null, char *error);

}

my_bool pointinpoly_init(UDF_INIT *initid, UDF_ARGS *args, char *message){
if (args->arg_count != 2 || args->arg_type[0] != STRING_RESULT || args->arg_type[1] != STRING_RESULT){
strmov(message ,"Usage: pointinpoly(point_wkb, poly_wkb)");
return 1;
}
return 0;
}

void pointinpoly_deinit(UDF_INIT *initid){ }

long long pointinpoly(UDF_INIT *initid, UDF_ARGS *args, char *is_null, char *error){
char* wkbp = args->args[0];
char* wkb = args->args[1];
double x;
double y;
double xi, yi, xj, yj;
int status = 0;

uint n_linear_rings;
uint n_points;
uint wkb_type;

wkbp += 10;
wkb_type = uint4korr(wkbp);
if(wkb_type != 1){
return 0;
}
wkbp += 4;
float8get(x, wkbp);
wkbp += 8;
float8get(y, wkbp);


// wkb geometry collection: byte_order
++wkb;
// wkb geometry collection: wkbType //wkb_type = uint4korr(wkb);
wkb += 4;
// wkb geometry collection: num_wkbGeometries //wkb_geom_num = uint4korr(wkb);
wkb += 4;
// wkb geometry: byte_order
++wkb;
// wkb geometry: wkbType
wkb_type = uint4korr(wkb);
if(wkb_type != 3){
return 0;
}
wkb += 4;

n_linear_rings = uint4korr(wkb);
wkb += 4;
for (; n_linear_rings > 0; --n_linear_rings){
n_points = uint4korr(wkb);
wkb += 4;

// Imame paskutinio tasko koordinate
wkb += 2*8*(n_points-1);
float8get(xj, wkb);
wkb += 8;
float8get(yj, wkb);
wkb -= 2*8*(n_points-1) + 8;

for(int i = 0; i < n_points; i++){
float8get(xi, wkb);
wkb += 8;
float8get(yi, wkb);
wkb += 8;
if((((yi <= y) && (y < yj)) || ((yj <= y) && (y < yi))) && (x < (xj - xi) * (y - yi) / (yj - yi) + xi)){
status = !status;
}
xj = xi;
yj = yi;
}
}
return status;
}
//#####################################################