1
2 /*
3 * Modified for R by Nicholas Lewin-Koh also made modifications for
4 * multipart polygons. Sept 29, 2000
5 */
6
7
8 #include <R.h>
9 #include <Rinternals.h>
10 #include <Rmath.h>
11 #include "maptools.h"
12
13 SEXP RshpCentrd_2d (SEXP);
14 SEXP R_RingCentrd_2d (int , SEXP, double *);
15
16
17 /* **************************************************************************
18 * RshpCentrd_2d
19 *
20 * Return the single mathematical / geometric centroid of a potentially
21 * complex/compound RShapeObject
22 *
23 * reject non area SHP Types
24 *
25 * **************************************************************************/
RshpCentrd_2d(SEXP call)26 SEXP RshpCentrd_2d (SEXP call) {
27 int ring, ringPrev, ring_nVertices, rStart, nprts;
28 int i,j,totvert;
29 double Area, ringArea;
30 SEXP ringCentrd, Cent, shape, flag, ringVerts;
31 shape = CADR(call);
32 flag = CADDR(call);
33
34 /* if ( !(SHPDimension(INTEGER(getAttrib(shape,install("shp.type")))[0]) */
35 /* & SHPD_AREA) ) */
36 /* error("Not a class of shape with defined 2d area"); */
37
38 nprts = INTEGER(getAttrib(shape, install("nParts")))[0];
39 Area = 0;
40 if(INTEGER(flag)[0]==0 ||nprts==1){
41 PROTECT(Cent=allocVector(REALSXP, 2));
42 REAL(Cent)[0] = 0.0;
43 REAL(Cent)[1] = 0.0;
44 }
45 else{
46 PROTECT(Cent=allocMatrix(REALSXP, nprts, 2));
47 }
48 /* for each ring in compound / complex object calc the ring cntrd */
49
50 ringPrev = INTEGER(getAttrib(shape, install("nVerts")))[0];
51 totvert = INTEGER(getAttrib(shape, install("nVerts")))[0];
52
53 if(nprts==0) nprts=1;
54 for ( ring = nprts-1; ring >= 0; ring-- ) {
55 rStart = INTEGER(VECTOR_ELT(shape,0))[ring];
56 ring_nVertices = ringPrev - rStart;
57 /* Rprintf("ringPrev= %d, rStart=%d, ring_nVertices=%d \n", */
58 /* ringPrev, rStart, ring_nVertices); */
59
60 PROTECT(ringVerts=allocMatrix(REALSXP, ring_nVertices, 2));
61 for(i=rStart,j=0;i<ringPrev ;i++,j++){
62 REAL(ringVerts)[j]=REAL(VECTOR_ELT(shape,1))[i];
63 REAL(ringVerts)[j+ring_nVertices]=REAL(VECTOR_ELT(shape,1))[i+totvert];
64 }
65 /* Rprintf(" matrix begin %f, matrix end: %f \n", */
66 /* REAL(ringVerts)[0],REAL(ringVerts)[(2*ring_nVertices)-1]); */
67 PROTECT(ringCentrd =
68 R_RingCentrd_2d (ring_nVertices, ringVerts, &ringArea));
69 /* Rprintf("xcent: %f, ycent: %f, area: %f\n ", */
70 /* REAL(ringCentrd)[0],REAL(ringCentrd)[1],ringArea ); */
71
72 /* use Superposition of these rings to build a composite Centroid */
73 /* sum the ring centrds * ringAreas, at the end divide by total area */
74 if(INTEGER(flag)[0]==0 ||nprts==1){
75 REAL(Cent)[0] += REAL(ringCentrd)[0] * ringArea;
76 REAL(Cent)[1] += REAL(ringCentrd)[1] * ringArea;
77 }
78 else{
79 REAL(Cent)[ring]= REAL(ringCentrd)[0];
80 REAL(Cent)[ring+nprts]= REAL(ringCentrd)[1];
81 }
82 Area += ringArea;
83 ringPrev = rStart;
84 UNPROTECT(2);
85 }
86
87 /* hold on the division by AREA until were at the end */
88 if(INTEGER(flag)[0]==0 ||nprts==1){
89 REAL(Cent)[0] = REAL(Cent)[0] / Area;
90 REAL(Cent)[1] = REAL(Cent)[1] / Area;
91 UNPROTECT(1);
92 return ( Cent );
93 }
94 else{
95 UNPROTECT(1);
96 return ( Cent );
97 }
98 }
99
100
101 /* **************************************************************************
102 * RingCentroid_2d
103 * Copyright (c) 1999, Carl Anderson
104 *
105 * This code is based in part on the earlier work of Frank Warmerdam
106 *
107 *
108 * Return the mathematical / geometric centroid of a single closed ring
109 *
110 * **************************************************************************/
R_RingCentrd_2d(int nVert,SEXP xy,double * Area)111 SEXP R_RingCentrd_2d (int nVert, SEXP xy, double *Area ) {
112 int iv /*, jv */;
113 /* int sign_x, sign_y; */
114 double /* dy_Area, */ dx_Area, Cx_accum, Cy_accum, ppx, ppy;
115 double x_base, y_base, x, y;
116 SEXP RingCent;
117 /* the centroid of a closed Ring is defined as
118 *
119 * Cx = sum (cx * dArea ) / Total Area
120 * and
121 * Cy = sum (cy * dArea ) / Total Area
122 */
123
124 x_base = REAL(xy)[0];
125 y_base = REAL(xy)[nVert];
126
127 Cy_accum = 0.0;
128 Cx_accum = 0.0;
129
130 ppx = REAL(xy)[1] - x_base;
131 ppy = REAL(xy)[nVert + 1] - y_base;
132 *Area = 0;
133
134 /* Skip the closing vector */
135 for ( iv = 2; iv <= nVert - 2; iv++ ) {
136 x = REAL(xy)[iv] - x_base;
137 y = REAL(xy)[nVert + iv] - y_base;
138
139 /* calc the area and centroid of triangle built out of an arbitrary */
140 /* base_point on the ring and each successive pair on the ring */
141
142 /* Area of a triangle is the cross product of its defining vectors */
143 /* Centroid of a triangle is the average of its vertices */
144
145 dx_Area = ((x * ppy) - (y * ppx)) * 0.5;
146 *Area += dx_Area;
147
148 Cx_accum += ( ppx + x ) * dx_Area;
149 Cy_accum += ( ppy + y ) * dx_Area;
150 /* #ifdef DEBUG2 */
151 /* printf("(ringcentrd_2d) Pp( %f, %f), P(%f, %f)\n", ppx, ppy, x, y); */
152 /* printf("(ringcentrd_2d) dA: %f, sA: %f, Cx: %f, Cy: %f \n", */
153 /* dx_Area, *Area, Cx_accum, Cy_accum); */
154 /* #endif */
155 ppx = x;
156 ppy = y;
157 }
158
159 /* #ifdef DEBUG2 */
160 /* printf("(ringcentrd_2d) Cx: %f, Cy: %f \n", */
161 /* ( Cx_accum / ( *Area * 3) ), ( Cy_accum / (*Area * 3) )); */
162 /* #endif */
163
164 /* adjust back to world coords
165 */
166 PROTECT(RingCent=allocVector(REALSXP,2));
167 REAL(RingCent)[0] = ( Cx_accum / ( *Area * 3)) + x_base;
168 REAL(RingCent)[1] = ( Cy_accum / ( *Area * 3)) + y_base;
169 UNPROTECT(1);
170 return (RingCent);
171 }
172
173
174 /* Modified 24 May 2005 Roger S. Bivand for maptools
175 Written by Joseph O'Rourke
176 orourke@cs.smith.edu
177 October 27, 1995
178
179 Computes the centroid (center of gravity) of an arbitrary
180 simple polygon via a weighted sum of signed triangle areas,
181 weighted by the centroid of each triangle.
182 Reads x,y coordinates from stdin.
183 NB: Assumes points are entered in ccw order!
184 E.g., input for square:
185 0 0
186 10 0
187 10 10
188 0 10
189 This solves Exercise 12, p.47, of my text,
190 Computational Geometry in C. See the book for an explanation
191 of why this works. Follow links from
192 http://cs.smith.edu/~orourke/
193
194 */
195
196 #define DIM 2 /* Dimension of points */
197 typedef double tPointd[DIM]; /* type double point */
198
199 /*#define PMAX 1000 Max # of pts in polygon */
200 /* typedef tPointd *tPolygond; */ /* type double polygon */
201
202 double Area2( tPointd a, tPointd b, tPointd c );
203 void FindCG( int n, tPointd *P, tPointd CG, double *Areasum2 );
204 void Centroid3( tPointd p1, tPointd p2, tPointd p3, tPointd c );
205
RFindCG(int * n,double * x,double * y,double * xc,double * yc,double * area)206 void RFindCG( int *n, double *x, double *y, double *xc, double *yc,
207 double *area ) {
208
209 int i, nn;
210 tPointd *P;
211 tPointd CG;
212 double Areasum2;
213 nn = n[0];
214 P = (tPointd *) R_alloc((size_t) nn, sizeof(tPointd));
215 for (i=0; i<nn; i++) {
216 P[i][0] = x[i];
217 P[i][1] = y[i];
218 }
219 FindCG(nn, P, CG, &Areasum2);
220 xc[0] = CG[0];
221 yc[0] = CG[1];
222 area[0] = Areasum2/2;
223 return;
224 }
225
226 /*
227 Returns the cg in CG. Computes the weighted sum of
228 each triangle's area times its centroid. Twice area
229 and three times centroid is used to avoid division
230 until the last moment.
231 */
FindCG(int n,tPointd * P,tPointd CG,double * Areasum2)232 void FindCG( int n, tPointd *P, tPointd CG, double *Areasum2)
233 {
234 int i;
235 double A2; /* Partial area sum */
236 tPointd Cent3;
237
238 CG[0] = 0;
239 CG[1] = 0;
240 Areasum2[0] = 0;
241 for (i = 1; i < n-1; i++) {
242 Centroid3( P[0], P[i], P[i+1], Cent3 );
243 A2 = Area2( P[0], P[i], P[i+1]);
244 CG[0] += A2 * Cent3[0];
245 CG[1] += A2 * Cent3[1];
246 Areasum2[0] += A2;
247 }
248 CG[0] /= 3 * Areasum2[0];
249 CG[1] /= 3 * Areasum2[0];
250 return;
251 }
252 /*
253 Returns three times the centroid. The factor of 3 is
254 left in to permit division to be avoided until later.
255 */
Centroid3(tPointd p1,tPointd p2,tPointd p3,tPointd c)256 void Centroid3( tPointd p1, tPointd p2, tPointd p3, tPointd c )
257 {
258 c[0] = p1[0] + p2[0] + p3[0];
259 c[1] = p1[1] + p2[1] + p3[1];
260 return;
261 }
262 /*
263 Returns twice the signed area of the triangle determined by a,b,c,
264 positive if a,b,c are oriented ccw, and negative if cw.
265 */
Area2(tPointd a,tPointd b,tPointd c)266 double Area2( tPointd a, tPointd b, tPointd c )
267 {
268 double area;
269 area = (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]);
270 return(area);
271 }
272
273