-
Notifications
You must be signed in to change notification settings - Fork 2
/
iscwmex.c
114 lines (95 loc) · 2.62 KB
/
iscwmex.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
111
112
113
114
/*
* A mex function to compute the orientation of polygons
*
* [cw] = iscw(pa);
*
* pa : cell array of polygons (nx2 matrices)
* cw : 1xn logical vector of orientations
*
* Ulf Griesmann, NIST, November 2012
*/
#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
clock_wise(double * RESTRICT p, int N);
/*-------------------------------------------------------------*/
void
mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
mxArray *par; /* pointer to array structure */
double *pda; /* polynom data */
mxLogical *pout; /* pointer to output data */
int ma;
int k, Na;
/* check argument pa */
if ( !mxIsCell(prhs[0]) ) {
mexErrMsgTxt("iscwmex : argument must be a cell array.");
}
Na = mxGetNumberOfElements(prhs[0]);
if (!Na) {
mexErrMsgTxt("iscwmex : no input polygons pa.");
}
/* create output array */
plhs[0] = mxCreateLogicalMatrix(1, Na);
pout = (mxLogical *)mxGetData(plhs[0]);
/*
* copy and prepare data
*/
/* pa */
for (k=0; k<Na; k++) {
/* get the next polygon from the cell array */
par = mxGetCell(prhs[0], k); /* ptr to mxArray */
pda = mxGetData(par); /* ptr to a data */
ma = mxGetM(par); /* rows = vertex number */
/* calculate orientation */
pout[k] = (mxLogical)clock_wise(pda, ma);
}
}
/* Return the clockwise status of a curve, clockwise or counterclockwise
* N vertices making up curve p.
* It is assumed that
* - the polygon is closed
* - the last point is not repeated.
* - the polygon is simple (does not intersect itself or have holes)
*
* Algorithm from Paul Bourke, http://paulbourke.net/geometry/polygonmesh/
*/
INLINE int
clock_wise(double * RESTRICT p, int N)
{
int i,j,k,n;
int count = 0;
register double z;
/* check if too few vertices */
if (N < 3)
mexErrMsgTxt("iscwmex : polygons must have at lest 3 vertices.");
/* check if the last vertex repeats the first */
n = N;
if ((p[0] == p[N-1]) && (p[N] == p[2*N-1]))
n = N-1;
/* calculate cross products */
for (i=0;i<n;i++) {
j = (i + 1) % n;
k = (i + 2) % n;
z = (p[j] - p[i]) * (p[N+k] - p[N+j]) -
(p[k] - p[j]) * (p[N+j] - p[N+i]);
if (z < 0)
count--;
else if (z > 0)
count++;
}
if (count > 0)
return 1;
else
return 0;
}