1 /************************************************************************
2 *
3 *
4 * Test for C-Wrapper of GEOS library
5 *
6 * Copyright (C) 2005 Refractions Research Inc.
7 *
8 * This is free software; you can redistribute and/or modify it under
9 * the terms of the GNU Lesser General Public Licence as published
10 * by the Free Software Foundation.
11 * See the COPYING file for more information.
12 *
13 * Author: Sandro Santilli <strk@kbt.io>
14 *
15 ***********************************************************************/
16
17 #define _GNU_SOURCE
18
19 #include <unistd.h>
20 #include <pthread.h>
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <stdarg.h>
24
25 #include "geos_c.h"
26
27 #define MAXWKTLEN 1047551
28
29 void
usage(char * me)30 usage(char *me)
31 {
32 fprintf(stderr, "Usage: %s <wktfile>\n", me);
33 exit(1);
34 }
35
36 void
notice1(const char * fmt,...)37 notice1(const char *fmt, ...) {
38 va_list ap;
39
40 fprintf( stdout, "NOTICE1: ");
41
42 va_start (ap, fmt);
43 vfprintf( stdout, fmt, ap);
44 va_end(ap);
45 fprintf( stdout, "\n" );
46 }
47
48 void
notice2(const char * fmt,...)49 notice2(const char *fmt, ...) {
50 va_list ap;
51
52 fprintf( stdout, "NOTICE2: ");
53
54 va_start (ap, fmt);
55 vfprintf( stdout, fmt, ap);
56 va_end(ap);
57 fprintf( stdout, "\n" );
58 }
59
60 void
notice(const char * fmt,...)61 notice(const char *fmt, ...) {
62 va_list ap;
63
64 fprintf( stdout, "NOTICE: ");
65
66 va_start (ap, fmt);
67 vfprintf( stdout, fmt, ap);
68 va_end(ap);
69 fprintf( stdout, "\n" );
70 }
71
72 void
log_and_exit(const char * fmt,...)73 log_and_exit(const char *fmt, ...) {
74 va_list ap;
75
76 fprintf( stdout, "ERROR: ");
77
78 va_start (ap, fmt);
79 vfprintf( stdout, fmt, ap);
80 va_end(ap);
81 fprintf( stdout, "\n" );
82 exit(1);
83 }
84
85 void
log_and_exit1(const char * fmt,...)86 log_and_exit1(const char *fmt, ...) {
87 va_list ap;
88
89 fprintf( stdout, "ERROR1: ");
90
91 va_start (ap, fmt);
92 vfprintf( stdout, fmt, ap);
93 va_end(ap);
94 fprintf( stdout, "\n" );
95 pthread_exit(NULL);
96 }
97
98 void
log_and_exit2(const char * fmt,...)99 log_and_exit2(const char *fmt, ...) {
100 va_list ap;
101
102 fprintf( stdout, "ERROR2: ");
103
104 va_start (ap, fmt);
105 vfprintf( stdout, fmt, ap);
106 va_end(ap);
107 fprintf( stdout, "\n" );
108 pthread_exit(NULL);
109 }
110
111 GEOSGeometry*
fineGrainedReconstructionTest(const GEOSGeometry * g1)112 fineGrainedReconstructionTest(const GEOSGeometry* g1)
113 {
114 GEOSCoordSequence* cs;
115 GEOSGeometry* g2;
116 GEOSGeometry* shell;
117 const GEOSGeometry* gtmp;
118 GEOSGeometry**geoms;
119 unsigned int ngeoms, i;
120 int type;
121
122 /* Geometry reconstruction from CoordSeq */
123 type = GEOSGeomTypeId(g1);
124 switch ( type )
125 {
126 case GEOS_POINT:
127 cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(g1));
128 g2 = GEOSGeom_createPoint(cs);
129 return g2;
130 break;
131 case GEOS_LINESTRING:
132 cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(g1));
133 g2 = GEOSGeom_createLineString(cs);
134 return g2;
135 break;
136 case GEOS_LINEARRING:
137 cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(g1));
138 g2 = GEOSGeom_createLinearRing(cs);
139 return g2;
140 break;
141 case GEOS_POLYGON:
142 gtmp = GEOSGetExteriorRing(g1);
143 cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(gtmp));
144 shell = GEOSGeom_createLinearRing(cs);
145 ngeoms = GEOSGetNumInteriorRings(g1);
146 geoms = (GEOSGeometry**)malloc(ngeoms*sizeof(GEOSGeometry*));
147 for (i=0; i<ngeoms; i++)
148 {
149 gtmp = GEOSGetInteriorRingN(g1, i);
150 cs = GEOSCoordSeq_clone(GEOSGeom_getCoordSeq(gtmp));
151 geoms[i] = GEOSGeom_createLinearRing(cs);
152 }
153 g2 = GEOSGeom_createPolygon(shell, geoms, ngeoms);
154 free(geoms);
155 return g2;
156 break;
157 case GEOS_MULTIPOINT:
158 case GEOS_MULTILINESTRING:
159 case GEOS_MULTIPOLYGON:
160 case GEOS_GEOMETRYCOLLECTION:
161 ngeoms = GEOSGetNumGeometries(g1);
162 geoms = (GEOSGeometry**)malloc(ngeoms*sizeof(GEOSGeometry*));
163 for (i=0; i<ngeoms; i++)
164 {
165 gtmp = GEOSGetGeometryN(g1, i);
166 geoms[i] = fineGrainedReconstructionTest(gtmp);
167 }
168 g2 = GEOSGeom_createCollection(type, geoms, ngeoms);
169 free(geoms);
170 return g2;
171 break;
172 default:
173 log_and_exit("Unknown geometry type %d\n", type);
174 return NULL;
175 }
176 }
177
178 void
printHEX(FILE * where,const unsigned char * bytes,size_t n)179 printHEX(FILE *where, const unsigned char *bytes, size_t n)
180 {
181 static char hex[] = "0123456789ABCDEF";
182 int i;
183
184 for (i=0; i<n; i++)
185 {
186 fprintf(where, "%c%c", hex[bytes[i]>>4], hex[bytes[i]&0x0F]);
187 }
188 }
189
190 int
do_all(char * inputfile)191 do_all(char *inputfile)
192 {
193 GEOSGeometry* g1;
194 GEOSGeometry* g2;
195 GEOSGeometry* g3;
196 GEOSGeometry* g4;
197 const GEOSGeometry **gg;
198 unsigned int npoints, ndims;
199 static char wkt[MAXWKTLEN];
200 FILE *input;
201 char *ptr;
202 unsigned char* uptr;
203 size_t size;
204 double dist, area;
205
206 input = fopen(inputfile, "r");
207 if ( ! input ) { perror("fopen"); exit(1); }
208
209 size = fread(wkt, 1, MAXWKTLEN-1, input);
210 fclose(input);
211 if ( ! size ) { perror("fread"); exit(1); }
212 if ( size == MAXWKTLEN-1 ) { perror("WKT input too big!"); exit(1); }
213 wkt[size] = '\0'; /* ensure it is null terminated */
214
215 /* WKT input */
216 g1 = GEOSGeomFromWKT(wkt);
217
218 /* WKT output */
219 ptr = GEOSGeomToWKT(g1);
220 printf("Input (WKT): %s\n", ptr);
221 free(ptr);
222
223 /* WKB output */
224 uptr = GEOSGeomToWKB_buf(g1, &size);
225 printf("Input (WKB): "); printHEX(stdout, uptr, size); putchar('\n');
226
227 /* WKB input */
228 g2 = GEOSGeomFromWKB_buf(uptr, size); free(uptr);
229 if ( ! GEOSEquals(g1, g2) ) log_and_exit("Round WKB conversion failed");
230 GEOSGeom_destroy(g2);
231
232 /* Size and dimension */
233 npoints = GEOSGetNumCoordinates(g1);
234 ndims = GEOSGeom_getDimensions(g1);
235 printf("Geometry coordinates: %dx%d\n", npoints, ndims);
236
237 /* Geometry fine-grained deconstruction/reconstruction test */
238 g2 = fineGrainedReconstructionTest(g1);
239 if ( ! GEOSEquals(g1, g2) )
240 {
241 log_and_exit("Reconstruction test failed\n");
242 }
243 GEOSGeom_destroy(g2);
244
245 /* Unary predicates */
246 if ( GEOSisEmpty(g1) ) printf("isEmpty\n");
247 if ( GEOSisValid(g1) ) printf("isValid\n");
248 if ( GEOSisSimple(g1) ) printf("isSimple\n");
249 if ( GEOSisRing(g1) ) printf("isRing\n");
250
251 /* Convex Hull */
252 g2 = GEOSConvexHull(g1);
253 if ( ! g2 )
254 {
255 log_and_exit("GEOSConvexHull() raised an exception");
256 }
257 ptr = GEOSGeomToWKT(g2);
258 printf("ConvexHull: %s\n", ptr);
259 free(ptr);
260
261 /* Buffer */
262 GEOSGeom_destroy(g1);
263 g1 = GEOSBuffer(g2, 100, 30);
264 if ( ! g1 )
265 {
266 log_and_exit("GEOSBuffer() raised an exception");
267 }
268 ptr = GEOSGeomToWKT(g1);
269 printf("Buffer: %s\n", ptr);
270 free(ptr);
271
272
273 /* Intersection */
274 g3 = GEOSIntersection(g1, g2);
275 if ( ! GEOSEquals(g3, g2) )
276 {
277 GEOSGeom_destroy(g1);
278 GEOSGeom_destroy(g2);
279 GEOSGeom_destroy(g3);
280 log_and_exit("Intersection(g, Buffer(g)) didn't return g");
281 }
282 ptr = GEOSGeomToWKT(g3);
283 printf("Intersection: %s\n", ptr);
284 GEOSGeom_destroy(g3);
285 free(ptr);
286
287 /* Difference */
288 g3 = GEOSDifference(g1, g2);
289 ptr = GEOSGeomToWKT(g3);
290 printf("Difference: %s\n", ptr);
291 GEOSGeom_destroy(g3);
292 free(ptr);
293
294 /* SymDifference */
295 g3 = GEOSSymDifference(g1, g2);
296 ptr = GEOSGeomToWKT(g3);
297 printf("SymDifference: %s\n", ptr);
298 free(ptr);
299
300 /* Boundary */
301 g4 = GEOSBoundary(g3);
302 ptr = GEOSGeomToWKT(g4);
303 printf("Boundary: %s\n", ptr);
304 GEOSGeom_destroy(g3);
305 GEOSGeom_destroy(g4);
306 free(ptr);
307
308 /* Union */
309 g3 = GEOSUnion(g1, g2);
310 if ( ! GEOSEquals(g3, g1) )
311 {
312 GEOSGeom_destroy(g1);
313 GEOSGeom_destroy(g2);
314 GEOSGeom_destroy(g3);
315 log_and_exit("Union(g, Buffer(g)) didn't return Buffer(g)");
316 }
317 ptr = GEOSGeomToWKT(g3);
318 printf("Union: %s\n", ptr);
319 free(ptr);
320
321 /* PointOnSurcace */
322 g4 = GEOSPointOnSurface(g3);
323 ptr = GEOSGeomToWKT(g4);
324 printf("PointOnSurface: %s\n", ptr);
325 GEOSGeom_destroy(g3);
326 GEOSGeom_destroy(g4);
327 free(ptr);
328
329 /* Centroid */
330 g3 = GEOSGetCentroid(g2);
331 ptr = GEOSGeomToWKT(g3);
332 printf("Centroid: %s\n", ptr);
333 GEOSGeom_destroy(g3);
334 free(ptr);
335
336 /* Relate (and RelatePattern )*/
337 ptr = GEOSRelate(g1, g2);
338 if ( ! GEOSRelatePattern(g1, g2, ptr) )
339 {
340 GEOSGeom_destroy(g1);
341 GEOSGeom_destroy(g2);
342 free(ptr);
343 log_and_exit("! RelatePattern(g1, g2, Relate(g1, g2))");
344 }
345 printf("Relate: %s\n", ptr);
346 free(ptr);
347
348 /* Polygonize */
349 gg = (const GEOSGeometry**)malloc(2*sizeof(GEOSGeometry*));
350 gg[0] = g1;
351 gg[1] = g2;
352 g3 = GEOSPolygonize(gg, 2);
353 free(gg);
354 if ( ! g3 )
355 {
356 log_and_exit("Exception running GEOSPolygonize");
357 }
358 ptr = GEOSGeomToWKT(g3);
359 GEOSGeom_destroy(g3);
360 printf("Polygonize: %s\n", ptr);
361 free(ptr);
362
363 /* LineMerge */
364 g3 = GEOSLineMerge(g1);
365 if ( ! g3 )
366 {
367 log_and_exit("Exception running GEOSLineMerge");
368 }
369 ptr = GEOSGeomToWKT(g3);
370 printf("LineMerge: %s\n", ptr);
371 free(ptr);
372 GEOSGeom_destroy(g3);
373
374 /* Binary predicates */
375 if ( GEOSIntersects(g1, g2) ) printf("Intersect\n");
376 if ( GEOSDisjoint(g1, g2) ) printf("Disjoint\n");
377 if ( GEOSTouches(g1, g2) ) printf("Touches\n");
378 if ( GEOSCrosses(g1, g2) ) printf("Crosses\n");
379 if ( GEOSWithin(g1, g2) ) printf("Within\n");
380 if ( GEOSContains(g1, g2) ) printf("Contains\n");
381 if ( GEOSOverlaps(g1, g2) ) printf("Overlaps\n");
382
383 /* Distance */
384 if ( GEOSDistance(g1, g2, &dist) ) printf("Distance: %g\n", dist);
385
386 /* Area */
387 if ( GEOSArea(g1, &area) ) printf("Area 1: %g\n", area);
388 if ( GEOSArea(g2, &area) ) printf("Area 2: %g\n", area);
389
390 GEOSGeom_destroy(g2);
391
392 /* Simplify */
393 g3 = GEOSSimplify(g1, 0.5);
394 ptr = GEOSGeomToWKT(g3);
395 printf("Simplify: %s\n", ptr);
396 free(ptr);
397 GEOSGeom_destroy(g3);
398
399 /* Topology Preserve Simplify */
400 g3 = GEOSTopologyPreserveSimplify(g1, 0.5);
401 ptr = GEOSGeomToWKT(g3);
402 printf("Simplify: %s\n", ptr);
403 free(ptr);
404 GEOSGeom_destroy(g3);
405
406 GEOSGeom_destroy(g1);
407
408 return 0;
409 }
410
threadfunc1(void * arg)411 void *threadfunc1( void *arg )
412 {
413 initGEOS( notice1, log_and_exit1 );
414 printf("GEOS version %s\n", GEOSversion());
415 putc('.', stderr); fflush(stderr);
416 do_all((char*)arg);
417 putc('+', stderr); fflush(stderr);
418 finishGEOS();
419
420 pthread_exit(NULL);
421 }
422
threadfunc2(void * arg)423 void *threadfunc2( void *arg )
424 {
425 initGEOS( notice2, log_and_exit2 );
426 printf("GEOS version %s\n", GEOSversion());
427 putc('.', stderr); fflush(stderr);
428 do_all((char *)arg);
429 putc('+', stderr); fflush(stderr);
430 finishGEOS();
431
432 pthread_exit(NULL);
433 }
434
435 int
main(int argc,char ** argv)436 main(int argc, char **argv)
437 {
438 pthread_t thread1, thread2;
439
440 if ( argc < 2 ) usage(argv[0]);
441 pthread_create( &thread1, NULL, threadfunc1, argv[1] );
442 pthread_create( &thread2, NULL, threadfunc2, argv[1] );
443
444 pthread_join( thread1, NULL );
445 pthread_join( thread2, NULL );
446
447 return EXIT_SUCCESS;
448 }
449
450