PostGIS  2.3.7dev-r@@SVN_REVISION@@
cu_algorithm.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  * Copyright 2008 Paul Ramsey
6  *
7  * This is free software; you can redistribute and/or modify it under
8  * the terms of the GNU General Public Licence. See the COPYING file.
9  *
10  **********************************************************************/
11 
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include "CUnit/Basic.h"
16 
17 #include "liblwgeom_internal.h"
18 #include "cu_tester.h"
19 
20 /*
21 ** Global variables used by tests below
22 */
23 
24 /* Two-point objects */
25 POINTARRAY *pa21 = NULL;
26 POINTARRAY *pa22 = NULL;
27 LWLINE *l21 = NULL;
28 LWLINE *l22 = NULL;
29 /* Parsing support */
31 
32 
33 /*
34 ** The suite initialization function.
35 ** Create any re-used objects.
36 */
37 static int init_cg_suite(void)
38 {
39  pa21 = ptarray_construct(0, 0, 2);
40  pa22 = ptarray_construct(0, 0, 2);
41  l21 = lwline_construct(SRID_UNKNOWN, NULL, pa21);
42  l22 = lwline_construct(SRID_UNKNOWN, NULL, pa22);
43  return 0;
44 
45 }
46 
47 /*
48 ** The suite cleanup function.
49 ** Frees any global objects.
50 */
51 static int clean_cg_suite(void)
52 {
53  if ( l21 ) lwline_free(l21);
54  if ( l22 ) lwline_free(l22);
55  return 0;
56 }
57 
58 /*
59 ** Test left/right side.
60 */
61 static void test_lw_segment_side(void)
62 {
63  int rv = 0;
64  POINT2D p1, p2, q;
65 
66  /* Vertical line at x=0 */
67  p1.x = 0.0;
68  p1.y = 0.0;
69  p2.x = 0.0;
70  p2.y = 1.0;
71 
72  /* On the left */
73  q.x = -2.0;
74  q.y = 1.5;
75  rv = lw_segment_side(&p1, &p2, &q);
76  //printf("left %g\n",rv);
77  CU_ASSERT(rv < 0);
78 
79  /* On the right */
80  q.x = 2.0;
81  rv = lw_segment_side(&p1, &p2, &q);
82  //printf("right %g\n",rv);
83  CU_ASSERT(rv > 0);
84 
85  /* On the line */
86  q.x = 0.0;
87  rv = lw_segment_side(&p1, &p2, &q);
88  //printf("on line %g\n",rv);
89  CU_ASSERT_EQUAL(rv, 0);
90 
91 }
92 
93 static void test_lw_arc_center(void)
94 {
95 /* double lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result); */
96  POINT2D c1;
97  double d1;
98  POINT2D p1, p2, p3;
99 
100  p1.x = 2047538.600;
101  p1.y = 7268770.435;
102  p2.x = 2047538.598;
103  p2.y = 7268770.435;
104  p3.x = 2047538.596;
105  p3.y = 7268770.436;
106 
107  d1 = lw_arc_center(&p1, &p2, &p3, &c1);
108 
109  CU_ASSERT_DOUBLE_EQUAL(d1, 0.0046097720751, 0.0001);
110  CU_ASSERT_DOUBLE_EQUAL(c1.x, 2047538.599, 0.001);
111  CU_ASSERT_DOUBLE_EQUAL(c1.y, 7268770.4395, 0.001);
112 
113  // printf("lw_arc_center: (%12.12g, %12.12g) %12.12g\n", c1.x, c1.y, d1);
114 }
115 
116 /*
117 ** Test crossings side.
118 */
119 static void test_lw_segment_intersects(void)
120 {
121 
122 #define setpoint(p, x1, y1) {(p).x = (x1); (p).y = (y1);}
123 
124  POINT2D p1, p2, q1, q2;
125 
126  /* P: Vertical line at x=0 */
127  setpoint(p1, 0.0, 0.0);
128  p1.x = 0.0;
129  p1.y = 0.0;
130  p2.x = 0.0;
131  p2.y = 1.0;
132 
133  /* Q: Horizontal line crossing left to right */
134  q1.x = -0.5;
135  q1.y = 0.5;
136  q2.x = 0.5;
137  q2.y = 0.5;
138  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_RIGHT );
139 
140  /* Q: Horizontal line crossing right to left */
141  q1.x = 0.5;
142  q1.y = 0.5;
143  q2.x = -0.5;
144  q2.y = 0.5;
145  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_LEFT );
146 
147  /* Q: Horizontal line not crossing right to left */
148  q1.x = 0.5;
149  q1.y = 1.5;
150  q2.x = -0.5;
151  q2.y = 1.5;
152  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
153 
154  /* Q: Horizontal line crossing at second vertex right to left */
155  q1.x = 0.5;
156  q1.y = 1.0;
157  q2.x = -0.5;
158  q2.y = 1.0;
159  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
160 
161  /* Q: Horizontal line crossing at first vertex right to left */
162  q1.x = 0.5;
163  q1.y = 0.0;
164  q2.x = -0.5;
165  q2.y = 0.0;
166  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_LEFT );
167 
168  /* Q: Diagonal line with large range crossing at first vertex right to left */
169  q1.x = 0.5;
170  q1.y = 10.0;
171  q2.x = -0.5;
172  q2.y = -10.0;
173  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_LEFT );
174 
175  /* Q: Diagonal line with large range crossing at second vertex right to left */
176  q1.x = 0.5;
177  q1.y = 11.0;
178  q2.x = -0.5;
179  q2.y = -9.0;
180  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
181 
182  /* Q: Horizontal touching from left at second vertex*/
183  q1.x = -0.5;
184  q1.y = 0.5;
185  q2.x = 0.0;
186  q2.y = 0.5;
187  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
188 
189  /* Q: Horizontal touching from right at first vertex */
190  q1.x = 0.0;
191  q1.y = 0.5;
192  q2.x = 0.5;
193  q2.y = 0.5;
194  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_RIGHT );
195 
196  /* Q: Horizontal touching from left and far below on second vertex */
197  q1.x = -0.5;
198  q1.y = -10.5;
199  q2.x = 0.0;
200  q2.y = 0.5;
201  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
202 
203  /* Q: Horizontal touching from right and far above on second vertex */
204  q1.x = 0.5;
205  q1.y = 10.5;
206  q2.x = 0.0;
207  q2.y = 0.5;
208  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
209 
210  /* Q: Co-linear from top */
211  q1.x = 0.0;
212  q1.y = 10.0;
213  q2.x = 0.0;
214  q2.y = 0.5;
215  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_COLINEAR );
216 
217  /* Q: Co-linear from bottom */
218  q1.x = 0.0;
219  q1.y = -10.0;
220  q2.x = 0.0;
221  q2.y = 0.5;
222  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_COLINEAR );
223 
224  /* Q: Co-linear contained */
225  q1.x = 0.0;
226  q1.y = 0.4;
227  q2.x = 0.0;
228  q2.y = 0.5;
229  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_COLINEAR );
230 
231  /* Q: Horizontal touching at end point from left */
232  q1.x = -0.5;
233  q1.y = 1.0;
234  q2.x = 0.0;
235  q2.y = 1.0;
236  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
237 
238  /* Q: Horizontal touching at end point from right */
239  q1.x = 0.0;
240  q1.y = 1.0;
241  q2.x = 0.0;
242  q2.y = 0.5;
243  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_COLINEAR );
244 
245  /* Q: Horizontal touching at start point from left */
246  q1.x = 0.0;
247  q1.y = 0.0;
248  q2.x = -0.5;
249  q2.y = 0.0;
250  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_LEFT );
251 
252  /* Q: Horizontal touching at start point from right */
253  q1.x = 0.0;
254  q1.y = 0.0;
255  q2.x = 0.5;
256  q2.y = 0.0;
257  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_RIGHT );
258 
259 }
260 
262 {
263 
264  POINT4D p;
265 
266  /*
267  ** Simple test, two two-point lines
268  */
269 
270  /* Vertical line from 0,0 to 1,1 */
271  p.x = 0.0;
272  p.y = 0.0;
273  ptarray_set_point4d(pa21, 0, &p);
274  p.y = 1.0;
275  ptarray_set_point4d(pa21, 1, &p);
276 
277  /* Horizontal, crossing mid-segment */
278  p.x = -0.5;
279  p.y = 0.5;
280  ptarray_set_point4d(pa22, 0, &p);
281  p.x = 0.5;
282  ptarray_set_point4d(pa22, 1, &p);
283 
284  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_CROSS_RIGHT );
285 
286  /* Horizontal, crossing at top end vertex (end crossings don't count) */
287  p.x = -0.5;
288  p.y = 1.0;
289  ptarray_set_point4d(pa22, 0, &p);
290  p.x = 0.5;
291  ptarray_set_point4d(pa22, 1, &p);
292 
293  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_NO_CROSS );
294 
295  /* Horizontal, crossing at bottom end vertex */
296  p.x = -0.5;
297  p.y = 0.0;
298  ptarray_set_point4d(pa22, 0, &p);
299  p.x = 0.5;
300  ptarray_set_point4d(pa22, 1, &p);
301 
302  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_CROSS_RIGHT );
303 
304  /* Horizontal, no crossing */
305  p.x = -0.5;
306  p.y = 2.0;
307  ptarray_set_point4d(pa22, 0, &p);
308  p.x = 0.5;
309  ptarray_set_point4d(pa22, 1, &p);
310 
311  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_NO_CROSS );
312 
313  /* Vertical, no crossing */
314  p.x = -0.5;
315  p.y = 0.0;
316  ptarray_set_point4d(pa22, 0, &p);
317  p.y = 1.0;
318  ptarray_set_point4d(pa22, 1, &p);
319 
320  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_NO_CROSS );
321 
322 }
323 
325 {
326  LWLINE *l51;
327  LWLINE *l52;
328  /*
329  ** More complex test, longer lines and multiple crossings
330  */
331  /* Vertical line with vertices at y integers */
332  l51 = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
333 
334  /* Two crossings at segment midpoints */
335  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, -1 1.5, 1 3, 1 4, 1 5)", LW_PARSER_CHECK_NONE);
337  lwline_free(l52);
338 
339  /* One crossing at interior vertex */
340  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 0 1, -1 1, -1 2, -1 3)", LW_PARSER_CHECK_NONE);
341  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_CROSS_LEFT );
342  lwline_free(l52);
343 
344  /* Two crossings at interior vertices */
345  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 0 1, -1 1, 0 3, 1 3)", LW_PARSER_CHECK_NONE);
347  lwline_free(l52);
348 
349  /* Two crossings, one at the first vertex on at interior vertex */
350  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 0, 0 0, -1 1, 0 3, 1 3)", LW_PARSER_CHECK_NONE);
352  lwline_free(l52);
353 
354  /* Two crossings, one at the first vertex on the next interior vertex */
355  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 0, 0 0, -1 1, 0 1, 1 2)", LW_PARSER_CHECK_NONE);
357  lwline_free(l52);
358 
359  /* Three crossings, two at midpoints, one at vertex */
360  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(0.5 1, -1 0.5, 1 2, -1 2, -1 3)", LW_PARSER_CHECK_NONE);
361  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_MULTICROSS_END_LEFT );
362  lwline_free(l52);
363 
364  /* One mid-point co-linear crossing */
365  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 0 1.5, 0 2.5, -1 3, -1 4)", LW_PARSER_CHECK_NONE);
366  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_CROSS_LEFT );
367  lwline_free(l52);
368 
369  /* One on-vertices co-linear crossing */
370  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 0 1, 0 2, -1 4, -1 4)", LW_PARSER_CHECK_NONE);
371  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_CROSS_LEFT );
372  lwline_free(l52);
373 
374  /* No crossing, but end on a co-linearity. */
375  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 1 2, 1 3, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
376  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_NO_CROSS );
377  lwline_free(l52);
378 
379  lwline_free(l51);
380 
381 }
382 
383 
384 static void test_lwline_crossing_bugs(void)
385 {
386  LWLINE *l1;
387  LWLINE *l2;
388 
389  l1 = (LWLINE*)lwgeom_from_wkt("LINESTRING(2.99 90.16,71 74,20 140,171 154)", LW_PARSER_CHECK_NONE);
390  l2 = (LWLINE*)lwgeom_from_wkt("LINESTRING(25 169,89 114,40 70,86 43)", LW_PARSER_CHECK_NONE);
391 
393  lwline_free(l1);
394  lwline_free(l2);
395 
396 }
397 
398 static void test_lwpoint_set_ordinate(void)
399 {
400  POINT4D p;
401 
402  p.x = 0.0;
403  p.y = 0.0;
404  p.z = 0.0;
405  p.m = 0.0;
406 
407  lwpoint_set_ordinate(&p, 'X', 1.5);
408  CU_ASSERT_EQUAL( p.x, 1.5 );
409 
410  lwpoint_set_ordinate(&p, 'M', 2.5);
411  CU_ASSERT_EQUAL( p.m, 2.5 );
412 
413  lwpoint_set_ordinate(&p, 'Z', 3.5);
414  CU_ASSERT_EQUAL( p.z, 3.5 );
415 
416 }
417 
418 static void test_lwpoint_get_ordinate(void)
419 {
420  POINT4D p;
421 
422  p.x = 10.0;
423  p.y = 20.0;
424  p.z = 30.0;
425  p.m = 40.0;
426 
427  CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 'X'), 10.0 );
428  CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 'Y'), 20.0 );
429  CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 'Z'), 30.0 );
430  CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 'M'), 40.0 );
431 
432 }
433 
434 static void test_point_interpolate(void)
435 {
436  POINT4D p, q, r;
437  int rv = 0;
438 
439  p.x = 10.0;
440  p.y = 20.0;
441  p.z = 30.0;
442  p.m = 40.0;
443 
444  q.x = 20.0;
445  q.y = 30.0;
446  q.z = 40.0;
447  q.m = 50.0;
448 
449  rv = point_interpolate(&p, &q, &r, 1, 1, 'Z', 35.0);
450  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
451  CU_ASSERT_EQUAL( r.x, 15.0);
452 
453  rv = point_interpolate(&p, &q, &r, 1, 1, 'M', 41.0);
454  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
455  CU_ASSERT_EQUAL( r.y, 21.0);
456 
457  rv = point_interpolate(&p, &q, &r, 1, 1, 'M', 50.0);
458  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
459  CU_ASSERT_EQUAL( r.y, 30.0);
460 
461  rv = point_interpolate(&p, &q, &r, 1, 1, 'M', 40.0);
462  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
463  CU_ASSERT_EQUAL( r.y, 20.0);
464 
465 }
466 
467 static void test_lwline_clip(void)
468 {
469  LWCOLLECTION *c;
470  LWLINE *line = NULL;
471  LWLINE *l51 = NULL;
472  char *ewkt;
473 
474  /* Vertical line with vertices at y integers */
475  l51 = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
476 
477  /* Clip in the middle, mid-range. */
478  c = lwline_clip_to_ordinate_range(l51, 'Y', 1.5, 2.5);
479  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
480  //printf("c = %s\n", ewkt);
481  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
482  lwfree(ewkt);
484 
485  /* Clip off the top. */
486  c = lwline_clip_to_ordinate_range(l51, 'Y', 3.5, 5.5);
487  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
488  //printf("c = %s\n", ewkt);
489  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 3.5,0 4))");
490  lwfree(ewkt);
492 
493  /* Clip off the bottom. */
494  c = lwline_clip_to_ordinate_range(l51, 'Y', -1.5, 2.5);
495  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
496  //printf("c = %s\n", ewkt);
497  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 2.5))" );
498  lwfree(ewkt);
500 
501  /* Range holds entire object. */
502  c = lwline_clip_to_ordinate_range(l51, 'Y', -1.5, 5.5);
503  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
504  //printf("c = %s\n", ewkt);
505  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))" );
506  lwfree(ewkt);
508 
509  /* Clip on vertices. */
510  c = lwline_clip_to_ordinate_range(l51, 'Y', 1.0, 2.0);
511  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
512  //printf("c = %s\n", ewkt);
513  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1,0 2))" );
514  lwfree(ewkt);
516 
517  /* Clip on vertices off the bottom. */
518  c = lwline_clip_to_ordinate_range(l51, 'Y', -1.0, 2.0);
519  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
520  //printf("c = %s\n", ewkt);
521  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2))" );
522  lwfree(ewkt);
524 
525  /* Clip on top. */
526  c = lwline_clip_to_ordinate_range(l51, 'Y', -1.0, 0.0);
527  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
528  //printf("c = %s\n", ewkt);
529  CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(0 0))" );
530  lwfree(ewkt);
532 
533  /* ST_LocateBetweenElevations(ST_GeomFromEWKT('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)'), 1, 2)) */
534  line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
535  c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 2.0);
536  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
537  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))" );
538  lwfree(ewkt);
540  lwline_free(line);
541 
542  /* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 2)) */
543  line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
544  c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 2.0);
545  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
546  //printf("a = %s\n", ewkt);
547  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))" );
548  lwfree(ewkt);
550  lwline_free(line);
551 
552  /* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 1)) */
553  line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
554  c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 1.0);
555  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
556  //printf("b = %s\n", ewkt);
557  CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))" );
558  lwfree(ewkt);
560  lwline_free(line);
561 
562  /* ST_LocateBetweenElevations('LINESTRING(1 1 1, 1 2 2)', 1,1) */
563  line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1 1, 1 2 2)", LW_PARSER_CHECK_NONE);
564  c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 1.0);
565  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
566  //printf("c = %s\n", ewkt);
567  CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))" );
568  lwfree(ewkt);
570  lwline_free(line);
571 
572  lwline_free(l51);
573 
574 }
575 
576 static void test_lwmline_clip(void)
577 {
578  LWCOLLECTION *c;
579  char *ewkt;
580  LWMLINE *mline = NULL;
581  LWLINE *line = NULL;
582 
583  /*
584  ** Set up the input line. Trivial one-member case.
585  */
586  mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
587 
588  /* Clip in the middle, mid-range. */
589  c = lwmline_clip_to_ordinate_range(mline, 'Y', 1.5, 2.5);
590  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
591  //printf("c = %s\n", ewkt);
592  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
593  lwfree(ewkt);
595 
596  lwmline_free(mline);
597 
598  /*
599  ** Set up the input line. Two-member case.
600  */
601  mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((1 0,1 1,1 2,1 3,1 4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
602 
603  /* Clip off the top. */
604  c = lwmline_clip_to_ordinate_range(mline, 'Y', 3.5, 5.5);
605  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
606  //printf("c = %s\n", ewkt);
607  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((1 3.5,1 4),(0 3.5,0 4))");
608  lwfree(ewkt);
610 
611  lwmline_free(mline);
612 
613  /*
614  ** Set up staggered input line to create multi-type output.
615  */
616  mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((1 0,1 -1,1 -2,1 -3,1 -4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
617 
618  /* Clip from 0 upwards.. */
619  c = lwmline_clip_to_ordinate_range(mline, 'Y', 0.0, 2.5);
620  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
621  //printf("c = %s\n", ewkt);
622  CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 0),LINESTRING(0 0,0 1,0 2,0 2.5))");
623  lwfree(ewkt);
625 
626  lwmline_free(mline);
627 
628  /*
629  ** Set up input line from MAC
630  */
631  line = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0 0 0,1 1 1 1,2 2 2 2,3 3 3 3,4 4 4 4,3 3 3 5,2 2 2 6,1 1 1 7,0 0 0 8)", LW_PARSER_CHECK_NONE);
632 
633  /* Clip from 3 to 3.5 */
634  c = lwline_clip_to_ordinate_range(line, 'Z', 3.0, 3.5);
635  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
636  //printf("c = %s\n", ewkt);
637  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,3.5 3.5 3.5 3.5),(3.5 3.5 3.5 4.5,3 3 3 5))");
638  lwfree(ewkt);
640 
641  /* Clip from 2 to 3.5 */
642  c = lwline_clip_to_ordinate_range(line, 'Z', 2.0, 3.5);
643  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
644  //printf("c = %s\n", ewkt);
645  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3,3.5 3.5 3.5 3.5),(3.5 3.5 3.5 4.5,3 3 3 5,2 2 2 6))");
646  lwfree(ewkt);
648 
649  /* Clip from 3 to 4 */
650  c = lwline_clip_to_ordinate_range(line, 'Z', 3.0, 4.0);
651  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
652  //printf("c = %s\n", ewkt);
653  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,4 4 4 4,3 3 3 5))");
654  lwfree(ewkt);
656 
657  /* Clip from 2 to 3 */
658  c = lwline_clip_to_ordinate_range(line, 'Z', 2.0, 3.0);
659  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
660  //printf("c = %s\n", ewkt);
661  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3),(3 3 3 5,2 2 2 6))");
662  lwfree(ewkt);
664 
665 
666  lwline_free(line);
667 
668 }
669 
670 
671 
672 static void test_lwline_clip_big(void)
673 {
674  POINTARRAY *pa = ptarray_construct(1, 0, 3);
675  LWLINE *line = lwline_construct(SRID_UNKNOWN, NULL, pa);
676  LWCOLLECTION *c;
677  char *ewkt;
678  POINT4D p;
679 
680  p.x = 0.0;
681  p.y = 0.0;
682  p.z = 0.0;
683  ptarray_set_point4d(pa, 0, &p);
684 
685  p.x = 1.0;
686  p.y = 1.0;
687  p.z = 1.0;
688  ptarray_set_point4d(pa, 1, &p);
689 
690  p.x = 2.0;
691  p.y = 2.0;
692  p.z = 2.0;
693  ptarray_set_point4d(pa, 2, &p);
694 
695  c = lwline_clip_to_ordinate_range(line, 'Z', 0.5, 1.5);
696  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
697  //printf("c = %s\n", ewkt);
698  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0.5 0.5 0.5,1 1 1,1.5 1.5 1.5))" );
699 
700  lwfree(ewkt);
702  lwline_free(line);
703 }
704 
705 static void test_geohash_precision(void)
706 {
707  GBOX bbox;
708  GBOX bounds;
709  int precision = 0;
710  gbox_init(&bbox);
711  gbox_init(&bounds);
712 
713  bbox.xmin = 23.0;
714  bbox.xmax = 23.0;
715  bbox.ymin = 25.2;
716  bbox.ymax = 25.2;
717  precision = lwgeom_geohash_precision(bbox, &bounds);
718  //printf("\nprecision %d\n",precision);
719  CU_ASSERT_EQUAL(precision, 20);
720 
721  bbox.xmin = 23.0;
722  bbox.ymin = 23.0;
723  bbox.xmax = 23.1;
724  bbox.ymax = 23.1;
725  precision = lwgeom_geohash_precision(bbox, &bounds);
726  //printf("precision %d\n",precision);
727  CU_ASSERT_EQUAL(precision, 3);
728 
729  bbox.xmin = 23.0;
730  bbox.ymin = 23.0;
731  bbox.xmax = 23.0001;
732  bbox.ymax = 23.0001;
733  precision = lwgeom_geohash_precision(bbox, &bounds);
734  //printf("precision %d\n",precision);
735  CU_ASSERT_EQUAL(precision, 7);
736 
737 }
738 
739 static void test_geohash_point(void)
740 {
741  char *geohash;
742 
743  geohash = geohash_point(0, 0, 16);
744  //printf("\ngeohash %s\n",geohash);
745  CU_ASSERT_STRING_EQUAL(geohash, "s000000000000000");
746  lwfree(geohash);
747 
748  geohash = geohash_point(90, 0, 16);
749  //printf("\ngeohash %s\n",geohash);
750  CU_ASSERT_STRING_EQUAL(geohash, "w000000000000000");
751  lwfree(geohash);
752 
753  geohash = geohash_point(20.012345, -20.012345, 15);
754  //printf("\ngeohash %s\n",geohash);
755  CU_ASSERT_STRING_EQUAL(geohash, "kkqnpkue9ktbpe5");
756  lwfree(geohash);
757 
758 }
759 
760 static void test_geohash(void)
761 {
762  LWPOINT *lwpoint = NULL;
763  LWLINE *lwline = NULL;
764  LWMLINE *lwmline = NULL;
765  char *geohash = NULL;
766 
767  lwpoint = (LWPOINT*)lwgeom_from_wkt("POINT(23.0 25.2)", LW_PARSER_CHECK_NONE);
768  geohash = lwgeom_geohash((LWGEOM*)lwpoint,0);
769  //printf("\ngeohash %s\n",geohash);
770  CU_ASSERT_STRING_EQUAL(geohash, "ss2r77s0du7p2ewb8hmx");
771  lwpoint_free(lwpoint);
772  lwfree(geohash);
773 
774  lwpoint = (LWPOINT*)lwgeom_from_wkt("POINT(23.0 25.2 2.0)", LW_PARSER_CHECK_NONE);
775  geohash = lwgeom_geohash((LWGEOM*)lwpoint,0);
776  //printf("geohash %s\n",geohash);
777  CU_ASSERT_STRING_EQUAL(geohash, "ss2r77s0du7p2ewb8hmx");
778  lwpoint_free(lwpoint);
779  lwfree(geohash);
780 
781  lwline = (LWLINE*)lwgeom_from_wkt("LINESTRING(23.0 23.0,23.1 23.1)", LW_PARSER_CHECK_NONE);
782  geohash = lwgeom_geohash((LWGEOM*)lwline,0);
783  //printf("geohash %s\n",geohash);
784  CU_ASSERT_STRING_EQUAL(geohash, "ss0");
785  lwline_free(lwline);
786  lwfree(geohash);
787 
788  lwline = (LWLINE*)lwgeom_from_wkt("LINESTRING(23.0 23.0,23.001 23.001)", LW_PARSER_CHECK_NONE);
789  geohash = lwgeom_geohash((LWGEOM*)lwline,0);
790  //printf("geohash %s\n",geohash);
791  CU_ASSERT_STRING_EQUAL(geohash, "ss06g7");
792  lwline_free(lwline);
793  lwfree(geohash);
794 
795  lwmline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((23.0 23.0,23.1 23.1),(23.0 23.0,23.1 23.1))", LW_PARSER_CHECK_NONE);
796  geohash = lwgeom_geohash((LWGEOM*)lwmline,0);
797  //printf("geohash %s\n",geohash);
798  CU_ASSERT_STRING_EQUAL(geohash, "ss0");
799  lwmline_free(lwmline);
800  lwfree(geohash);
801 }
802 
803 static void test_isclosed(void)
804 {
805  LWGEOM *geom;
806 
807  /* LINESTRING */
808 
809  /* Not Closed on 2D */
810  geom = lwgeom_from_wkt("LINESTRING(1 2,3 4)", LW_PARSER_CHECK_NONE);
811  CU_ASSERT(!lwline_is_closed((LWLINE *) geom));
812  lwgeom_free(geom);
813 
814  /* Closed on 2D */
815  geom = lwgeom_from_wkt("LINESTRING(1 2,3 4,1 2)", LW_PARSER_CHECK_NONE);
816  CU_ASSERT(lwline_is_closed((LWLINE *) geom));
817  lwgeom_free(geom);
818 
819  /* Not closed on 3D */
820  geom = lwgeom_from_wkt("LINESTRING(1 2 3,4 5 6)", LW_PARSER_CHECK_NONE);
821  CU_ASSERT(!lwline_is_closed((LWLINE *) geom));
822  lwgeom_free(geom);
823 
824  /* Closed on 3D */
825  geom = lwgeom_from_wkt("LINESTRING(1 2 3,4 5 6,1 2 3)", LW_PARSER_CHECK_NONE);
826  CU_ASSERT(lwline_is_closed((LWLINE *) geom));
827  lwgeom_free(geom);
828 
829  /* Closed on 4D, even if M is not the same */
830  geom = lwgeom_from_wkt("LINESTRING(1 2 3 4,5 6 7 8,1 2 3 0)", LW_PARSER_CHECK_NONE);
831  CU_ASSERT(lwline_is_closed((LWLINE *) geom));
832  lwgeom_free(geom);
833 
834 
835  /* CIRCULARSTRING */
836 
837  /* Not Closed on 2D */
838  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2,3 4,5 6)", LW_PARSER_CHECK_NONE);
839  CU_ASSERT(!lwcircstring_is_closed((LWCIRCSTRING *) geom));
840  lwgeom_free(geom);
841 
842  /* Closed on 2D */
843  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2,3 4,1 2)", LW_PARSER_CHECK_NONE);
844  CU_ASSERT(lwcircstring_is_closed((LWCIRCSTRING *) geom));
845  lwgeom_free(geom);
846 
847  /* Not closed on 3D */
848  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2 3,4 5 6,7 8 9)", LW_PARSER_CHECK_NONE);
849  CU_ASSERT(!lwcircstring_is_closed((LWCIRCSTRING *) geom));
850  lwgeom_free(geom);
851 
852  /* Closed on 3D */
853  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2 3,4 5 6,1 2 3)", LW_PARSER_CHECK_NONE);
854  CU_ASSERT(lwcircstring_is_closed((LWCIRCSTRING *) geom));
855  lwgeom_free(geom);
856 
857  /* Closed on 4D, even if M is not the same */
858  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2 3 4,5 6 7 8,1 2 3 0)", LW_PARSER_CHECK_NONE);
859  CU_ASSERT(lwcircstring_is_closed((LWCIRCSTRING *) geom));
860  lwgeom_free(geom);
861 
862 
863  /* COMPOUNDCURVE */
864 
865  /* Not Closed on 2D */
866  geom = lwgeom_from_wkt("COMPOUNDCURVE(CIRCULARSTRING(1 2,3 4,1 2),(1 2,7 8,5 6))", LW_PARSER_CHECK_NONE);
867  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
868  lwgeom_free(geom);
869 
870  geom = lwgeom_from_wkt("COMPOUNDCURVE((1 2,3 4,1 2),CIRCULARSTRING(1 2,7 8,5 6))", LW_PARSER_CHECK_NONE);
871  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
872  lwgeom_free(geom);
873 
874  /* Closed on 2D */
875  geom = lwgeom_from_wkt("COMPOUNDCURVE(CIRCULARSTRING(1 2,3 4,5 6), (5 6,7 8,1 2))", LW_PARSER_CHECK_NONE);
876  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
877  lwgeom_free(geom);
878 
879  geom = lwgeom_from_wkt("COMPOUNDCURVE((1 2,3 4,5 6),CIRCULARSTRING(5 6,7 8,1 2))", LW_PARSER_CHECK_NONE);
880  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
881  lwgeom_free(geom);
882 
883  /* Not Closed on 3D */
884  geom = lwgeom_from_wkt("COMPOUNDCURVE(CIRCULARSTRING(1 2 3,4 5 6,1 2 3),(1 2 3,7 8 9,10 11 12))", LW_PARSER_CHECK_NONE);
885  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
886  lwgeom_free(geom);
887 
888  geom = lwgeom_from_wkt("COMPOUNDCURVE((1 2 3,4 5 6,1 2 3),CIRCULARSTRING(1 2 3,7 8 9,10 11 12))", LW_PARSER_CHECK_NONE);
889  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
890  lwgeom_free(geom);
891 
892  /* Closed on 3D */
893  geom = lwgeom_from_wkt("COMPOUNDCURVE(CIRCULARSTRING(1 2 3,4 5 6,7 8 9),(7 8 9,10 11 12,1 2 3))", LW_PARSER_CHECK_NONE);
894  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
895  lwgeom_free(geom);
896 
897  geom = lwgeom_from_wkt("COMPOUNDCURVE((1 2 3,4 5 6,7 8 9),CIRCULARSTRING(7 8 9,10 11 12,1 2 3))", LW_PARSER_CHECK_NONE);
898  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
899  lwgeom_free(geom);
900 
901  /* Closed on 4D, even if M is not the same */
902  geom = lwgeom_from_wkt("COMPOUNDCURVE((1 2 3 4,5 6 7 8,9 10 11 12),CIRCULARSTRING(9 10 11 12,13 14 15 16,1 2 3 0))", LW_PARSER_CHECK_NONE);
903  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
904  lwgeom_free(geom);
905 }
906 
907 
908 static void test_geohash_point_as_int(void)
909 {
910  unsigned int gh;
911  POINT2D p;
912  unsigned long long rs;
913 
914  p.x = 50; p.y = 35;
915  gh = geohash_point_as_int(&p);
916  rs = 3440103613;
917  CU_ASSERT_EQUAL(gh, rs);
918  p.x = 140; p.y = 45;
919  gh = geohash_point_as_int(&p);
920  rs = 3982480893;
921  CU_ASSERT_EQUAL(gh, rs);
922  p.x = 140; p.y = 55;
923  gh = geohash_point_as_int(&p);
924  rs = 4166944232;
925  CU_ASSERT_EQUAL(gh, rs);
926 }
927 
928 static void test_lwgeom_simplify(void)
929 {
930  LWGEOM *l;
931  LWGEOM *g;
932  char *ewkt;
933 
934  /* Simplify but only so far... */
935  g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
936  l = lwgeom_simplify(g, 10, LW_TRUE);
937  ewkt = lwgeom_to_ewkt(l);
938  CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)");
939  lwgeom_free(g);
940  lwgeom_free(l);
941  lwfree(ewkt);
942 
943  /* Simplify but only so far... */
944  g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
945  l = lwgeom_simplify(g, 10, LW_TRUE);
946  ewkt = lwgeom_to_ewkt(l);
947  CU_ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))");
948  lwgeom_free(g);
949  lwgeom_free(l);
950  lwfree(ewkt);
951 
952  /* Simplify and collapse */
953  g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
954  l = lwgeom_simplify(g, 10, LW_FALSE);
955  CU_ASSERT_EQUAL(l, NULL);
956  lwgeom_free(g);
957  lwgeom_free(l);
958 
959  /* Simplify and collapse */
960  g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
961  l = lwgeom_simplify(g, 10, LW_FALSE);
962  CU_ASSERT_EQUAL(l, NULL);
963  lwgeom_free(g);
964  lwgeom_free(l);
965 
966  /* Not simplifiable */
967  g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE);
968  l = lwgeom_simplify(g, 1.0, LW_FALSE);
969  ewkt = lwgeom_to_ewkt(l);
970  CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)");
971  lwgeom_free(g);
972  lwgeom_free(l);
973  lwfree(ewkt);
974 
975  /* Simplifiable */
976  g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE);
977  l = lwgeom_simplify(g, 1.0, LW_FALSE);
978  ewkt = lwgeom_to_ewkt(l);
979  CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)");
980  lwgeom_free(g);
981  lwgeom_free(l);
982  lwfree(ewkt);
983 }
984 
985 
986 static void do_median_dims_check(char* wkt, int expected_dims)
987 {
989  LWPOINT* result = lwgeom_median(g, 1e-8, 100, LW_FALSE);
990 
991  CU_ASSERT_EQUAL(expected_dims, lwgeom_ndims((LWGEOM*) result));
992 
993  lwgeom_free(g);
994  lwpoint_free(result);
995 }
996 
998 {
999  do_median_dims_check("MULTIPOINT ((1 3), (4 7), (2 9), (0 4), (2 2))", 2);
1000  do_median_dims_check("MULTIPOINT Z ((1 3 4), (4 7 8), (2 9 1), (0 4 4), (2 2 3))", 3);
1001  do_median_dims_check("MULTIPOINT M ((1 3 4), (4 7 8), (2 9 1), (0 4 4), (2 2 3))", 2);
1002  do_median_dims_check("MULTIPOINT ZM ((1 3 4 5), (4 7 8 6), (2 9 1 7), (0 4 4 8), (2 2 3 9))", 3);
1003 }
1004 
1005 static void do_median_test(char* input, char* expected)
1006 {
1008  LWPOINT* expected_result = lwgeom_as_lwpoint(lwgeom_from_wkt(expected, LW_PARSER_CHECK_NONE));
1009  POINT3DZ actual_pt;
1010  POINT3DZ expected_pt;
1011 
1012  LWPOINT* result = lwgeom_median(g, FP_TOLERANCE / 10.0, 1000, LW_TRUE);
1013  int passed = LW_TRUE;
1014 
1015  lwpoint_getPoint3dz_p(result, &actual_pt);
1016  lwpoint_getPoint3dz_p(expected_result, &expected_pt);
1017 
1018  passed = passed && lwgeom_is_empty((LWGEOM*) expected_result) == lwgeom_is_empty((LWGEOM*) result);
1019  passed = passed && (lwgeom_has_z((LWGEOM*) expected_result) == lwgeom_has_z((LWGEOM*) result));
1020 
1021  if (!lwgeom_is_empty((LWGEOM*) result))
1022  {
1023  passed = passed && FP_EQUALS(actual_pt.x, expected_pt.x);
1024  passed = passed && FP_EQUALS(actual_pt.y, expected_pt.y);
1025  passed = passed && (!lwgeom_has_z((LWGEOM*) expected_result) || FP_EQUALS(actual_pt.z, expected_pt.z));
1026  }
1027 
1028  if (!passed)
1029  printf("median_test expected %s got %s\n", lwgeom_to_ewkt((LWGEOM*) expected_result), lwgeom_to_ewkt((LWGEOM*) result));
1030 
1031  CU_ASSERT_TRUE(passed);
1032 
1033  lwgeom_free(g);
1034  lwpoint_free(expected_result);
1035  lwpoint_free(result);
1036 }
1037 
1038 static void test_median_robustness(void)
1039 {
1040  /* A simple implementation of Weiszfeld's algorithm will fail if the median is equal
1041  * to any one of the inputs, during any iteration of the algorithm.
1042  *
1043  * Because the algorithm uses the centroid as a starting point, this situation will
1044  * occur in the test case below.
1045  */
1046  do_median_test("MULTIPOINT ((0 -1), (0 0), (0 1))", "POINT (0 0)");
1047 
1048  /* Same as above but 3D, and shifter */
1049  do_median_test("MULTIPOINT ((1 -1 3), (1 0 2), (2 1 1))", "POINT (1 0 2)");
1050 
1051  /* Starting point is duplicated */
1052  do_median_test("MULTIPOINT ((0 -1), (0 0), (0 0), (0 1))", "POINT (0 0)");
1053 
1054  /* Cube */
1055  do_median_test("MULTIPOINT ((10 10 10), (10 20 10), (20 10 10), (20 20 10), (10 10 20), (10 20 20), (20 10 20), (20 20 20))",
1056  "POINT (15 15 15)");
1057 
1058  /* Some edge cases */
1059  do_median_test("MULTIPOINT EMPTY", "POINT EMPTY");
1060  do_median_test("MULTIPOINT (EMPTY)", "POINT EMPTY");
1061  do_median_test("POINT (7 6)", "POINT (7 6)");
1062  do_median_test("POINT (7 6 2)", "POINT (7 6 2)");
1063  do_median_test("MULTIPOINT ((7 6 2), EMPTY)", "POINT (7 6 2)");
1064 }
1065 
1066 static void test_point_density(void)
1067 {
1068  LWGEOM *geom;
1069  LWMPOINT *mpt;
1070  // char *ewkt;
1071 
1072  /* POLYGON */
1073  geom = lwgeom_from_wkt("POLYGON((1 0,0 1,1 2,2 1,1 0))", LW_PARSER_CHECK_NONE);
1074  mpt = lwgeom_to_points(geom, 100);
1075  CU_ASSERT_EQUAL(mpt->ngeoms,100);
1076  // ewkt = lwgeom_to_ewkt((LWGEOM*)mpt);
1077  // printf("%s\n", ewkt);
1078  // lwfree(ewkt);
1079  lwmpoint_free(mpt);
1080 
1081  mpt = lwgeom_to_points(geom, 1);
1082  CU_ASSERT_EQUAL(mpt->ngeoms,1);
1083  lwmpoint_free(mpt);
1084 
1085  mpt = lwgeom_to_points(geom, 0);
1086  CU_ASSERT_EQUAL(mpt, NULL);
1087  lwmpoint_free(mpt);
1088 
1089  lwgeom_free(geom);
1090 
1091  /* MULTIPOLYGON */
1092  geom = lwgeom_from_wkt("MULTIPOLYGON(((10 0,0 10,10 20,20 10,10 0)),((0 0,5 0,5 5,0 5,0 0)))", LW_PARSER_CHECK_NONE);
1093 
1094  mpt = lwgeom_to_points(geom, 1000);
1095  CU_ASSERT_EQUAL(mpt->ngeoms,1000);
1096  lwmpoint_free(mpt);
1097 
1098  mpt = lwgeom_to_points(geom, 1);
1099  CU_ASSERT_EQUAL(mpt->ngeoms,1);
1100  lwmpoint_free(mpt);
1101 
1102  lwgeom_free(geom);
1103 }
1104 
1106 {
1107  LWPOLY* p;
1108  GBOX* g;
1109  const int srid = 4326;
1110  const int segments_per_quad = 17;
1111  const int x = 10;
1112  const int y = 20;
1113  const int r = 5;
1114 
1115  /* With normal arguments you should get something circle-y */
1116  p = lwpoly_construct_circle(srid, x, y, r, segments_per_quad, LW_TRUE);
1117 
1118  ASSERT_INT_EQUAL(lwgeom_count_vertices(p), segments_per_quad * 4 + 1);
1120 
1122  CU_ASSERT_DOUBLE_EQUAL(g->xmin, x-r, 0.1);
1123  CU_ASSERT_DOUBLE_EQUAL(g->xmax, x+r, 0.1);
1124  CU_ASSERT_DOUBLE_EQUAL(g->ymin, y-r, 0.1);
1125  CU_ASSERT_DOUBLE_EQUAL(g->ymax, y+r, 0.1);
1126 
1127  CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(p), M_PI*5*5, 0.1);
1128 
1129  lwpoly_free(p);
1130 
1131  /* No segments? No circle. */
1132  p = lwpoly_construct_circle(srid, x, y, r, 0, LW_TRUE);
1133  CU_ASSERT_TRUE(p == NULL);
1134 
1135  /* Negative radius? No circle. */
1136  p = lwpoly_construct_circle(srid, x, y, -1e-3, segments_per_quad, LW_TRUE);
1137  CU_ASSERT_TRUE(p == NULL);
1138 
1139  /* Zero radius? Invalid circle */
1140  p = lwpoly_construct_circle(srid, x, y, 0, segments_per_quad, LW_TRUE);
1141  CU_ASSERT_TRUE(p != NULL);
1142  lwpoly_free(p);
1143 }
1144 
1145 static void test_kmeans(void)
1146 {
1147  static int cluster_size = 100;
1148  static int num_clusters = 4;
1149  static double spread = 1.5;
1150  int N = cluster_size * num_clusters;
1151  LWGEOM **geoms;
1152  int i, j, k=0;
1153  int *r;
1154 
1155  geoms = lwalloc(sizeof(LWGEOM*) * N);
1156 
1157  for (j = 0; j < num_clusters; j++) {
1158  for (i = 0; i < cluster_size; i++)
1159  {
1160  double u1 = 1.0 * rand() / RAND_MAX;
1161  double u2 = 1.0 * rand() / RAND_MAX;
1162  double z1 = spread * j + sqrt(-2*log2(u1))*cos(2*M_PI*u2);
1163  double z2 = spread * j + sqrt(-2*log2(u1))*sin(2*M_PI*u2);
1164 
1165  LWPOINT *lwp = lwpoint_make2d(SRID_UNKNOWN, z1, z2);
1166  geoms[k++] = lwpoint_as_lwgeom(lwp);
1167  }
1168  }
1169 
1170  r = lwgeom_cluster_2d_kmeans((const LWGEOM **)geoms, N, num_clusters);
1171 
1172  // for (i = 0; i < k; i++)
1173  // {
1174  // printf("[%d] %s\n", r[i], lwgeom_to_ewkt(geoms[i]));
1175  // }
1176 
1177  /* Clean up */
1178  lwfree(r);
1179  for (i = 0; i < k; i++)
1180  lwgeom_free(geoms[i]);
1181  lwfree(geoms);
1182 
1183  return;
1184 }
1185 
1186 /*
1187 ** Used by test harness to register the tests in this file.
1188 */
1189 void algorithms_suite_setup(void);
1191 {
1192  CU_pSuite suite = CU_add_suite("computational_geometry", init_cg_suite, clean_cg_suite);
1206  PG_ADD_TEST(suite,test_geohash);
1208  PG_ADD_TEST(suite,test_isclosed);
1212  PG_ADD_TEST(suite,test_kmeans);
1216 }
void ptarray_set_point4d(POINTARRAY *pa, int n, const POINT4D *p4d)
Definition: lwgeom_api.c:549
double x
Definition: liblwgeom.h:351
double lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)
Determines the center of the circle defined by the three given points.
Definition: lwalgorithm.c:227
double z
Definition: liblwgeom.h:333
static void test_lwline_crossing_short_lines(void)
Definition: cu_algorithm.c:261
double y
Definition: liblwgeom.h:333
POINTARRAY * ptarray_construct(char hasz, char hasm, uint32_t npoints)
Construct an empty pointarray, allocating storage and setting the npoints, but not filling in any inf...
Definition: ptarray.c:62
static void test_lwline_clip_big(void)
Definition: cu_algorithm.c:672
static void test_lw_segment_side(void)
Definition: cu_algorithm.c:61
double m
Definition: liblwgeom.h:351
static void test_isclosed(void)
Definition: cu_algorithm.c:803
static void test_geohash_point_as_int(void)
Definition: cu_algorithm.c:908
double x
Definition: liblwgeom.h:333
char * r
Definition: cu_in_wkt.c:24
void lwmline_free(LWMLINE *mline)
Definition: lwmline.c:112
void lwfree(void *mem)
Definition: lwutil.c:242
int lwline_crossing_direction(const LWLINE *l1, const LWLINE *l2)
Given two lines, characterize how (and if) they cross each other.
Definition: lwalgorithm.c:460
POINTARRAY * pa21
Definition: cu_algorithm.c:25
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:145
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:482
double xmax
Definition: liblwgeom.h:292
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:195
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1063
static void test_lwmline_clip(void)
Definition: cu_algorithm.c:576
static void do_median_dims_check(char *wkt, int expected_dims)
Definition: cu_algorithm.c:986
#define LW_SUCCESS
Definition: liblwgeom.h:79
void lwline_free(LWLINE *line)
Definition: lwline.c:76
static void test_lw_segment_intersects(void)
Definition: cu_algorithm.c:119
static void test_lwpoint_get_ordinate(void)
Definition: cu_algorithm.c:418
double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
Given a POINT4D and an ordinate number, return the value of the ordinate.
static void test_median_robustness(void)
static void test_lwpoly_construct_circle(void)
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:835
static void test_geohash(void)
Definition: cu_algorithm.c:760
LWLINE * l22
Definition: cu_algorithm.c:28
static void test_lwline_clip(void)
Definition: cu_algorithm.c:467
LWLINE * l21
Definition: cu_algorithm.c:27
int lwcompound_is_closed(const LWCOMPOUND *curve)
Definition: lwcompound.c:35
static int clean_cg_suite(void)
Definition: cu_algorithm.c:51
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition: lwin_wkt.c:904
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:93
LWCOLLECTION * lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to)
Clip a line based on the from/to range of one of its ordinates.
static void test_geohash_point(void)
Definition: cu_algorithm.c:739
void algorithms_suite_setup(void)
static void test_lw_arc_center(void)
Definition: cu_algorithm.c:93
static void test_lwline_crossing_bugs(void)
Definition: cu_algorithm.c:384
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:849
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:252
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:863
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:1985
double x
Definition: liblwgeom.h:327
static void test_lwline_crossing_long_lines(void)
Definition: cu_algorithm.c:324
void lwmpoint_free(LWMPOINT *mpt)
Definition: lwmpoint.c:72
double ymin
Definition: liblwgeom.h:293
void gbox_init(GBOX *gbox)
Zero out all the entries in the GBOX.
Definition: g_box.c:51
double xmin
Definition: liblwgeom.h:291
static int init_cg_suite(void)
Definition: cu_algorithm.c:37
#define LW_FALSE
Definition: liblwgeom.h:76
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:150
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM...
Definition: liblwgeom.h:1992
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
int lwpoint_getPoint3dz_p(const LWPOINT *point, POINT3DZ *out)
Definition: lwpoint.c:47
int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value)
Given two points, a dimensionality, an ordinate, and an interpolation value generate a new point that...
static void test_kmeans(void)
LWGEOM * lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
Definition: lwgeom.c:1571
int lwcircstring_is_closed(const LWCIRCSTRING *curve)
Definition: lwcircstring.c:268
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:187
#define PG_ADD_TEST(suite, testfunc)
LWCOLLECTION * lwmline_clip_to_ordinate_range(const LWMLINE *mline, char ordinate, double from, double to)
Clip a multi-line based on the from/to range of one of its ordinates.
#define ASSERT_INT_EQUAL(o, e)
LWGEOM_PARSER_RESULT parse_result
Definition: cu_algorithm.c:30
static void do_median_test(char *input, char *expected)
int lwgeom_geohash_precision(GBOX bbox, GBOX *bounds)
Definition: lwalgorithm.c:742
const GBOX * lwgeom_get_bbox(const LWGEOM *lwgeom)
Get a non-empty geometry bounding box, computing and caching it if not already there.
Definition: lwgeom.c:653
uint8_t precision
Definition: cu_in_twkb.c:25
int * lwgeom_cluster_2d_kmeans(const LWGEOM **geoms, int ngeoms, int k)
Take a list of LWGEOMs and a number of clusters and return an integer array indicating which cluster ...
Definition: lwkmeans.c:77
#define FP_TOLERANCE
Floating point comparators.
void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
Given a point, ordinate number and value, set that ordinate on the point.
double ymax
Definition: liblwgeom.h:294
double y
Definition: liblwgeom.h:327
char * geohash_point(double longitude, double latitude, int precision)
Definition: lwalgorithm.c:581
double z
Definition: liblwgeom.h:351
static void test_median_handles_3d_correctly(void)
Definition: cu_algorithm.c:997
tuple x
Definition: pixval.py:53
static void test_geohash_precision(void)
Definition: cu_algorithm.c:705
unsigned int geohash_point_as_int(POINT2D *pt)
Definition: lwalgorithm.c:647
#define setpoint(p, x1, y1)
static void test_lwpoint_set_ordinate(void)
Definition: cu_algorithm.c:398
char * lwgeom_geohash(const LWGEOM *lwgeom, int precision)
Calculate the GeoHash (http://geohash.org) string for a geometry.
Definition: lwalgorithm.c:827
double lwgeom_area(const LWGEOM *geom)
Definition: lwgeom.c:1592
int lwline_is_closed(const LWLINE *line)
Definition: lwline.c:468
static void test_point_interpolate(void)
Definition: cu_algorithm.c:434
static void test_point_density(void)
LWMPOINT * lwgeom_to_points(const LWGEOM *lwgeom, int npoints)
#define FP_EQUALS(A, B)
void lwcollection_free(LWCOLLECTION *col)
Definition: lwcollection.c:339
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:267
LWPOLY * lwpoly_construct_circle(int srid, double x, double y, double radius, uint32_t segments_per_quarter, char exterior)
Definition: lwpoly.c:98
static void test_lwgeom_simplify(void)
Definition: cu_algorithm.c:928
LWPOINT * lwgeom_median(const LWGEOM *g, double tol, uint32_t maxiter, char fail_if_not_converged)
void * lwalloc(size_t size)
Definition: lwutil.c:227
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members) ...
Definition: lwgeom.c:1310
int lwgeom_count_vertices(const LWGEOM *geom)
Count the total number of vertices in any LWGEOM.
Definition: lwgeom.c:1153
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:64
double y
Definition: liblwgeom.h:351
POINTARRAY * pa22
Definition: cu_algorithm.c:26
int ngeoms
Definition: liblwgeom.h:467
tuple y
Definition: pixval.py:54
int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
returns the kind of CG_SEGMENT_INTERSECTION_TYPE behavior of lineseg 1 (constructed from p1 and p2) a...
Definition: lwalgorithm.c:371