src/Rgpc.c

   1 /*
   2   gpclib:  General Polygon Clipping library for R
   3   Copyright (C) 2003 Roger D. Peng <rpeng@jhsph.edu>
   4 */
   5
   6 #include <R.h>
   7 #include <Rinternals.h>
   8 #include "gpc.h"
   9
  10 static int compute_polygon_size(gpc_polygon *p);
  11 static void gpc_polygon_to_double(double *a, int na, gpc_polygon *p);
  12 static void double_to_gpc_polygon(gpc_polygon *p, double *a, int na);
  13
  14
  15
  16 SEXP Rgpc_polygon_clip(SEXP subjpoly, SEXP clippoly, SEXP op) {
  17     gpc_polygon subject, clip, result;
  18     int polysize, nsubj, nclip, iop;
  19     SEXP returnval;
  20     double *xsubjpoly, *xclippoly, *xreturnval;
  21
  22     PROTECT(subjpoly = coerceVector(subjpoly, REALSXP));
  23     PROTECT(clippoly = coerceVector(clippoly, REALSXP));
  24     PROTECT(op = coerceVector(op, INTSXP));
  25
  26     nsubj = length(subjpoly);
  27     nclip = length(clippoly);
  28     xsubjpoly = REAL(subjpoly);
  29     xclippoly = REAL(clippoly);
  30     iop = INTEGER(op)[0];
  31
  32     double_to_gpc_polygon(&subject, xsubjpoly, nsubj);
  33     double_to_gpc_polygon(&clip, xclippoly, nclip);
  34
  35     if(iop == 1)
  36         gpc_polygon_clip(GPC_INT, &subject, &clip, &result);
  37     else if(iop == 2)
  38         gpc_polygon_clip(GPC_DIFF, &subject, &clip, &result);
  39     else
  40         gpc_polygon_clip(GPC_UNION, &subject, &clip, &result);
  41
  42     polysize = compute_polygon_size(&result);
  43
  44     PROTECT(returnval = allocVector(REALSXP, polysize));
  45     xreturnval = REAL(returnval);
  46
  47     gpc_polygon_to_double(xreturnval, polysize, &result);
  48
  49     /*
  50       gpc_free_polygon(&subject);
  51       gpc_free_polygon(&clip);
  52     */
  53     gpc_free_polygon(&result);
  54
  55
  56     UNPROTECT(4);
  57
  58     return(returnval);
  59
  60 }
  61
  62
  63 /* unserialize the polygon */
  64
  65 static void double_to_gpc_polygon(gpc_polygon *p, double *a, int na)
  66 {
  67     int i, j, k;
  68
  69     p->num_contours = a[0];
  70     p->hole = (int *)R_alloc(p->num_contours, sizeof(int));
  71     p->contour = (gpc_vertex_list *)R_alloc(p->num_contours, sizeof(gpc_vertex_list));
  72     i = 1;
  73
  74     for(j=0; j < p->num_contours; j++) {
  75         p->contour[j].num_vertices = a[i++];
  76         p->contour[j].vertex = (gpc_vertex *)R_alloc(p->contour[j].num_vertices,
  77                                                      sizeof(gpc_vertex));
  78         p->hole[j] = (int) a[i++];
  79
  80         for(k=0; k < p->contour[j].num_vertices; k++) {
  81             p->contour[j].vertex[k].x = a[i++];
  82             p->contour[j].vertex[k].y = a[i++];
  83         }
  84         if(i > na) {
  85             Rprintf("index out of range: %d\n", i);
  86             return;
  87         }
  88     }
  89 }
  90
  91 /* serialize polygon to vector */
  92
  93 static void gpc_polygon_to_double(double *a, int na, gpc_polygon *p)
  94 {
  95     int i, j, k;
  96
  97     a[0] = p->num_contours;
  98     i = 1;
  99
 100     for(j=0; j < p->num_contours; j++) {
 101         a[i++] = p->contour[j].num_vertices;
 102         a[i++] = (double) p->hole[j];
 103
 104         if(i > na) {
 105             Rprintf("index out of range: %d\n", i);
 106             return;
 107         }
 108         for(k=0; k < p->contour[j].num_vertices; k++) {
 109             a[i++] = p->contour[j].vertex[k].x;
 110
 111             if(i > na) {
 112                 Rprintf("index out of range: %d\n", i);
 113                 return;
 114             }
 115             a[i++] = p->contour[j].vertex[k].y;
 116
 117             if(i > na) {
 118                 Rprintf("index out of range: %d\n", i);
 119                 return;
 120             }
 121         }
 122     }
 123 }
 124
 125
 126 static int compute_polygon_size(gpc_polygon *p)
 127 {
 128     int psize = 1, i;
 129
 130     psize += p->num_contours;
 131     psize += p->num_contours; /* for the hole flags */
 132
 133     for(i=0; i < p->num_contours; i++) {
 134         psize += 2 * p->contour[i].num_vertices;
 135     }
 136     return(psize);
 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    Older code had separate functions for intersect/union/diff.  These
 168    are now done with a single function + flag (duh!).  But I'll save these
 169    functions just in case....
 170 */
 171
 172 /*********************************************************************
 173
 174 SEXP gpc_polygon_intersect(SEXP subjpoly, SEXP clippoly)
 175 {
 176     gpc_polygon subject, clip, result;
 177     int polysize, nsubj, nclip;
 178     SEXP returnval;
 179     double *xreturnval;
 180     double *xsubjpoly, *xclippoly;
 181
 182     PROTECT(subjpoly = AS_NUMERIC(subjpoly));
 183     PROTECT(clippoly = AS_NUMERIC(clippoly));
 184     nsubj = LENGTH(subjpoly);
 185     nclip = LENGTH(clippoly);
 186
 187     xsubjpoly = NUMERIC_POINTER(subjpoly);
 188     xclippoly = NUMERIC_POINTER(clippoly);
 189
 190     double_to_gpc_polygon(&subject, xsubjpoly, nsubj);
 191     double_to_gpc_polygon(&clip, xclippoly, nclip);
 192     gpc_polygon_clip(GPC_INT, &subject, &clip, &result);
 193
 194     polysize = compute_polygon_size(&result);
 195
 196     PROTECT(returnval = NEW_NUMERIC(polysize));
 197     xreturnval = NUMERIC_POINTER(returnval);
 198
 199     gpc_polygon_to_double(xreturnval, polysize, &result);
 200
 201     gpc_free_polygon(&subject);
 202     gpc_free_polygon(&clip);
 203     gpc_free_polygon(&result);
 204
 205     UNPROTECT(3);
 206
 207     return(returnval);
 208 }
 209
 210
 211
 212 SEXP gpc_polygon_difference(SEXP subjpoly, SEXP clippoly)
 213 {
 214     gpc_polygon subject, clip, result;
 215     int polysize, nsubj, nclip;
 216     SEXP returnval;
 217     double *xreturnval;
 218     double *xsubjpoly, *xclippoly;
 219
 220     PROTECT(subjpoly = AS_NUMERIC(subjpoly));
 221     PROTECT(clippoly = AS_NUMERIC(clippoly));
 222     nsubj = LENGTH(subjpoly);
 223     nclip = LENGTH(clippoly);
 224
 225     xsubjpoly = NUMERIC_POINTER(subjpoly);
 226     xclippoly = NUMERIC_POINTER(clippoly);
 227
 228     double_to_gpc_polygon(&subject, xsubjpoly, nsubj);
 229     double_to_gpc_polygon(&clip, xclippoly, nclip);
 230     gpc_polygon_clip(GPC_DIFF, &subject, &clip, &result);
 231
 232     polysize = compute_polygon_size(&result);
 233
 234     PROTECT(returnval = NEW_NUMERIC(polysize));
 235     xreturnval = NUMERIC_POINTER(returnval);
 236
 237     gpc_polygon_to_double(xreturnval, polysize, &result);
 238
 239     gpc_free_polygon(&subject);
 240     gpc_free_polygon(&clip);
 241     gpc_free_polygon(&result);
 242
 243     UNPROTECT(3);
 244
 245     return(returnval);
 246 }
 247
 248
 249 SEXP gpc_polygon_union(SEXP subjpoly, SEXP clippoly)
 250 {
 251     gpc_polygon subject, clip, result;
 252     int polysize, nsubj, nclip;
 253     SEXP returnval;
 254     double *xreturnval;
 255     double *xsubjpoly, *xclippoly;
 256
 257     PROTECT(subjpoly = AS_NUMERIC(subjpoly));
 258     PROTECT(clippoly = AS_NUMERIC(clippoly));
 259     nsubj = LENGTH(subjpoly);
 260     nclip = LENGTH(clippoly);
 261
 262     xsubjpoly = NUMERIC_POINTER(subjpoly);
 263     xclippoly = NUMERIC_POINTER(clippoly);
 264
 265     double_to_gpc_polygon(&subject, xsubjpoly, nsubj);
 266     double_to_gpc_polygon(&clip, xclippoly, nclip);
 267     gpc_polygon_clip(GPC_UNION, &subject, &clip, &result);
 268
 269     polysize = compute_polygon_size(&result);
 270
 271     PROTECT(returnval = NEW_NUMERIC(polysize));
 272     xreturnval = NUMERIC_POINTER(returnval);
 273
 274     gpc_polygon_to_double(xreturnval, polysize, &result);
 275
 276     gpc_free_polygon(&subject);
 277     gpc_free_polygon(&clip);
 278     gpc_free_polygon(&result);
 279
 280     UNPROTECT(3);
 281
 282     return(returnval);
 283 }
 284 ********************************************************************/