/
chull.c
426 lines (390 loc) · 10.4 KB
/
chull.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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
/*
* chull finds the convex hull of a set of points in the plane.
*
* It is based on a C translation (by f2c) of
* ACM TOMS algorithm 523 by W. F. Eddy, vol 3 (1977), 398-403, 411-2.
*
* converted to double precision, output order altered
* by B.D. Ripley, March 1999
*
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "R_ext/Boolean.h"/* TRUE,... */
#include "R_ext/Applic.h"
static void split(int n, double *x,
int m, int *in,
int ii, int jj,
int s,
int *iabv, int *na, int *maxa,
int *ibel, int *nb, int *maxb)
{
/* split() takes the m points of array x whose
subscripts are in array in and partitions them by the
line joining the two points in array x whose subscripts are ii and jj.
The subscripts of the points above the line are put into array
iabv, and the subscripts of the points below are put into array ibel.
na and nb are, respectively, the number of points
above the line and the number below.
maxa and maxb are the subscripts for array
x of the point furthest above the line and the point
furthest below, respectively. if either subset is null
the corresponding subscript (maxa or maxb) is set to zero.
formal parameters
INPUT
n integer total number of data points
x real array (2,n) (x,y) co-ordinates of the data
m integer number of points in input subset
in integer array (m) subscripts for array x of the
points in the input subset
ii integer subscript for array x of one point
on the partitioning line
jj integer subscript for array x of another
point on the partitioning line
s integer switch to determine output.
refer to comments below
OUTPUT
iabv integer array (m) subscripts for array x of the
points above the partitioning line
na integer number of elements in iabv
maxa integer subscript for array x of point
furthest above the line.
set to zero if na is zero
ibel integer array (m) subscripts for array x of the
points below the partitioning line
nb integer number of elements in ibel
maxb integer subscript for array x of point
furthest below the line.
set to zero if nb is zero
if s = 2 dont save ibel,nb,maxb.
if s =-2 dont save iabv,na,maxa.
otherwise save everything
if s is positive the array being partitioned is above
the initial partitioning line.
if it is negative, then the set of points is below.
*/
/* Local variables (=0 : -Wall) */
double a=0, b=0, down, d1, up, xt, z;
int i, is;
Rboolean vert, neg_dir=0;
/* Parameter adjustments */
--x;
xt = x[ii];
/* Check to see if the line is vertical */
vert = (x[jj] == xt);
d1 = x[jj + n] - x[ii + n];
if (vert) {
neg_dir = ((s > 0 && d1 < 0.) || (s < 0 && d1 > 0.));
} else {
a = d1 / (x[jj] - xt);
b = x[ii + n] - a * xt;
}
up = 0.; *na = 0; *maxa = 0;
down = 0.; *nb = 0; *maxb = 0;
for (i = 0; i < m; ++i) {
is = in[i];
if (vert) {
if(neg_dir) z = xt - x[is];
else z = x[is] - xt;
} else {
z = x[is + n] - a * x[is] - b;
}
if (z > 0.) { /* the point is ABOVE the line */
if (s == -2) continue;
iabv[*na] = is;
++(*na);
if (z >= up) {
up = z;
*maxa = *na;
}
}
else if (s != 2 && z < 0.) { /* the point is BELOW the line */
ibel[*nb] = is;
++(*nb);
if (z <= down) {
down = z;
*maxb = *nb;
}
}
}
}
void R_chull(int *n, double *x, int *m, int *in,
int *ia, int *ib,
int *ih, int *nh, int *il)
{
/* this subroutine determines which of the m points of array
x whose subscripts are in array in are vertices of the
minimum area convex polygon containing the m points. the
subscripts of the vertices are placed in array ih in the
order they are found. nh is the number of elements in
array ih and array il. array il is a linked list giving
the order of the elements of array ih in a counter
clockwise direction. this algorithm corresponds to a
preorder traversal of a certain binary tree. each vertex
of the binary tree represents a subset of the m points.
at each step the subset of points corresponding to the
current vertex of the tree is partitioned by a line
joining two vertices of the convex polygon. the left son
vertex in the binary tree represents the subset of points
above the partitioning line and the right son vertex, the
subset below the line. the leaves of the tree represent
either null subsets or subsets inside a triangle whose
vertices coincide with vertices of the convex polygon.
formal parameters
INPUT
n integer total number of data points (= nrow(x))
x real array (2,n) (x,y) co-ordinates of the data
m integer number of points in the input subset
in integer array (m) subscripts for array x of the points
in the input subset
work area
ia integer array (m) subscripts for array x of left son subsets.
see comments after dimension statements
ib integer array (m) subscripts for array x of right son subsets
OUTPUT
ih integer array (m) subscripts for array x of the
vertices of the convex hull
nh integer number of elements in arrays ih and il.
== number of vertices of the convex polygon
il is used internally here.
il integer array (m) a linked list giving in order in a
counter-clockwise direction the
elements of array ih
the upper end of array ia is used to store temporarily
the sizes of the subsets which correspond to right son
vertices, while traversing down the left sons when on the
left half of the tree, and to store the sizes of the left
sons while traversing the right sons(down the right half)
*/
#define y(k) x[k + x_dim1]
Rboolean mine, maxe;
int i, j, ilinh, ma, mb, kn, mm, kx, mx, mp1, mbb, nia, nib,
inh, min, mxa, mxb, mxbb;
int x_dim1, x_offset;
double d1;
/* Parameter adjustments */
x_dim1 = *n;
x_offset = 1;
x -= x_offset;
--il;
--ih;
--ib;
--ia;
--in;
if (*m == 1) {
goto L_1pt;
}
il[1] = 2;
il[2] = 1;
kn = in[1];
kx = in[2];
if (*m == 2) {
goto L_2pts;
}
mp1 = *m + 1;
min = 1;
mx = 1;
kx = in[1];
maxe = FALSE;
mine = FALSE;
/* find two vertices of the convex hull for the initial partition */
for (i = 2; i <= *m; ++i) {
j = in[i];
if ((d1 = x[j] - x[kx]) < 0.) {
} else if (d1 == 0) {
maxe = TRUE;
} else {
maxe = FALSE;
mx = i;
kx = j;
}
if ((d1 = x[j] - x[kn]) < 0.) {
mine = FALSE;
min = i;
kn = j;
} else if (d1 == 0) {
mine = TRUE;
}
}
if (kx == kn) { /* if the max and min are equal,
* all m points lie on a vertical line */
goto L_vertical;
}
if (maxe || mine) {/* if maxe (or mine) is TRUE, there are several
maxima (or minima) with equal first coordinates */
if (maxe) {/* have several points with the (same) largest x[] */
for (i = 1; i <= *m; ++i) {
j = in[i];
if (x[j] != x[kx]) continue;
if (y(j) <= y(kx)) continue;
mx = i;
kx = j;
}
}
if (mine) {/* have several points with the (same) smallest x[] */
for (i = 1; i <= *m; ++i) {
j = in[i];
if (x[j] != x[kn]) continue;
if (y(j) >= y(kn)) continue;
min = i;
kn = j;
}
}
}
/* L7:*/
ih[1] = kx;
ih[2] = kn;
*nh = 3;
inh = 1;
nib = 1;
ma = *m;
in[mx] = in[*m];
in[*m] = kx;
mm = *m - 2;
if (min == *m) {
min = mx;
}
in[min] = in[*m - 1];
in[*m - 1] = kn;
/* begin by partitioning the root of the tree */
split(*n, &x[x_offset], mm, &in[1],
ih[1], ih[2],
0,
&ia[1], &mb, &mxa,
&ib[1], &ia[ma], &mxbb);
/* first traverse the LEFT HALF of the tree */
/* start with the left son */
L8:
nib += ia[ma];
--ma;
do {
if (mxa != 0) {
il[*nh] = il[inh];
il[inh] = *nh;
ih[*nh] = ia[mxa];
ia[mxa] = ia[mb];
--mb;
++(*nh);
if (mb != 0) {
ilinh = il[inh];
split(*n, &x[x_offset], mb, &ia[1],
ih[inh], ih[ilinh],
1,
&ia[1], &mbb, &mxa,
&ib[nib], &ia[ma], &mxb);
mb = mbb;
goto L8;
}
/* then the right son */
inh = il[inh];
}
do {
inh = il[inh];
++ma;
nib -= ia[ma];
if (ma >= *m) goto L12;
} while(ia[ma] == 0);
ilinh = il[inh];
/* on the left side of the tree, the right son of a right son */
/* must represent a subset of points which is inside a */
/* triangle with vertices which are also vertices of the */
/* convex polygon and hence the subset may be neglected. */
split(*n, &x[x_offset], ia[ma], &ib[nib],
ih[inh], ih[ilinh],
2,
&ia[1], &mb, &mxa,
&ib[nib], &mbb, &mxb);
ia[ma] = mbb;
} while(TRUE);
/* now traverse the RIGHT HALF of the tree */
L12:
mxb = mxbb;
ma = *m;
mb = ia[ma];
nia = 1;
ia[ma] = 0;
/* start with the right son */
L13:
nia += ia[ma];
--ma;
do {
if (mxb != 0) {
il[*nh] = il[inh];
il[inh] = *nh;
ih[*nh] = ib[mxb];
ib[mxb] = ib[mb];
--mb;
++(*nh);
if (mb != 0) {
ilinh = il[inh];
split(*n, &x[x_offset], mb, &ib[nib],
ih[inh], ih[ilinh],
-1,
&ia[nia], &ia[ma], &mxa,
&ib[nib], &mbb, &mxb);
mb = mbb;
goto L13;
}
/* then the left son */
inh = il[inh];
}
do {
inh = il[inh];
++ma;
nia -= ia[ma];
if (ma == mp1) goto Finis;
} while(ia[ma] == 0);
ilinh = il[inh];
/* on the right side of the tree, the left son of a left son */
/* must represent a subset of points which is inside a */
/* triangle with vertices which are also vertices of the */
/* convex polygon and hence the subset may be neglected. */
split(*n, &x[x_offset], ia[ma], &ia[nia],
ih[inh], ih[ilinh],
-2,
&ia[nia], &mbb, &mxa,
&ib[nib], &mb, &mxb);
} while(TRUE);
/* -------------------------------------------------------------- */
L_vertical:/* all the points lie on a vertical line */
kx = in[1];
kn = in[1];
for (i = 1; i <= *m; ++i) {
j = in[i];
if (y(j) > y(kx)) {
mx = i;
kx = j;
}
if (y(j) < y(kn)) {
min = i;
kn = j;
}
}
if (kx == kn) goto L_1pt;
L_2pts:/* only two points */
ih[1] = kx;
ih[2] = kn;
if (x[kn] == x[kx] && y(kn) == y(kx))
*nh = 2;
else
*nh = 3;
goto Finis;
L_1pt:/* only one point */
*nh = 2;
ih[1] = in[1];
il[1] = 1;
Finis:
--(*nh);
/* put the results in order, as given by IH */
for (i = 1; i <= *nh; ++i) {
ia[i] = ih[i];
}
j = il[1];
for (i = 2; i <= *nh; ++i) {
ih[i] = ia[j];
j = il[j];
}
return;
#undef y
} /* chull */