-
Notifications
You must be signed in to change notification settings - Fork 2
/
isinpolygon.c
110 lines (93 loc) · 2.39 KB
/
isinpolygon.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
/*
* A mex function to check if points are inside a polygon
*
* [inp] = inpolygon(polygon, xy);
*
* polygon : an m x 2 array with polygon vertices
* xy : m x 2 array of points, one point per row
* inp : m x 1 logical vector
*
* Ulf Griesmann, NIST, June 2013
*/
#include <math.h>
#include <string.h>
#include <mex.h>
#ifdef __GNUC__
#define RESTRICT __restrict
#define INLINE __inline__
#else
#define RESTRICT
#define INLINE
#endif
/*-- local prototypes -----------------------------------------*/
INLINE int
in_polygon(int N, double * RESTRICT pol, double X, double Y);
/*-------------------------------------------------------------*/
void
mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
double *poly; /* polygon data */
double *xy; /* array of points */
mxLogical *inp; /* pointer to output data */
int k, Np, Mp, Nxy, Mxy;
/*
* check number of arguments
*/
if (nrhs < 2) {
mexErrMsgTxt("isinpolygon : missing input argument(s).");
}
/*
* check polygon
*/
Np = mxGetN(prhs[0]);
Mp = mxGetM(prhs[0]);
if (Mp < 3 || Np != 2) {
mexErrMsgTxt("isinpolygon : empty input polygon or wrong layout.");
}
poly = mxGetData(prhs[0]);
/*
* check array of points
*/
Nxy = mxGetN(prhs[1]);
Mxy = mxGetM(prhs[1]);
if (Nxy != 2) {
mexErrMsgTxt("isinpolygon : polygon vertex array has wrong layout.");
}
xy = mxGetData(prhs[1]);
/*
* create output array
*/
plhs[0] = mxCreateLogicalMatrix(Mxy, 1);
inp = (mxLogical *)mxGetData(plhs[0]);
/*
* check points
*/
for (k=0; k<Mxy; k++) {
inp[k] = (mxLogical)in_polygon(Mp, poly, xy[k], xy[k+Mxy]);
}
}
/* Return != 0 if a point is inside the polygon
* N vertices make up the polygon pol.
*
* Algorithm from Paul Bourke, http://paulbourke.net/geometry
*/
INLINE int
in_polygon(int N, double * RESTRICT pol, double X, double Y)
{
int i, j, c = 0;
register double * RESTRICT xp, * RESTRICT yp;
/*
* point to x and y values of polygon vertices
* Note: xy is transposed !
*/
xp = pol;
yp = pol + N;
for (i = 0, j = N-1; i < N; j = i++) {
if (( ((yp[i] <= Y) && (Y < yp[j])) ||
((yp[j] <= Y) && (Y < yp[i]))) &&
(X < (xp[j] - xp[i]) * (Y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
c = !c;
}
return c;
}