PostGIS  3.1.6dev-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  * Copyright 2018 Darafei Praliaskouski, me@komzpa.net
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU General Public Licence. See the COPYING file.
10  *
11  **********************************************************************/
12 
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <string.h>
16 #include "CUnit/Basic.h"
17 
18 #include "liblwgeom_internal.h"
19 #include "cu_tester.h"
20 
21 /*
22 ** Global variables used by tests below
23 */
24 
25 /* Two-point objects */
26 POINTARRAY *pa21 = NULL;
27 POINTARRAY *pa22 = NULL;
28 LWLINE *l21 = NULL;
29 LWLINE *l22 = NULL;
30 /* Parsing support */
31 LWGEOM_PARSER_RESULT parse_result;
32 
33 
34 /*
35 ** The suite initialization function.
36 ** Create any re-used objects.
37 */
38 static int init_cg_suite(void)
39 {
40  pa21 = ptarray_construct(0, 0, 2);
41  pa22 = ptarray_construct(0, 0, 2);
42  l21 = lwline_construct(SRID_UNKNOWN, NULL, pa21);
43  l22 = lwline_construct(SRID_UNKNOWN, NULL, pa22);
44  return 0;
45 
46 }
47 
48 /*
49 ** The suite cleanup function.
50 ** Frees any global objects.
51 */
52 static int clean_cg_suite(void)
53 {
54  if ( l21 ) lwline_free(l21);
55  if ( l22 ) lwline_free(l22);
56  return 0;
57 }
58 
59 /*
60 ** Test left/right side.
61 */
62 static void test_lw_segment_side(void)
63 {
64  int rv = 0;
65  POINT2D p1, p2, q;
66 
67  /* Vertical line at x=0 */
68  p1.x = 0.0;
69  p1.y = 0.0;
70  p2.x = 0.0;
71  p2.y = 1.0;
72 
73  /* On the left */
74  q.x = -2.0;
75  q.y = 1.5;
76  rv = lw_segment_side(&p1, &p2, &q);
77  //printf("left %g\n",rv);
78  CU_ASSERT(rv < 0);
79 
80  /* On the right */
81  q.x = 2.0;
82  rv = lw_segment_side(&p1, &p2, &q);
83  //printf("right %g\n",rv);
84  CU_ASSERT(rv > 0);
85 
86  /* On the line */
87  q.x = 0.0;
88  rv = lw_segment_side(&p1, &p2, &q);
89  //printf("on line %g\n",rv);
90  CU_ASSERT_EQUAL(rv, 0);
91 
92 }
93 
94 static void test_lw_arc_center(void)
95 {
96 /* double lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result); */
97  POINT2D c1;
98  double d1;
99  POINT2D p1, p2, p3;
100 
101  p1.x = 2047538.600;
102  p1.y = 7268770.435;
103  p2.x = 2047538.598;
104  p2.y = 7268770.435;
105  p3.x = 2047538.596;
106  p3.y = 7268770.436;
107 
108  d1 = lw_arc_center(&p1, &p2, &p3, &c1);
109 
110  CU_ASSERT_DOUBLE_EQUAL(d1, 0.0046097720751, 0.0001);
111  CU_ASSERT_DOUBLE_EQUAL(c1.x, 2047538.599, 0.001);
112  CU_ASSERT_DOUBLE_EQUAL(c1.y, 7268770.4395, 0.001);
113 
114  // printf("lw_arc_center: (%12.12g, %12.12g) %12.12g\n", c1.x, c1.y, d1);
115 }
116 
117 /*
118 ** Test crossings side.
119 */
120 static void test_lw_segment_intersects(void)
121 {
122 
123 #define setpoint(p, x1, y1) {(p).x = (x1); (p).y = (y1);}
124 
125  POINT2D p1, p2, q1, q2;
126 
127  /* P: Vertical line at x=0 */
128  setpoint(p1, 0.0, 0.0);
129  p1.x = 0.0;
130  p1.y = 0.0;
131  p2.x = 0.0;
132  p2.y = 1.0;
133 
134  /* Q: Horizontal line crossing left to right */
135  q1.x = -0.5;
136  q1.y = 0.5;
137  q2.x = 0.5;
138  q2.y = 0.5;
139  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_RIGHT );
140 
141  /* Q: Horizontal line crossing right to left */
142  q1.x = 0.5;
143  q1.y = 0.5;
144  q2.x = -0.5;
145  q2.y = 0.5;
146  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_LEFT );
147 
148  /* Q: Horizontal line not crossing right to left */
149  q1.x = 0.5;
150  q1.y = 1.5;
151  q2.x = -0.5;
152  q2.y = 1.5;
153  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
154 
155  /* Q: Horizontal line crossing at second vertex right to left */
156  q1.x = 0.5;
157  q1.y = 1.0;
158  q2.x = -0.5;
159  q2.y = 1.0;
160  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
161 
162  /* Q: Horizontal line crossing at first vertex right to left */
163  q1.x = 0.5;
164  q1.y = 0.0;
165  q2.x = -0.5;
166  q2.y = 0.0;
167  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_LEFT );
168 
169  /* Q: Diagonal line with large range crossing at first vertex right to left */
170  q1.x = 0.5;
171  q1.y = 10.0;
172  q2.x = -0.5;
173  q2.y = -10.0;
174  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_LEFT );
175 
176  /* Q: Diagonal line with large range crossing at second vertex right to left */
177  q1.x = 0.5;
178  q1.y = 11.0;
179  q2.x = -0.5;
180  q2.y = -9.0;
181  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
182 
183  /* Q: Horizontal touching from left at second vertex*/
184  q1.x = -0.5;
185  q1.y = 0.5;
186  q2.x = 0.0;
187  q2.y = 0.5;
188  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
189 
190  /* Q: Horizontal touching from right at first vertex */
191  q1.x = 0.0;
192  q1.y = 0.5;
193  q2.x = 0.5;
194  q2.y = 0.5;
195  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_RIGHT );
196 
197  /* Q: Horizontal touching from left and far below on second vertex */
198  q1.x = -0.5;
199  q1.y = -10.5;
200  q2.x = 0.0;
201  q2.y = 0.5;
202  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
203 
204  /* Q: Horizontal touching from right and far above on second vertex */
205  q1.x = 0.5;
206  q1.y = 10.5;
207  q2.x = 0.0;
208  q2.y = 0.5;
209  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
210 
211  /* Q: Co-linear from top */
212  q1.x = 0.0;
213  q1.y = 10.0;
214  q2.x = 0.0;
215  q2.y = 0.5;
216  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_COLINEAR );
217 
218  /* Q: Co-linear from bottom */
219  q1.x = 0.0;
220  q1.y = -10.0;
221  q2.x = 0.0;
222  q2.y = 0.5;
223  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_COLINEAR );
224 
225  /* Q: Co-linear contained */
226  q1.x = 0.0;
227  q1.y = 0.4;
228  q2.x = 0.0;
229  q2.y = 0.5;
230  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_COLINEAR );
231 
232  /* Q: Horizontal touching at end point from left */
233  q1.x = -0.5;
234  q1.y = 1.0;
235  q2.x = 0.0;
236  q2.y = 1.0;
237  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_NO_INTERSECTION );
238 
239  /* Q: Horizontal touching at end point from right */
240  q1.x = 0.0;
241  q1.y = 1.0;
242  q2.x = 0.0;
243  q2.y = 0.5;
244  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_COLINEAR );
245 
246  /* Q: Horizontal touching at start point from left */
247  q1.x = 0.0;
248  q1.y = 0.0;
249  q2.x = -0.5;
250  q2.y = 0.0;
251  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_LEFT );
252 
253  /* Q: Horizontal touching at start point from right */
254  q1.x = 0.0;
255  q1.y = 0.0;
256  q2.x = 0.5;
257  q2.y = 0.0;
258  CU_ASSERT( lw_segment_intersects(&p1, &p2, &q1, &q2) == SEG_CROSS_RIGHT );
259 
260 }
261 
262 static void test_lwline_crossing_short_lines(void)
263 {
264 
265  POINT4D p;
266 
267  /*
268  ** Simple test, two two-point lines
269  */
270 
271  /* Vertical line from 0,0 to 1,1 */
272  p.x = 0.0;
273  p.y = 0.0;
274  ptarray_set_point4d(pa21, 0, &p);
275  p.y = 1.0;
276  ptarray_set_point4d(pa21, 1, &p);
277 
278  /* Horizontal, crossing mid-segment */
279  p.x = -0.5;
280  p.y = 0.5;
281  ptarray_set_point4d(pa22, 0, &p);
282  p.x = 0.5;
283  ptarray_set_point4d(pa22, 1, &p);
284 
285  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_CROSS_RIGHT );
286 
287  /* Horizontal, crossing at top end vertex (end crossings don't count) */
288  p.x = -0.5;
289  p.y = 1.0;
290  ptarray_set_point4d(pa22, 0, &p);
291  p.x = 0.5;
292  ptarray_set_point4d(pa22, 1, &p);
293 
294  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_NO_CROSS );
295 
296  /* Horizontal, crossing at bottom end vertex */
297  p.x = -0.5;
298  p.y = 0.0;
299  ptarray_set_point4d(pa22, 0, &p);
300  p.x = 0.5;
301  ptarray_set_point4d(pa22, 1, &p);
302 
303  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_CROSS_RIGHT );
304 
305  /* Horizontal, no crossing */
306  p.x = -0.5;
307  p.y = 2.0;
308  ptarray_set_point4d(pa22, 0, &p);
309  p.x = 0.5;
310  ptarray_set_point4d(pa22, 1, &p);
311 
312  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_NO_CROSS );
313 
314  /* Vertical, no crossing */
315  p.x = -0.5;
316  p.y = 0.0;
317  ptarray_set_point4d(pa22, 0, &p);
318  p.y = 1.0;
319  ptarray_set_point4d(pa22, 1, &p);
320 
321  CU_ASSERT( lwline_crossing_direction(l21, l22) == LINE_NO_CROSS );
322 
323 }
324 
325 static void test_lwline_crossing_long_lines(void)
326 {
327  LWLINE *l51;
328  LWLINE *l52;
329  /*
330  ** More complex test, longer lines and multiple crossings
331  */
332  /* Vertical line with vertices at y integers */
333  l51 = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
334 
335  /* Two crossings at segment midpoints */
336  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, -1 1.5, 1 3, 1 4, 1 5)", LW_PARSER_CHECK_NONE);
337  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
338  lwline_free(l52);
339 
340  /* One crossing at interior vertex */
341  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 0 1, -1 1, -1 2, -1 3)", LW_PARSER_CHECK_NONE);
342  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_CROSS_LEFT );
343  lwline_free(l52);
344 
345  /* Two crossings at interior vertices */
346  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 0 1, -1 1, 0 3, 1 3)", LW_PARSER_CHECK_NONE);
347  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
348  lwline_free(l52);
349 
350  /* Two crossings, one at the first vertex on at interior vertex */
351  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 0, 0 0, -1 1, 0 3, 1 3)", LW_PARSER_CHECK_NONE);
352  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
353  lwline_free(l52);
354 
355  /* Two crossings, one at the first vertex on the next interior vertex */
356  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 0, 0 0, -1 1, 0 1, 1 2)", LW_PARSER_CHECK_NONE);
357  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_MULTICROSS_END_SAME_FIRST_LEFT );
358  lwline_free(l52);
359 
360  /* Three crossings, two at midpoints, one at vertex */
361  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(0.5 1, -1 0.5, 1 2, -1 2, -1 3)", LW_PARSER_CHECK_NONE);
362  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_MULTICROSS_END_LEFT );
363  lwline_free(l52);
364 
365  /* One mid-point co-linear crossing */
366  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 0 1.5, 0 2.5, -1 3, -1 4)", LW_PARSER_CHECK_NONE);
367  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_CROSS_LEFT );
368  lwline_free(l52);
369 
370  /* One on-vertices co-linear crossing */
371  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 0 1, 0 2, -1 4, -1 4)", LW_PARSER_CHECK_NONE);
372  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_CROSS_LEFT );
373  lwline_free(l52);
374 
375  /* No crossing, but end on a co-linearity. */
376  l52 = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1, 1 2, 1 3, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
377  CU_ASSERT( lwline_crossing_direction(l51, l52) == LINE_NO_CROSS );
378  lwline_free(l52);
379 
380  lwline_free(l51);
381 
382 }
383 
384 
385 static void test_lwline_crossing_bugs(void)
386 {
387  LWLINE *l1;
388  LWLINE *l2;
389 
390  l1 = (LWLINE*)lwgeom_from_wkt("LINESTRING(2.99 90.16,71 74,20 140,171 154)", LW_PARSER_CHECK_NONE);
391  l2 = (LWLINE*)lwgeom_from_wkt("LINESTRING(25 169,89 114,40 70,86 43)", LW_PARSER_CHECK_NONE);
392 
393  CU_ASSERT( lwline_crossing_direction(l1, l2) == LINE_MULTICROSS_END_RIGHT );
394  lwline_free(l1);
395  lwline_free(l2);
396 
397 }
398 
399 static void test_lwpoint_set_ordinate(void)
400 {
401  POINT4D p;
402 
403  p.x = 0.0;
404  p.y = 0.0;
405  p.z = 0.0;
406  p.m = 0.0;
407 
408  lwpoint_set_ordinate(&p, 'X', 1.5);
409  CU_ASSERT_EQUAL( p.x, 1.5 );
410 
411  lwpoint_set_ordinate(&p, 'M', 2.5);
412  CU_ASSERT_EQUAL( p.m, 2.5 );
413 
414  lwpoint_set_ordinate(&p, 'Z', 3.5);
415  CU_ASSERT_EQUAL( p.z, 3.5 );
416 
417 }
418 
419 static void test_lwpoint_get_ordinate(void)
420 {
421  POINT4D p;
422 
423  p.x = 10.0;
424  p.y = 20.0;
425  p.z = 30.0;
426  p.m = 40.0;
427 
428  CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 'X'), 10.0 );
429  CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 'Y'), 20.0 );
430  CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 'Z'), 30.0 );
431  CU_ASSERT_EQUAL( lwpoint_get_ordinate(&p, 'M'), 40.0 );
432 
433 }
434 
435 static void
436 test_point_interpolate(void)
437 {
438  POINT4D p, q, r = {0, 0, 0, 0};
439  int rv = 0;
440 
441  p.x = 10.0;
442  p.y = 20.0;
443  p.z = 30.0;
444  p.m = 40.0;
445 
446  q.x = 20.0;
447  q.y = 30.0;
448  q.z = 40.0;
449  q.m = 50.0;
450 
451  rv = point_interpolate(&p, &q, &r, 1, 1, 'Z', 35.0);
452  CU_ASSERT_EQUAL(rv, LW_SUCCESS);
453  CU_ASSERT_EQUAL(r.x, 15.0);
454 
455  rv = point_interpolate(&p, &q, &r, 1, 1, 'M', 41.0);
456  CU_ASSERT_EQUAL(rv, LW_SUCCESS);
457  CU_ASSERT_EQUAL(r.y, 21.0);
458 
459  rv = point_interpolate(&p, &q, &r, 1, 1, 'M', 50.0);
460  CU_ASSERT_EQUAL(rv, LW_SUCCESS);
461  CU_ASSERT_EQUAL(r.y, 30.0);
462 
463  rv = point_interpolate(&p, &q, &r, 1, 1, 'M', 40.0);
464  CU_ASSERT_EQUAL(rv, LW_SUCCESS);
465  CU_ASSERT_EQUAL(r.y, 20.0);
466 }
467 
468 static void test_lwline_interpolate_points(void)
469 {
470  LWLINE* line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING ZM (0 0 3 1, 1 1 2 2, 10 10 4 3)", LW_PARSER_CHECK_NONE));
471  LWLINE* line2d = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING (1 1, 3 7, 9 12)", LW_PARSER_CHECK_NONE));
472  LWLINE* empty_line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING EMPTY", LW_PARSER_CHECK_NONE));
473 
474  POINTARRAY* rpa;
475  POINT4D pta;
476  POINT4D ptb;
477  double eps = 1e-10;
478 
479  /* Empty line gives empty point */
480  rpa = lwline_interpolate_points(empty_line, 0.5, LW_TRUE);
481  ASSERT_INT_EQUAL(rpa->npoints, 0);
482  ptarray_free(rpa);
483 
484  /* Get first endpoint when fraction = 0 */
485  rpa = lwline_interpolate_points(line, 0, LW_TRUE);
486  ASSERT_INT_EQUAL(rpa->npoints, 1);
487  pta = getPoint4d(line->points, 0);
488  ptb = getPoint4d(rpa, 0);
489  ASSERT_POINT4D_EQUAL(pta, ptb, eps);
490  ptarray_free(rpa);
491 
492  /* Get last endpoint when fraction = 0 */
493  rpa = lwline_interpolate_points(line, 1, LW_TRUE);
494  ASSERT_INT_EQUAL(rpa->npoints, 1);
495  pta = getPoint4d(line->points, line->points->npoints - 1);
496  ptb = getPoint4d(rpa, 0);
497  ASSERT_POINT4D_EQUAL(pta, ptb, eps);
498  ptarray_free(rpa);
499 
500  /* Interpolate a single point */
501  /* First vertex is at 10% */
502  rpa = lwline_interpolate_points(line, 0.1, LW_FALSE);
503  ASSERT_INT_EQUAL(rpa->npoints, 1);
504  pta = getPoint4d(line->points, 1);
505  ptb = getPoint4d(rpa, 0);
506  ASSERT_POINT4D_EQUAL(pta, ptb, eps);
507  ptarray_free(rpa);
508 
509  /* 5% is halfway to first vertex */
510  rpa = lwline_interpolate_points(line, 0.05, LW_FALSE);
511  ASSERT_INT_EQUAL(rpa->npoints, 1);
512  pta.x = 0.5;
513  pta.y = 0.5;
514  pta.m = 1.5;
515  pta.z = 2.5;
516  ptb = getPoint4d(rpa, 0);
517  ASSERT_POINT4D_EQUAL(pta, ptb, eps);
518  ptarray_free(rpa);
519 
520  /* Now repeat points */
521  rpa = lwline_interpolate_points(line, 0.4, LW_TRUE);
522  ASSERT_INT_EQUAL(rpa->npoints, 2);
523  pta.x = 4;
524  pta.y = 4;
525  ptb = getPoint4d(rpa, 0);
526  ASSERT_POINT2D_EQUAL(pta, ptb, eps);
527 
528  pta.x = 8;
529  pta.y = 8;
530  ptb = getPoint4d(rpa, 1);
531  ASSERT_POINT2D_EQUAL(pta, ptb, eps);
532  ptarray_free(rpa);
533 
534  /* Make sure it works on 2D lines */
535  rpa = lwline_interpolate_points(line2d, 0.4, LW_TRUE);
536  ASSERT_INT_EQUAL(rpa->npoints, 2);
537  CU_ASSERT_FALSE(ptarray_has_z(rpa));
538  CU_ASSERT_FALSE(ptarray_has_m(rpa));
539  ptarray_free(rpa);
540 
541  lwgeom_free(lwline_as_lwgeom(line));
542  lwgeom_free(lwline_as_lwgeom(line2d));
543  lwgeom_free(lwline_as_lwgeom(empty_line));
544 }
545 
546 static void
547 test_lwline_interpolate_point_3d(void)
548 {
549  LWLINE *line;
550  POINT4D point;
551  LWPOINT *pt;
552 
553  /* Empty line -> Empty point*/
554  line = lwline_construct_empty(4326, LW_TRUE, LW_FALSE);
555  pt = lwline_interpolate_point_3d(line, 0.5);
556  CU_ASSERT(lwpoint_is_empty(pt));
557  CU_ASSERT(lwgeom_has_z(lwpoint_as_lwgeom(pt)));
558  CU_ASSERT_FALSE(lwgeom_has_m(lwpoint_as_lwgeom(pt)));
559  lwpoint_free(pt);
560  lwline_free(line);
561 
562  line = lwgeom_as_lwline(lwgeom_from_wkt("LINESTRING Z (0 0 0, 1 1 1, 2 2 2)", LW_PARSER_CHECK_NONE));
563 
564  /* distance = 0 -> first point */
565  pt = lwline_interpolate_point_3d(line, 0);
566  lwpoint_getPoint4d_p(pt, &point);
567  CU_ASSERT_DOUBLE_EQUAL(point.x, 0, 0.0001);
568  CU_ASSERT_DOUBLE_EQUAL(point.y, 0, 0.001);
569  CU_ASSERT_DOUBLE_EQUAL(point.z, 0, 0.001);
570  lwpoint_free(pt);
571 
572  /* distance = 1 -> last point */
573  pt = lwline_interpolate_point_3d(line, 1);
574  lwpoint_getPoint4d_p(pt, &point);
575  CU_ASSERT_DOUBLE_EQUAL(point.x, 2, 0.0001);
576  CU_ASSERT_DOUBLE_EQUAL(point.y, 2, 0.001);
577  CU_ASSERT_DOUBLE_EQUAL(point.z, 2, 0.001);
578  lwpoint_free(pt);
579 
580  /* simple where distance 50% -> second point */
581  pt = lwline_interpolate_point_3d(line, 0.5);
582  lwpoint_getPoint4d_p(pt, &point);
583  CU_ASSERT_DOUBLE_EQUAL(point.x, 1, 0.0001);
584  CU_ASSERT_DOUBLE_EQUAL(point.y, 1, 0.001);
585  CU_ASSERT_DOUBLE_EQUAL(point.z, 1, 0.001);
586  lwpoint_free(pt);
587 
588  /* simple where distance 80% -> between second and last point */
589  pt = lwline_interpolate_point_3d(line, 0.8);
590  lwpoint_getPoint4d_p(pt, &point);
591  CU_ASSERT_DOUBLE_EQUAL(point.x, 1.6, 0.0001);
592  CU_ASSERT_DOUBLE_EQUAL(point.y, 1.6, 0.001);
593  CU_ASSERT_DOUBLE_EQUAL(point.z, 1.6, 0.001);
594  lwpoint_free(pt);
595 
596  lwline_free(line);
597 }
598 
599 static void test_lwline_clip(void)
600 {
601  LWCOLLECTION *c;
602  LWGEOM *line = NULL;
603  LWGEOM *l51 = NULL;
604  char *ewkt;
605 
606  /* Vertical line with vertices at y integers */
607  l51 = lwgeom_from_wkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
608 
609  /* Clip in the middle, mid-range. */
610  c = lwgeom_clip_to_ordinate_range(l51, 'Y', 1.5, 2.5, 0);
611  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
612  //printf("c = %s\n", ewkt);
613  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
614  lwfree(ewkt);
615  lwcollection_free(c);
616 
617  /* Clip off the top. */
618  c = lwgeom_clip_to_ordinate_range(l51, 'Y', 3.5, 5.5, 0);
619  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
620  //printf("c = %s\n", ewkt);
621  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 3.5,0 4))");
622  lwfree(ewkt);
623  lwcollection_free(c);
624 
625  /* Clip off the bottom. */
626  c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.5, 2.5, 0);
627  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
628  //printf("c = %s\n", ewkt);
629  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 2.5))");
630  lwfree(ewkt);
631  lwcollection_free(c);
632 
633  /* Range holds entire object. */
634  c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.5, 5.5, 0);
635  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
636  //printf("c = %s\n", ewkt);
637  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))");
638  lwfree(ewkt);
639  lwcollection_free(c);
640 
641  /* Clip on vertices. */
642  c = lwgeom_clip_to_ordinate_range(l51, 'Y', 1.0, 2.0, 0);
643  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
644  //printf("c = %s\n", ewkt);
645  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1,0 2))");
646  lwfree(ewkt);
647  lwcollection_free(c);
648 
649  /* Clip on vertices off the bottom. */
650  c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.0, 2.0, 0);
651  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
652  //printf("c = %s\n", ewkt);
653  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2))");
654  lwfree(ewkt);
655  lwcollection_free(c);
656 
657  /* Clip on top. */
658  c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.0, 0.0, 0);
659  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
660  //printf("c = %s\n", ewkt);
661  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(0 0))");
662  lwfree(ewkt);
663  lwcollection_free(c);
664 
665  /* ST_LocateBetweenElevations(ST_GeomFromEWKT('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)'), 1, 2)) */
666  line = lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
667  c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 2.0, 0);
668  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
669  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))");
670  lwfree(ewkt);
671  lwcollection_free(c);
672  lwgeom_free(line);
673 
674  /* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 2)) */
675  line = lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
676  c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 2.0, 0);
677  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
678  //printf("a = %s\n", ewkt);
679  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))");
680  lwfree(ewkt);
681  lwcollection_free(c);
682  lwgeom_free(line);
683 
684  /* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 1)) */
685  line = lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
686  c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 1.0, 0);
687  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
688  //printf("b = %s\n", ewkt);
689  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))");
690  lwfree(ewkt);
691  lwcollection_free(c);
692  lwgeom_free(line);
693 
694  /* ST_LocateBetweenElevations('LINESTRING(1 1 1, 1 2 2)', 1,1) */
695  line = lwgeom_from_wkt("LINESTRING(1 1 1, 1 2 2)", LW_PARSER_CHECK_NONE);
696  c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 1.0, 0);
697  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
698  //printf("c = %s\n", ewkt);
699  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))");
700  lwfree(ewkt);
701  lwcollection_free(c);
702  lwgeom_free(line);
703 
704  lwgeom_free(l51);
705 }
706 
707 static void
708 test_lwpoly_clip(void)
709 {
710  LWCOLLECTION *c;
711  LWGEOM *g = NULL;
712  char *ewkt;
713 
714  g = lwgeom_from_wkt(
715  "POLYGON ((0.51 -0.25, 1.27 -0.14, 1.27 0.25, 0.6 0.3, 0.7 0.7, 1.2 0.7, 0.8 0.5, 1.3 0.4, 1.2 1.2, 0.5 1.2, 0.5 -0.1, 0.3 -0.1, 0.3 1.3, -0.18 1.25, -0.17 -0.25, 0.51 -0.25))",
716  LW_PARSER_CHECK_NONE);
717  c = lwgeom_clip_to_ordinate_range(g, 'X', 0.0, 1.0, 0);
718 
719  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
720  // printf("c = %s\n", ewkt);
721  ASSERT_STRING_EQUAL(
722  ewkt,
723  "MULTIPOLYGON(((0.51 -0.25,1 -0.179078947368,1 0.270149253731,0.6 0.3,0.7 0.7,1 0.7,1 0.6,0.8 0.5,1 0.46,1 1.2,0.5 1.2,0.5 -0.1,0.3 -0.1,0.3 1.3,0 1.26875,0 -0.25,0.51 -0.25)))");
724  lwfree(ewkt);
725  lwcollection_free(c);
726  lwgeom_free(g);
727 
728  g = lwgeom_from_wkt(
729  "MULTIPOLYGON(((0.51 -0.25,1 -0.179078947368,1 0.270149253731,0.6 0.3,0.7 0.7,1 0.7,1 0.6,0.8 0.5,1 0.46,1 1.2,0.5 1.2,0.5 -0.1,0.3 -0.1,0.3 1.3,0 1.26875,0 -0.25,0.51 -0.25)))",
730  LW_PARSER_CHECK_NONE);
731  c = lwgeom_clip_to_ordinate_range(g, 'Y', 0.0, 1.0, 0);
732 
733  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
734  //printf("c = %s\n", ewkt);
735  ASSERT_STRING_EQUAL(
736  ewkt,
737  "MULTIPOLYGON(((1 0,1 0.270149253731,0.6 0.3,0.7 0.7,1 0.7,1 0.6,0.8 0.5,1 0.46,1 1,0.5 1,0.5 0,0.3 0,0.3 1,0 1,0 0,1 0)))");
738  lwfree(ewkt);
739  lwcollection_free(c);
740  lwgeom_free(g);
741 }
742 
743 static void
744 test_lwtriangle_clip(void)
745 {
746  LWCOLLECTION *c;
747  LWGEOM *g = NULL;
748  char *ewkt;
749 
750  g = lwgeom_from_wkt("TRIANGLE((0 0 0, 1 1 1, 3 2 2, 0 0 0))", LW_PARSER_CHECK_NONE);
751  c = lwgeom_clip_to_ordinate_range(g, 'Z', -10.0, 4.0, 0);
752 
753  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
754  // printf("c = %s\n", ewkt);
755  ASSERT_STRING_EQUAL(ewkt, "TIN(((0 0 0,1 1 1,3 2 2,0 0 0)))");
756  lwfree(ewkt);
757  lwcollection_free(c);
758  lwgeom_free(g);
759 
760  g = lwgeom_from_wkt("TRIANGLE((0 0 0, 1 1 1, 3 2 2, 0 0 0))", LW_PARSER_CHECK_NONE);
761  c = lwgeom_clip_to_ordinate_range(g, 'Z', 0.0, 1.0, 0);
762 
763  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
764  // printf("c = %s\n", ewkt);
765  ASSERT_STRING_EQUAL(ewkt, "TIN(((0 0 0,1 1 1,1.5 1 1,0 0 0)))");
766  lwfree(ewkt);
767  lwcollection_free(c);
768  lwgeom_free(g);
769 
770  g = lwgeom_from_wkt("TRIANGLE((0 0 0, 1 1 1, 3 2 3, 0 0 0))", LW_PARSER_CHECK_NONE);
771  c = lwgeom_clip_to_ordinate_range(g, 'Z', 1.0, 2.0, 0);
772 
773  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
774  // printf("c = %s\n", ewkt);
775  ASSERT_STRING_EQUAL(
776  ewkt,
777  "TIN(((1 1 1,2 1.5 2,2 1.333333333333 2,1 1 1)),((1 1 1,2 1.333333333333 2,1 0.666666666667 1,1 1 1)))");
778  lwfree(ewkt);
779  lwcollection_free(c);
780  lwgeom_free(g);
781 
782  g = lwgeom_from_wkt(
783  "TIN Z (((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 -1 -1,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,-1 2 -1,-1 2 2)),((-1 -1 -1,-1 2 2,-1 2 -1,-1 -1 -1)),((-1 -1 2,-1 -1 -1,2 -1 -1,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 2 -1,2 2 2)),((-1 2 -1,2 2 2,2 2 -1,-1 2 -1)),((-1 -1 -1,-1 2 -1,2 2 -1,-1 -1 -1)),((2 -1 -1,-1 -1 -1,2 2 -1,2 -1 -1)),((-1 -1 2,2 -1 -1,2 -1 2,-1 -1 2)),((2 2 -1,2 -1 2,2 -1 -1,2 2 -1)))",
784  LW_PARSER_CHECK_NONE);
785  c = lwgeom_clip_to_ordinate_range(g, 'Z', 0.0, DBL_MAX, 0);
786 
787  /* Adjust for Rasberry Pi 64-bit (Debian Buster) */
788  gridspec grid = {0};
789  grid.xsize = grid.ysize = grid.zsize = grid.msize = 1;
790  lwgeom_grid_in_place((LWGEOM *)c, &grid);
791 
792  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
793  // printf("c = %s\n", ewkt);
794  ASSERT_STRING_EQUAL(
795  ewkt,
796  "TIN(((0 0 2,0 1 2,-1 2 2,0 0 2)),((0 1 2,0 0 2,0 1 4,0 1 2)),((-1 2 2,0 1 2,1 1 2,-1 2 2)),((0 0 2,-1 2 2,-1 -1 2,0 0 2)),((0 1 4,0 0 2,0 0 4,0 1 4)),((0 1 2,0 1 4,1 1 4,0 1 2)),((1 1 2,0 1 2,1 1 4,1 1 2)),((-1 2 2,1 1 2,2 2 2,-1 2 2)),((-1 -1 2,-1 2 2,-1 0 0,-1 -1 2)),((-1 -1 2,-1 0 0,-1 -1 0,-1 -1 2)),((0 0 2,-1 -1 2,1 0 2,0 0 2)),((0 0 4,0 0 2,1 0 2,0 0 4)),((0 1 4,0 0 4,1 0 4,0 1 4)),((1 1 4,0 1 4,1 0 4,1 1 4)),((1 1 2,1 1 4,1 0 4,1 1 2)),((2 2 2,1 1 2,2 -1 2,2 2 2)),((-1 2 2,2 2 2,0 2 0,-1 2 2)),((-1 2 2,0 2 0,-1 2 0,-1 2 2)),((-1 0 0,-1 2 2,-1 2 0,-1 0 0)),((-1 -1 2,-1 -1 0,1 -1 0,-1 -1 2)),((1 0 2,-1 -1 2,2 -1 2,1 0 2)),((0 0 4,1 0 2,1 0 4,0 0 4)),((1 1 2,1 0 4,1 0 2,1 1 2)),((2 -1 2,1 1 2,1 0 2,2 -1 2)),((2 2 2,2 -1 2,2 1 0,2 2 2)),((2 2 2,2 1 0,2 2 0,2 2 2)),((0 2 0,2 2 2,2 2 0,0 2 0)),((-1 -1 2,1 -1 0,2 -1 0,-1 -1 2)),((-1 -1 2,2 -1 0,2 -1 2,-1 -1 2)),((2 1 0,2 -1 2,2 -1 0,2 1 0)))");
797  lwfree(ewkt);
798  lwcollection_free(c);
799  lwgeom_free(g);
800 }
801 
802 static void test_lwmline_clip(void)
803 {
804  LWCOLLECTION *c;
805  char *ewkt;
806  LWGEOM *mline = NULL;
807  LWGEOM *line = NULL;
808 
809  /*
810  ** Set up the input line. Trivial one-member case.
811  */
812  mline = lwgeom_from_wkt("MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
813 
814  /* Clip in the middle, mid-range. */
815  c = lwgeom_clip_to_ordinate_range(mline, 'Y', 1.5, 2.5, 0);
816  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
817  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
818  lwfree(ewkt);
819  lwcollection_free(c);
820 
821  lwgeom_free(mline);
822 
823  /*
824  ** Set up the input line. Two-member case.
825  */
826  mline = 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);
827 
828  /* Clip off the top. */
829  c = lwgeom_clip_to_ordinate_range(mline, 'Y', 3.5, 5.5, 0);
830  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
831  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((1 3.5,1 4),(0 3.5,0 4))");
832  lwfree(ewkt);
833  lwcollection_free(c);
834 
835  lwgeom_free(mline);
836 
837  /*
838  ** Set up staggered input line to create multi-type output.
839  */
840  mline =
841  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);
842 
843  /* Clip from 0 upwards.. */
844  c = lwgeom_clip_to_ordinate_range(mline, 'Y', 0.0, 2.5, 0);
845  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
846  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 0),LINESTRING(0 0,0 1,0 2,0 2.5))");
847  lwfree(ewkt);
848  lwcollection_free(c);
849 
850  lwgeom_free(mline);
851 
852  /*
853  ** Set up input line from MAC
854  */
855  line = 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)",
856  LW_PARSER_CHECK_NONE);
857 
858  /* Clip from 3 to 3.5 */
859  c = lwgeom_clip_to_ordinate_range(line, 'Z', 3.0, 3.5, 0);
860  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
861  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))");
862  lwfree(ewkt);
863  lwcollection_free(c);
864 
865  /* Clip from 2 to 3.5 */
866  c = lwgeom_clip_to_ordinate_range(line, 'Z', 2.0, 3.5, 0);
867  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
868  ASSERT_STRING_EQUAL(ewkt,
869  "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))");
870  lwfree(ewkt);
871  lwcollection_free(c);
872 
873  /* Clip from 3 to 4 */
874  c = lwgeom_clip_to_ordinate_range(line, 'Z', 3.0, 4.0, 0);
875  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
876  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,4 4 4 4,3 3 3 5))");
877  lwfree(ewkt);
878  lwcollection_free(c);
879 
880  /* Clip from 2 to 3 */
881  c = lwgeom_clip_to_ordinate_range(line, 'Z', 2.0, 3.0, 0);
882  ewkt = lwgeom_to_ewkt((LWGEOM *)c);
883  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3),(3 3 3 5,2 2 2 6))");
884  lwfree(ewkt);
885  lwcollection_free(c);
886 
887  lwgeom_free(line);
888 }
889 
890 static void test_lwline_clip_big(void)
891 {
892  POINTARRAY *pa = ptarray_construct(1, 0, 3);
893  LWGEOM *line = (LWGEOM *)lwline_construct(SRID_UNKNOWN, NULL, pa);
894  LWCOLLECTION *c;
895  char *ewkt;
896  POINT4D p;
897 
898  p.x = 0.0;
899  p.y = 0.0;
900  p.z = 0.0;
901  ptarray_set_point4d(pa, 0, &p);
902 
903  p.x = 1.0;
904  p.y = 1.0;
905  p.z = 1.0;
906  ptarray_set_point4d(pa, 1, &p);
907 
908  p.x = 2.0;
909  p.y = 2.0;
910  p.z = 2.0;
911  ptarray_set_point4d(pa, 2, &p);
912 
913  c = lwgeom_clip_to_ordinate_range(line, 'Z', 0.5, 1.5, 0);
914  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
915  //printf("c = %s\n", ewkt);
916  ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0.5 0.5 0.5,1 1 1,1.5 1.5 1.5))");
917 
918  lwfree(ewkt);
919  lwcollection_free(c);
920  lwgeom_free(line);
921 }
922 
923 static void test_geohash_precision(void)
924 {
925  GBOX bbox;
926  GBOX bounds;
927  int precision = 0;
928  gbox_init(&bbox);
929  gbox_init(&bounds);
930 
931  bbox.xmin = 23.0;
932  bbox.xmax = 23.0;
933  bbox.ymin = 25.2;
934  bbox.ymax = 25.2;
935  precision = lwgeom_geohash_precision(bbox, &bounds);
936  //printf("\nprecision %d\n",precision);
937  CU_ASSERT_EQUAL(precision, 20);
938 
939  bbox.xmin = 23.0;
940  bbox.ymin = 23.0;
941  bbox.xmax = 23.1;
942  bbox.ymax = 23.1;
943  precision = lwgeom_geohash_precision(bbox, &bounds);
944  //printf("precision %d\n",precision);
945  CU_ASSERT_EQUAL(precision, 3);
946 
947  bbox.xmin = 23.0;
948  bbox.ymin = 23.0;
949  bbox.xmax = 23.0001;
950  bbox.ymax = 23.0001;
951  precision = lwgeom_geohash_precision(bbox, &bounds);
952  //printf("precision %d\n",precision);
953  CU_ASSERT_EQUAL(precision, 7);
954 
955 }
956 
957 static void test_geohash_point(void)
958 {
959  lwvarlena_t *geohash;
960 
961  geohash = geohash_point(0, 0, 16);
962  //printf("\ngeohash %s\n",geohash);
963  ASSERT_VARLENA_EQUAL(geohash, "s000000000000000");
964  lwfree(geohash);
965 
966  geohash = geohash_point(90, 0, 16);
967  //printf("\ngeohash %s\n",geohash);
968  ASSERT_VARLENA_EQUAL(geohash, "w000000000000000");
969  lwfree(geohash);
970 
971  geohash = geohash_point(20.012345, -20.012345, 15);
972  //printf("\ngeohash %s\n",geohash);
973  ASSERT_VARLENA_EQUAL(geohash, "kkqnpkue9ktbpe5");
974  lwfree(geohash);
975 
976 }
977 
978 static void test_geohash(void)
979 {
980  LWPOINT *lwpoint = NULL;
981  LWLINE *lwline = NULL;
982  LWMLINE *lwmline = NULL;
983  lwvarlena_t *geohash = NULL;
984 
985  lwpoint = (LWPOINT*)lwgeom_from_wkt("POINT(23.0 25.2)", LW_PARSER_CHECK_NONE);
986  geohash = lwgeom_geohash((LWGEOM*)lwpoint,0);
987  //printf("\ngeohash %s\n",geohash);
988  ASSERT_VARLENA_EQUAL(geohash, "ss2r77s0du7p2ewb8hmx");
989  lwpoint_free(lwpoint);
990  lwfree(geohash);
991 
992  lwpoint = (LWPOINT*)lwgeom_from_wkt("POINT(23.0 25.2 2.0)", LW_PARSER_CHECK_NONE);
993  geohash = lwgeom_geohash((LWGEOM*)lwpoint,0);
994  //printf("geohash %s\n",geohash);
995  ASSERT_VARLENA_EQUAL(geohash, "ss2r77s0du7p2ewb8hmx");
996  lwpoint_free(lwpoint);
997  lwfree(geohash);
998 
999  lwline = (LWLINE*)lwgeom_from_wkt("LINESTRING(23.0 23.0,23.1 23.1)", LW_PARSER_CHECK_NONE);
1000  geohash = lwgeom_geohash((LWGEOM*)lwline,0);
1001  //printf("geohash %s\n",geohash);
1002  ASSERT_VARLENA_EQUAL(geohash, "ss0");
1003  lwline_free(lwline);
1004  lwfree(geohash);
1005 
1006  lwline = (LWLINE*)lwgeom_from_wkt("LINESTRING(23.0 23.0,23.001 23.001)", LW_PARSER_CHECK_NONE);
1007  geohash = lwgeom_geohash((LWGEOM*)lwline,0);
1008  //printf("geohash %s\n",geohash);
1009  ASSERT_VARLENA_EQUAL(geohash, "ss06g7h");
1010  lwline_free(lwline);
1011  lwfree(geohash);
1012 
1013  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);
1014  geohash = lwgeom_geohash((LWGEOM*)lwmline,0);
1015  //printf("geohash %s\n",geohash);
1016  ASSERT_VARLENA_EQUAL(geohash, "ss0");
1017  lwmline_free(lwmline);
1018  lwfree(geohash);
1019 }
1020 
1021 static void test_isclosed(void)
1022 {
1023  LWGEOM *geom;
1024 
1025  /* LINESTRING */
1026 
1027  /* Not Closed on 2D */
1028  geom = lwgeom_from_wkt("LINESTRING(1 2,3 4)", LW_PARSER_CHECK_NONE);
1029  CU_ASSERT(!lwline_is_closed((LWLINE *) geom));
1030  lwgeom_free(geom);
1031 
1032  /* Closed on 2D */
1033  geom = lwgeom_from_wkt("LINESTRING(1 2,3 4,1 2)", LW_PARSER_CHECK_NONE);
1034  CU_ASSERT(lwline_is_closed((LWLINE *) geom));
1035  lwgeom_free(geom);
1036 
1037  /* Not closed on 3D */
1038  geom = lwgeom_from_wkt("LINESTRING(1 2 3,4 5 6)", LW_PARSER_CHECK_NONE);
1039  CU_ASSERT(!lwline_is_closed((LWLINE *) geom));
1040  lwgeom_free(geom);
1041 
1042  /* Closed on 3D */
1043  geom = lwgeom_from_wkt("LINESTRING(1 2 3,4 5 6,1 2 3)", LW_PARSER_CHECK_NONE);
1044  CU_ASSERT(lwline_is_closed((LWLINE *) geom));
1045  lwgeom_free(geom);
1046 
1047  /* Closed on 4D, even if M is not the same */
1048  geom = lwgeom_from_wkt("LINESTRING(1 2 3 4,5 6 7 8,1 2 3 0)", LW_PARSER_CHECK_NONE);
1049  CU_ASSERT(lwline_is_closed((LWLINE *) geom));
1050  lwgeom_free(geom);
1051 
1052 
1053  /* CIRCULARSTRING */
1054 
1055  /* Not Closed on 2D */
1056  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2,3 4,5 6)", LW_PARSER_CHECK_NONE);
1057  CU_ASSERT(!lwcircstring_is_closed((LWCIRCSTRING *) geom));
1058  lwgeom_free(geom);
1059 
1060  /* Closed on 2D */
1061  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2,3 4,1 2)", LW_PARSER_CHECK_NONE);
1062  CU_ASSERT(lwcircstring_is_closed((LWCIRCSTRING *) geom));
1063  lwgeom_free(geom);
1064 
1065  /* Not closed on 3D */
1066  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2 3,4 5 6,7 8 9)", LW_PARSER_CHECK_NONE);
1067  CU_ASSERT(!lwcircstring_is_closed((LWCIRCSTRING *) geom));
1068  lwgeom_free(geom);
1069 
1070  /* Closed on 3D */
1071  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2 3,4 5 6,1 2 3)", LW_PARSER_CHECK_NONE);
1072  CU_ASSERT(lwcircstring_is_closed((LWCIRCSTRING *) geom));
1073  lwgeom_free(geom);
1074 
1075  /* Closed on 4D, even if M is not the same */
1076  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2 3 4,5 6 7 8,1 2 3 0)", LW_PARSER_CHECK_NONE);
1077  CU_ASSERT(lwcircstring_is_closed((LWCIRCSTRING *) geom));
1078  lwgeom_free(geom);
1079 
1080 
1081  /* COMPOUNDCURVE */
1082 
1083  /* Not Closed on 2D */
1084  geom = lwgeom_from_wkt("COMPOUNDCURVE(CIRCULARSTRING(1 2,3 4,1 2),(1 2,7 8,5 6))", LW_PARSER_CHECK_NONE);
1085  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
1086  lwgeom_free(geom);
1087 
1088  geom = lwgeom_from_wkt("COMPOUNDCURVE((1 2,3 4,1 2),CIRCULARSTRING(1 2,7 8,5 6))", LW_PARSER_CHECK_NONE);
1089  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
1090  lwgeom_free(geom);
1091 
1092  /* Closed on 2D */
1093  geom = lwgeom_from_wkt("COMPOUNDCURVE(CIRCULARSTRING(1 2,3 4,5 6), (5 6,7 8,1 2))", LW_PARSER_CHECK_NONE);
1094  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
1095  lwgeom_free(geom);
1096 
1097  geom = lwgeom_from_wkt("COMPOUNDCURVE((1 2,3 4,5 6),CIRCULARSTRING(5 6,7 8,1 2))", LW_PARSER_CHECK_NONE);
1098  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
1099  lwgeom_free(geom);
1100 
1101  /* Not Closed on 3D */
1102  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);
1103  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
1104  lwgeom_free(geom);
1105 
1106  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);
1107  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
1108  lwgeom_free(geom);
1109 
1110  /* Closed on 3D */
1111  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);
1112  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
1113  lwgeom_free(geom);
1114 
1115  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);
1116  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
1117  lwgeom_free(geom);
1118 
1119  /* Closed on 4D, even if M is not the same */
1120  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);
1121  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
1122  lwgeom_free(geom);
1123 }
1124 
1125 
1126 static void test_geohash_point_as_int(void)
1127 {
1128  unsigned int gh;
1129  POINT2D p;
1130  unsigned long long rs;
1131 
1132  p.x = 50; p.y = 35;
1133  gh = geohash_point_as_int(&p);
1134  rs = 3440103613;
1135  CU_ASSERT_EQUAL(gh, rs);
1136  p.x = 140; p.y = 45;
1137  gh = geohash_point_as_int(&p);
1138  rs = 3982480893;
1139  CU_ASSERT_EQUAL(gh, rs);
1140  p.x = 140; p.y = 55;
1141  gh = geohash_point_as_int(&p);
1142  rs = 4166944232;
1143  CU_ASSERT_EQUAL(gh, rs);
1144 }
1145 
1146 static void
1147 test_geohash_bbox(void)
1148 {
1149  double lat[2], lon[2];
1150 
1151  /* SELECT ST_GeoHash(ST_SetSRID(ST_MakePoint(-126,48),4326)) */
1152  decode_geohash_bbox("c0w3hf1s70w3hf1s70w3", lat, lon, 100);
1153  CU_ASSERT_DOUBLE_EQUAL(lat[0], 48, 1e-11);
1154  CU_ASSERT_DOUBLE_EQUAL(lat[1], 48, 1e-11);
1155  CU_ASSERT_DOUBLE_EQUAL(lon[0], -126, 1e-11);
1156  CU_ASSERT_DOUBLE_EQUAL(lon[1], -126, 1e-11);
1157 
1158  cu_error_msg_reset();
1159  decode_geohash_bbox("@@@@@@", lat, lon, 100);
1160  ASSERT_STRING_EQUAL(cu_error_msg, "decode_geohash_bbox: Invalid character '@'");
1161 }
1162 
1163 static void test_lwgeom_remove_repeated_points(void)
1164 {
1165  LWGEOM *g;
1166  char *ewkt;
1167  int modified = LW_FALSE;
1168 
1169  g = lwgeom_from_wkt("MULTIPOINT(0 0, 10 0, 10 10, 10 10, 0 10, 0 10, 0 10, 0 0, 0 0, 0 0, 5 5, 0 0, 5 5)", LW_PARSER_CHECK_NONE);
1170  modified = lwgeom_remove_repeated_points_in_place(g, 1);
1171  ASSERT_INT_EQUAL(modified, LW_TRUE);
1172  ewkt = lwgeom_to_ewkt(g);
1173  ASSERT_STRING_EQUAL(ewkt, "MULTIPOINT(0 0,10 0,10 10,0 10,5 5)");
1174  lwgeom_free(g);
1175  lwfree(ewkt);
1176 
1177  g = lwgeom_from_wkt("LINESTRING(1612830.15445 4841287.12672,1612830.15824 4841287.12674,1612829.98813 4841274.56198)", LW_PARSER_CHECK_NONE);
1178  modified = lwgeom_remove_repeated_points_in_place(g, 0.01);
1179  ASSERT_INT_EQUAL(modified, LW_TRUE);
1180  ewkt = lwgeom_to_ewkt(g);
1181  ASSERT_STRING_EQUAL(ewkt, "LINESTRING(1612830.15445 4841287.12672,1612829.98813 4841274.56198)");
1182  lwgeom_free(g);
1183  lwfree(ewkt);
1184 
1185  g = lwgeom_from_wkt("MULTIPOINT(0 0,10 0,10 10,10 10,0 10,0 10,0 10,0 0,0 0,0 0,5 5,5 5,5 8,8 8,8 8,8 8,8 5,8 5,5 5,5 5,5 5,5 5,5 5,50 50,50 50,50 50,50 60,50 60,50 60,60 60,60 50,60 50,50 50,55 55,55 58,58 58,58 55,58 55,55 55)", LW_PARSER_CHECK_NONE);
1186  modified = lwgeom_remove_repeated_points_in_place(g, 1);
1187  ASSERT_INT_EQUAL(modified, LW_TRUE);
1188  ewkt = lwgeom_to_ewkt(g);
1189  ASSERT_STRING_EQUAL(
1190  ewkt, "MULTIPOINT(0 0,10 0,10 10,0 10,5 5,5 8,8 8,8 5,50 50,50 60,60 60,60 50,55 55,55 58,58 58,58 55)");
1191  lwgeom_free(g);
1192  lwfree(ewkt);
1193 
1194  g = lwgeom_from_wkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", LW_PARSER_CHECK_NONE);
1195  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1196  ASSERT_INT_EQUAL(modified, LW_TRUE);
1197  ewkt = lwgeom_to_ewkt(g);
1198  ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 1,1 0,0 0))");
1199  lwgeom_free(g);
1200  lwfree(ewkt);
1201 
1202  // Test the return value (modified or not)
1203  g = lwgeom_from_wkt("POINT(0 0)", LW_PARSER_CHECK_NONE);
1204  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1205  ASSERT_INT_EQUAL(modified, LW_FALSE);
1206  lwgeom_free(g);
1207 
1208  g = lwgeom_from_wkt("TRIANGLE((0 0, 5 0, 3 3, 0 0))", LW_PARSER_CHECK_NONE);
1209  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1210  ASSERT_INT_EQUAL(modified, LW_FALSE);
1211  lwgeom_free(g);
1212 
1213  g = lwgeom_from_wkt("POLYGON((0 0,0 1,1 1,1 0,0 0))", LW_PARSER_CHECK_NONE);
1214  modified = lwgeom_remove_repeated_points_in_place(g, 0.1);
1215  ASSERT_INT_EQUAL(modified, LW_FALSE);
1216  ewkt = lwgeom_to_ewkt(g);
1217  ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,0 1,1 1,1 0,0 0))");
1218  lwgeom_free(g);
1219  lwfree(ewkt);
1220 
1221  g = lwgeom_from_wkt("POLYGON((0 0,0 1,1 1,1 0,0 0), (0.4 0.4, 0.4 0.4, 0.4 0.5, 0.5 0.5, 0.5 0.4, 0.4 0.4))",
1222  LW_PARSER_CHECK_NONE);
1223  modified = lwgeom_remove_repeated_points_in_place(g, 0.1);
1224  ASSERT_INT_EQUAL(modified, LW_TRUE);
1225  ewkt = lwgeom_to_ewkt(g);
1226  ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,0 1,1 1,1 0,0 0),(0.4 0.4,0.5 0.5,0.5 0.4,0.4 0.4))");
1227  lwgeom_free(g);
1228  lwfree(ewkt);
1229 
1230  g = lwgeom_from_wkt("GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0,1 0,1 1,0 1,0 0)))", LW_PARSER_CHECK_NONE);
1231  modified = lwgeom_remove_repeated_points_in_place(g, 0.1);
1232  ASSERT_INT_EQUAL(modified, LW_FALSE);
1233  ewkt = lwgeom_to_ewkt(g);
1234  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0,1 0,1 1,0 1,0 0)))");
1235  lwgeom_free(g);
1236  lwfree(ewkt);
1237 
1238  g = lwgeom_from_wkt("GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0, 0 1, 1 1, 1 0, 0 0)))", LW_PARSER_CHECK_NONE);
1239  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1240  ASSERT_INT_EQUAL(modified, LW_TRUE);
1241  ewkt = lwgeom_to_ewkt(g);
1242  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0,1 1,1 0,0 0)))");
1243  lwgeom_free(g);
1244  lwfree(ewkt);
1245 
1246  g = lwgeom_from_wkt("GEOMETRYCOLLECTION(POLYGON((0 0, 0 1, 1 1, 1 0, 0 0)),POINT(2 0))", LW_PARSER_CHECK_NONE);
1247  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1248  ASSERT_INT_EQUAL(modified, LW_TRUE);
1249  ewkt = lwgeom_to_ewkt(g);
1250  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POLYGON((0 0,1 1,1 0,0 0)),POINT(2 0))");
1251  lwgeom_free(g);
1252  lwfree(ewkt);
1253 }
1254 
1255 static void test_lwgeom_simplify(void)
1256 {
1257  LWGEOM *l;
1258  LWGEOM *g;
1259  char *ewkt;
1260 
1261  /* Simplify but only so far... */
1262  g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
1263  l = lwgeom_simplify(g, 10, LW_TRUE);
1264  ewkt = lwgeom_to_ewkt(l);
1265  ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)");
1266  lwgeom_free(g);
1267  lwgeom_free(l);
1268  lwfree(ewkt);
1269 
1270  /* Simplify but only so far... */
1271  g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
1272  l = lwgeom_simplify(g, 10, LW_TRUE);
1273  ewkt = lwgeom_to_ewkt(l);
1274  ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))");
1275  lwgeom_free(g);
1276  lwgeom_free(l);
1277  lwfree(ewkt);
1278 
1279  /* Simplify and collapse */
1280  g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
1281  l = lwgeom_simplify(g, 10, LW_FALSE);
1282  CU_ASSERT_EQUAL(l, NULL);
1283  lwgeom_free(g);
1284  lwgeom_free(l);
1285 
1286  /* Simplify and collapse */
1287  g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
1288  l = lwgeom_simplify(g, 10, LW_FALSE);
1289  CU_ASSERT_EQUAL(l, NULL);
1290  lwgeom_free(g);
1291  lwgeom_free(l);
1292 
1293  /* Not simplifiable */
1294  g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE);
1295  l = lwgeom_simplify(g, 1.0, LW_FALSE);
1296  ewkt = lwgeom_to_ewkt(l);
1297  ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)");
1298  lwgeom_free(g);
1299  lwgeom_free(l);
1300  lwfree(ewkt);
1301 
1302  /* Simplifiable */
1303  g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE);
1304  l = lwgeom_simplify(g, 1.0, LW_FALSE);
1305  ewkt = lwgeom_to_ewkt(l);
1306  ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)");
1307  lwgeom_free(g);
1308  lwgeom_free(l);
1309  lwfree(ewkt);
1310 
1311  /* POLYGON with multiple inner rings*/
1312  g = lwgeom_from_wkt(
1313  "POLYGON("
1314  "(0 0, 100 0, 100 100, 0 100, 0 0),"
1315  "(1 1, 1 5, 5 5, 5 1, 1 1),"
1316  "(20 20, 20 40, 40 40, 40 20, 20 20)"
1317  ")",
1318  LW_PARSER_CHECK_NONE);
1319  l = lwgeom_simplify(g, 10, LW_FALSE);
1320  ewkt = lwgeom_to_ewkt(l);
1321  ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,100 0,100 100,0 100,0 0),(20 20,20 40,40 40,40 20,20 20))");
1322  lwgeom_free(g);
1323  lwgeom_free(l);
1324  lwfree(ewkt);
1325 
1326  /* Reorder inner rings: Same result */
1327  g = lwgeom_from_wkt(
1328  "POLYGON("
1329  "(0 0, 100 0, 100 100, 0 100, 0 0),"
1330  "(20 20, 20 40, 40 40, 40 20, 20 20),"
1331  "(1 1, 1 5, 5 5, 5 1, 1 1)"
1332  ")",
1333  LW_PARSER_CHECK_NONE);
1334  l = lwgeom_simplify(g, 10, LW_FALSE);
1335  ewkt = lwgeom_to_ewkt(l);
1336  ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,100 0,100 100,0 100,0 0),(20 20,20 40,40 40,40 20,20 20))");
1337  lwgeom_free(g);
1338  lwgeom_free(l);
1339  lwfree(ewkt);
1340 
1341  g = lwgeom_from_wkt(
1342  "POLYGON("
1343  "(0 0, 100 0, 100 100, 0 100, 0 0),"
1344  "(20 20, 20 40, 40 40, 40 20, 20 20),"
1345  "(1 1, 1 5, 5 5, 5 1, 1 1)"
1346  ")",
1347  LW_PARSER_CHECK_NONE);
1348  l = lwgeom_simplify(g, 100, LW_FALSE);
1349  CU_ASSERT_EQUAL(l, NULL);
1350  lwgeom_free(g);
1351  lwgeom_free(l);
1352 }
1353 
1354 
1355 static void do_median_dims_check(char* wkt, int expected_dims)
1356 {
1357  LWGEOM* g = lwgeom_from_wkt(wkt, LW_PARSER_CHECK_NONE);
1358  LWPOINT* result = lwgeom_median(g, 1e-8, 100, LW_FALSE);
1359 
1360  CU_ASSERT_EQUAL(expected_dims, lwgeom_ndims((LWGEOM*) result));
1361 
1362  lwgeom_free(g);
1363  lwpoint_free(result);
1364 }
1365 
1366 static void test_median_handles_3d_correctly(void)
1367 {
1368  do_median_dims_check("MULTIPOINT ((1 3), (4 7), (2 9), (0 4), (2 2))", 2);
1369  do_median_dims_check("MULTIPOINT Z ((1 3 4), (4 7 8), (2 9 1), (0 4 4), (2 2 3))", 3);
1370  do_median_dims_check("MULTIPOINT M ((1 3 4), (4 7 8), (2 9 1), (0 4 4), (2 2 3))", 2);
1371  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);
1372 }
1373 
1374 static double
1375 test_weighted_distance(const POINT4D* curr, const POINT4D* points, uint32_t npoints)
1376 {
1377  double distance = 0.0;
1378  uint32_t i;
1379  for (i = 0; i < npoints; i++)
1380  {
1381  distance += distance3d_pt_pt((POINT3D*)curr, (POINT3D*)&points[i]) * points[i].m;
1382  }
1383 
1384  return distance;
1385 }
1386 
1387 static void do_median_test(char* input, char* expected, int fail_if_not_converged, int iter_count)
1388 {
1389  cu_error_msg_reset();
1390  LWGEOM* g = lwgeom_from_wkt(input, LW_PARSER_CHECK_NONE);
1391  LWPOINT* expected_result = NULL;
1392  POINT4D actual_pt;
1393  POINT4D expected_pt;
1394  const double tolerance = FP_TOLERANCE / 10.0;
1395 
1396  LWPOINT* result = lwgeom_median(g, tolerance, iter_count, fail_if_not_converged);
1397  int passed = LW_FALSE;
1398 
1399  if (expected != NULL)
1400  {
1401  expected_result = lwgeom_as_lwpoint(lwgeom_from_wkt(expected, LW_PARSER_CHECK_NONE));
1402  lwpoint_getPoint4d_p(expected_result, &expected_pt);
1403  }
1404  if (result != NULL)
1405  {
1406  lwpoint_getPoint4d_p(result, &actual_pt);
1407  }
1408 
1409  if (result != NULL && expected != NULL) /* got something, expecting something */
1410  {
1411  passed = LW_TRUE;
1412  passed = passed && lwgeom_is_empty((LWGEOM*) expected_result) == lwgeom_is_empty((LWGEOM*) result);
1413  passed = passed && (lwgeom_has_z((LWGEOM*) expected_result) == lwgeom_has_z((LWGEOM*) result));
1414 
1415  if (passed && !lwgeom_is_empty((LWGEOM*) result))
1416  {
1417  if (g->type == POINTTYPE)
1418  {
1419  passed &= fabs(actual_pt.x - expected_pt.x) < tolerance;
1420  passed &= fabs(actual_pt.y - expected_pt.y) < tolerance;
1421  passed &= (!lwgeom_has_z((LWGEOM*) expected_result) || fabs(actual_pt.z - expected_pt.z) < tolerance);
1422  passed &= (!lwgeom_has_m((LWGEOM*) expected_result) || fabs(actual_pt.m - expected_pt.m) < tolerance);
1423  }
1424  else
1425  {
1426  /* Check that the difference between the obtained geometric
1427  median and the expected point is within tolerance */
1428  uint32_t npoints = 1;
1429  int input_empty = LW_TRUE;
1430  POINT4D* points = lwmpoint_extract_points_4d(lwgeom_as_lwmpoint(g), &npoints, &input_empty);
1431  double distance_expected = test_weighted_distance(&expected_pt, points, npoints);
1432  double distance_result = test_weighted_distance(&actual_pt, points, npoints);
1433 
1434  passed = distance_result <= (1.0 + tolerance) * distance_expected;
1435  if (!passed)
1436  {
1437  printf("Diff: Got %.10f Expected %.10f\n", distance_result, distance_expected);
1438  }
1439  lwfree(points);
1440  }
1441  }
1442 
1443  if (!passed)
1444  {
1445  printf("median_test input %s (parsed %s) expected %s got %s\n",
1446  input, lwgeom_to_ewkt(g),
1447  lwgeom_to_ewkt((LWGEOM*) expected_result),
1448  lwgeom_to_ewkt((LWGEOM*) result));
1449  }
1450 
1451  }
1452  else if (result == NULL && expected == NULL) /* got nothing, expecting nothing */
1453  {
1454  passed = LW_TRUE;
1455  }
1456  else if (result != NULL && expected == NULL) /* got something, expecting nothing */
1457  {
1458  passed = LW_FALSE;
1459  printf("median_test input %s (parsed %s) expected NULL got %s\n", input, lwgeom_to_ewkt(g), lwgeom_to_ewkt((LWGEOM*) result));
1460  }
1461  else if (result == NULL && expected != NULL) /* got nothing, expecting something */
1462  {
1463  passed = LW_FALSE;
1464  printf("%s", cu_error_msg);
1465  printf("median_test input %s (parsed %s) expected %s got NULL\n", input, lwgeom_to_ewkt(g), lwgeom_to_ewkt((LWGEOM*) expected_result));
1466  }
1467 
1468  CU_ASSERT_TRUE(passed);
1469 
1470  lwgeom_free(g);
1471  lwpoint_free(expected_result);
1472  lwpoint_free(result);
1473 }
1474 
1475 static void test_median_robustness(void)
1476 {
1477  /* A simple implementation of Weiszfeld's algorithm will fail if the median is equal
1478  * to any one of the inputs, during any iteration of the algorithm.
1479  *
1480  * Because the algorithm uses the centroid as a starting point, this situation will
1481  * occur in the test case below.
1482  */
1483  do_median_test("MULTIPOINT ((0 -1), (0 0), (0 1))", "POINT (0 0)", LW_TRUE, 1000);
1484 
1485  /* Same as above but 3D, and shifter */
1486  do_median_test("MULTIPOINT ((1 -1 3), (1 0 2), (2 1 1))", "POINT (1 0 2)", LW_TRUE, 1000);
1487 
1488  /* Starting point is duplicated */
1489  do_median_test("MULTIPOINT ((0 -1), (0 0), (0 0), (0 1))", "POINT (0 0)", LW_TRUE, 1000);
1490 
1491  /* Cube */
1492  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))",
1493  "POINT (15 15 15)", LW_TRUE, 1000);
1494 
1495  /* Some edge cases */
1496  do_median_test("POINT (7 6)", "POINT (7 6)", LW_TRUE, 1000);
1497  do_median_test("POINT (7 6 2)", "POINT (7 6 2)", LW_TRUE, 1000);
1498  do_median_test("MULTIPOINT ((7 6 2), EMPTY)", "POINT (7 6 2)", LW_TRUE, 1000);
1499 
1500  /* Empty input */
1501  do_median_test("MULTIPOINT EMPTY", "POINT EMPTY", LW_FALSE, 1000);
1502  do_median_test("MULTIPOINT (EMPTY)", "POINT EMPTY", LW_FALSE, 1000);
1503  do_median_test("MULTIPOINT EMPTY", "POINT EMPTY", LW_TRUE, 1000);
1504  do_median_test("MULTIPOINT (EMPTY)", "POINT EMPTY", LW_TRUE, 1000);
1505  do_median_test("MULTIPOINT ZM (1 -1 3 1, 1 0 2 7, 2 1 1 1, EMPTY)", "POINT (1 0 2)", LW_TRUE, 1000);
1506 
1507  /* Weighted input */
1508  do_median_test("MULTIPOINT ZM (1 -1 3 1, 1 0 2 7, 2 1 1 1)", "POINT (1 0 2)", LW_TRUE, 1000);
1509  do_median_test("MULTIPOINT ZM (-1 1 -3 1, -1 0 -2 7, -2 -1 -1 1)", "POINT (-1 0 -2)", LW_TRUE, 1000);
1510  do_median_test("MULTIPOINT ZM (-1 1 -3 1, -1 0 -2 7, -2 -1 -1 0.5, -2 -1 -1 0.5)", "POINT (-1 0 -2)", LW_TRUE, 1000);
1511 
1512  /* Point that is replaced by two half-weighted */
1513  do_median_test("MULTIPOINT ZM ((0 -1 0 1), (0 0 0 1), (0 1 0 0.5), (0 1 0 0.5))", "POINT (0 0 0)", LW_TRUE, 1000);
1514  /* Point is doubled and then erased by negative weight */
1515  do_median_test("MULTIPOINT ZM ((1 -1 3 1), (1 0 2 7), (2 1 1 2), (2 1 1 -1))", NULL, LW_TRUE, 1000);
1516  do_median_test("MULTIPOINT ZM ((1 -1 3 1), (1 0 2 7), (2 1 1 2), (2 1 1 -1))", NULL, LW_FALSE, 1000);
1517  /* Weightless input won't converge */
1518  do_median_test("MULTIPOINT ZM ((0 -1 0 0), (0 0 0 0), (0 0 0 0), (0 1 0 0))", NULL, LW_FALSE, 1000);
1519  do_median_test("MULTIPOINT ZM ((0 -1 0 0), (0 0 0 0), (0 0 0 0), (0 1 0 0))", NULL, LW_TRUE, 1000);
1520  /* Negative weight won't converge */
1521  do_median_test("MULTIPOINT ZM ((0 -1 0 -1), (0 0 0 -1), (0 1 0 -1))", NULL, LW_FALSE, 1000);
1522  do_median_test("MULTIPOINT ZM ((0 -1 0 -1), (0 0 0 -1), (0 1 0 -1))", NULL, LW_TRUE, 1000);
1523 
1524  /* Bind convergence too tightly */
1525  do_median_test("MULTIPOINT ((0 0), (1 1), (0 1), (2 2))", "POINT(0.75 1.0)", LW_FALSE, 0);
1526  do_median_test("MULTIPOINT ((0 0), (1 1), (0 1), (2 2))", NULL, LW_TRUE, 1);
1527  /* Unsupported geometry type */
1528  do_median_test("POLYGON((1 0,0 1,1 2,2 1,1 0))", NULL, LW_TRUE, 1000);
1529  do_median_test("POLYGON((1 0,0 1,1 2,2 1,1 0))", NULL, LW_FALSE, 1000);
1530 
1531  /* Median point is included */
1532  do_median_test("MULTIPOINT ZM ("
1533  "(1480 0 200 100),"
1534  "(620 0 200 100),"
1535  "(1000 0 -200 100),"
1536  "(1000 0 -590 100),"
1537  "(1025 0 65 100),"
1538  "(1025 0 -65 100)"
1539  ")",
1540  "POINT (1025 0 -65)", LW_TRUE, 10000);
1541 
1542 #if 0
1543  /* Leads to invalid result (0 0 0) with 80bit (fmulp + faddp) precision. ok with 64 bit float ops */
1544  do_median_test("MULTIPOINT ZM ("
1545  "(0 0 20000 0.5),"
1546  "(0 0 59000 0.5),"
1547  "(0 -3000 -3472.22222222222262644208967685699462890625 1),"
1548  "(0 3000 3472.22222222222262644208967685699462890625 1),"
1549  "(0 0 -1644.736842105263121993630193173885345458984375 1),"
1550  "(0 0 1644.736842105263121993630193173885345458984375 1),"
1551  "(0 48000 -20000 1.3),"
1552  "(0 -48000 -20000 1.3)"
1553  ")",
1554  "POINT (0 0 0)", LW_TRUE, 10000);
1555 #endif
1556 
1557 #if 0
1558  /* Leads to invalid result (0 0 0) with 64bit (vfmadd231sd) precision. Ok with 80 bit float ops */
1559  do_median_test("MULTIPOINT ZM ("
1560  "(0 0 20000 0.5),"
1561  "(0 0 59000 0.5),"
1562  "(0 -3000 -3472.22222222222262644208967685699462890625 1),"
1563  "(0 3000 3472.22222222222262644208967685699462890625 1),"
1564  "(0 -0.00000000000028047739569477638384522295466033823196 -1644.736842105263121993630193173885345458984375 1),"
1565  "(0 0.00000000000028047739569477638384522295466033823196 1644.736842105263121993630193173885345458984375 1),"
1566  "(0 48000 -20000 1.3),"
1567  "(0 -48000 -20000 1.3)"
1568  ")",
1569  "POINT (0 0 0)", LW_TRUE, 10000);
1570 #endif
1571 }
1572 
1573 static void test_point_density(void)
1574 {
1575  LWGEOM *geom;
1576  LWMPOINT *mpt;
1577  LWMPOINT *mpt2;
1578  LWPOINT *pt;
1579  LWPOINT *pt2;
1580  int eq, i;
1581  // char *ewkt;
1582 
1583  /* POLYGON */
1584  geom = lwgeom_from_wkt("POLYGON((0 0,1 0,1 1,0 1,0 0))", LW_PARSER_CHECK_NONE);
1585  mpt = lwgeom_to_points(geom, 100, 0); /* Set a zero seed to base it on Unix time and process ID */
1586  CU_ASSERT_EQUAL(mpt->ngeoms,100);
1587 
1588  /* Run a second time with a zero seed to get a different multipoint sequence */
1589  mpt2 = lwgeom_to_points(geom, 100, 0);
1590  eq = 0;
1591  for (i = 0; i < 100; i++)
1592  {
1593  pt = (LWPOINT*)mpt->geoms[i];
1594  pt2 = (LWPOINT*)mpt2->geoms[i];
1595  if (lwpoint_get_x(pt) == lwpoint_get_x(pt2) && lwpoint_get_y(pt) == lwpoint_get_y(pt2))
1596  eq++;
1597  }
1598  CU_ASSERT_EQUAL(eq, 0);
1599  lwmpoint_free(mpt);
1600  lwmpoint_free(mpt2);
1601  pt = NULL;
1602  pt2 = NULL;
1603 
1604  /* Set seed to get a deterministic sequence */
1605  mpt = lwgeom_to_points(geom, 1000, 12345);
1606 
1607  /* Check to find a different multipoint sequence with different seed */
1608  mpt2 = lwgeom_to_points(geom, 1000, 54321);
1609  eq = 0;
1610  for (i = 0; i < 1000; i++)
1611  {
1612  pt = (LWPOINT*)mpt->geoms[i];
1613  pt2 = (LWPOINT*)mpt2->geoms[i];
1614  if (lwpoint_get_x(pt) == lwpoint_get_x(pt2) && lwpoint_get_y(pt) == lwpoint_get_y(pt2))
1615  eq++;
1616  }
1617  CU_ASSERT_EQUAL(eq, 0);
1618  lwmpoint_free(mpt2);
1619  pt = NULL;
1620  pt2 = NULL;
1621 
1622  /* Check to find an identical multipoint sequence with same seed */
1623  mpt2 = lwgeom_to_points(geom, 1000, 12345);
1624  eq = 0;
1625  for (i = 0; i < 1000; i++)
1626  {
1627  pt = (LWPOINT*)mpt->geoms[i];
1628  pt2 = (LWPOINT*)mpt2->geoms[i];
1629  if (lwpoint_get_x(pt) == lwpoint_get_x(pt2) && lwpoint_get_y(pt) == lwpoint_get_y(pt2))
1630  eq++;
1631  }
1632  CU_ASSERT_EQUAL(eq, 1000);
1633  lwmpoint_free(mpt2);
1634  pt = NULL;
1635  pt2 = NULL;
1636 
1637 
1638  /* Check if the 1000th point is the expected value.
1639  * Note that if the RNG is not portable, this test may fail. */
1640  pt = (LWPOINT*)mpt->geoms[999];
1641  // ewkt = lwgeom_to_ewkt((LWGEOM*)pt);
1642  // printf("%s\n", ewkt);
1643  // lwfree(ewkt);
1644  CU_ASSERT_DOUBLE_EQUAL(lwpoint_get_x(pt), 0.801167838758, 1e-11);
1645  CU_ASSERT_DOUBLE_EQUAL(lwpoint_get_y(pt), 0.345281131175, 1e-11);
1646  lwmpoint_free(mpt);
1647  pt = NULL;
1648 
1649  mpt = lwgeom_to_points(geom, 0, 0);
1650  CU_ASSERT_EQUAL(mpt, NULL);
1651  lwmpoint_free(mpt);
1652 
1653  lwgeom_free(geom);
1654 
1655  /* MULTIPOLYGON */
1656  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);
1657 
1658  mpt = lwgeom_to_points(geom, 1000, 0);
1659  CU_ASSERT_EQUAL(mpt->ngeoms,1000);
1660  lwmpoint_free(mpt);
1661 
1662  mpt = lwgeom_to_points(geom, 1, 0);
1663  CU_ASSERT_EQUAL(mpt->ngeoms,1);
1664  lwmpoint_free(mpt);
1665 
1666  lwgeom_free(geom);
1667 }
1668 
1669 static void test_lwpoly_construct_circle(void)
1670 {
1671  LWPOLY* p;
1672  const GBOX* g;
1673  const int32_t srid = 4326;
1674  const uint32_t segments_per_quad = 17;
1675  const int x = 10;
1676  const int y = 20;
1677  const int r = 5;
1678 
1679  /* With normal arguments you should get something circle-y */
1680  p = lwpoly_construct_circle(srid, x, y, r, segments_per_quad, LW_TRUE);
1681 
1682  ASSERT_INT_EQUAL(lwgeom_count_vertices(lwpoly_as_lwgeom(p)), segments_per_quad * 4 + 1);
1683  ASSERT_INT_EQUAL(lwgeom_get_srid(lwpoly_as_lwgeom(p)), srid)
1684 
1685  g = lwgeom_get_bbox(lwpoly_as_lwgeom(p));
1686  CU_ASSERT_DOUBLE_EQUAL(g->xmin, x-r, 0.1);
1687  CU_ASSERT_DOUBLE_EQUAL(g->xmax, x+r, 0.1);
1688  CU_ASSERT_DOUBLE_EQUAL(g->ymin, y-r, 0.1);
1689  CU_ASSERT_DOUBLE_EQUAL(g->ymax, y+r, 0.1);
1690 
1691  CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(lwpoly_as_lwgeom(p)), M_PI*5*5, 0.1);
1692 
1693  lwpoly_free(p);
1694 
1695  /* No segments? No circle. */
1696  p = lwpoly_construct_circle(srid, x, y, r, 0, LW_TRUE);
1697  CU_ASSERT_TRUE(p == NULL);
1698 
1699  /* Negative radius? No circle. */
1700  p = lwpoly_construct_circle(srid, x, y, -1e-3, segments_per_quad, LW_TRUE);
1701  CU_ASSERT_TRUE(p == NULL);
1702 
1703  /* Zero radius? Invalid circle */
1704  p = lwpoly_construct_circle(srid, x, y, 0, segments_per_quad, LW_TRUE);
1705  CU_ASSERT_TRUE(p != NULL);
1706  lwpoly_free(p);
1707 }
1708 
1709 static void test_kmeans(void)
1710 {
1711  static int cluster_size = 100;
1712  static int num_clusters = 4;
1713  static double spread = 1.5;
1714  int N = cluster_size * num_clusters;
1715  LWGEOM **geoms;
1716  int i, j, k=0;
1717  int *r;
1718 
1719  geoms = lwalloc(sizeof(LWGEOM*) * N);
1720 
1721  for (j = 0; j < num_clusters; j++) {
1722  for (i = 0; i < cluster_size; i++)
1723  {
1724  double u1 = 1.0 * rand() / RAND_MAX;
1725  double u2 = 1.0 * rand() / RAND_MAX;
1726  double z1 = spread * j + sqrt(-2*log2(u1))*cos(2*M_PI*u2);
1727  double z2 = spread * j + sqrt(-2*log2(u1))*sin(2*M_PI*u2);
1728 
1729  LWPOINT *lwp = lwpoint_make2d(SRID_UNKNOWN, z1, z2);
1730  geoms[k++] = lwpoint_as_lwgeom(lwp);
1731  }
1732  }
1733 
1734  r = lwgeom_cluster_kmeans((const LWGEOM **)geoms, N, num_clusters);
1735 
1736  // for (i = 0; i < k; i++)
1737  // {
1738  // printf("[%d] %s\n", r[i], lwgeom_to_ewkt(geoms[i]));
1739  // }
1740 
1741  /* Clean up */
1742  lwfree(r);
1743  for (i = 0; i < k; i++)
1744  lwgeom_free(geoms[i]);
1745  lwfree(geoms);
1746 
1747  return;
1748 }
1749 
1750 static void test_trim_bits(void)
1751 {
1752  POINTARRAY *pta = ptarray_construct_empty(LW_TRUE, LW_TRUE, 2);
1753  POINT4D pt;
1754  LWLINE *line;
1755  int precision;
1756  uint32_t i;
1757 
1758  pt.x = 1.2345678901234;
1759  pt.y = 2.3456789012345;
1760  pt.z = 3.4567890123456;
1761  pt.m = 4.5678901234567;
1762 
1763  ptarray_insert_point(pta, &pt, 0);
1764 
1765  pt.x *= 3;
1766  pt.y *= 3;
1767  pt.y *= 3;
1768  pt.z *= 3;
1769 
1770  ptarray_insert_point(pta, &pt, 1);
1771 
1772  line = lwline_construct(0, NULL, pta);
1773 
1774  for (precision = -15; precision <= 15; precision++)
1775  {
1776  LWLINE *line2 = (LWLINE*) lwgeom_clone_deep((LWGEOM*) line);
1777  lwgeom_trim_bits_in_place((LWGEOM*) line2, precision, precision, precision, precision);
1778 
1779  for (i = 0; i < line->points->npoints; i++)
1780  {
1781  POINT4D pt1 = getPoint4d(line->points, i);
1782  POINT4D pt2 = getPoint4d(line2->points, i);
1783 
1784  CU_ASSERT_DOUBLE_EQUAL(pt1.x, pt2.x, pow(10, -1*precision));
1785  CU_ASSERT_DOUBLE_EQUAL(pt1.y, pt2.y, pow(10, -1*precision));
1786  CU_ASSERT_DOUBLE_EQUAL(pt1.z, pt2.z, pow(10, -1*precision));
1787  CU_ASSERT_DOUBLE_EQUAL(pt1.m, pt2.m, pow(10, -1*precision));
1788  }
1789 
1790  lwline_free(line2);
1791  }
1792 
1793  lwline_free(line);
1794 }
1795 
1796 /*
1797 ** Used by test harness to register the tests in this file.
1798 */
1799 void algorithms_suite_setup(void);
1800 void algorithms_suite_setup(void)
1801 {
1802  CU_pSuite suite = CU_add_suite("algorithm", init_cg_suite, clean_cg_suite);
1803  PG_ADD_TEST(suite,test_lw_segment_side);
1804  PG_ADD_TEST(suite,test_lw_segment_intersects);
1805  PG_ADD_TEST(suite,test_lwline_crossing_short_lines);
1806  PG_ADD_TEST(suite,test_lwline_crossing_long_lines);
1807  PG_ADD_TEST(suite,test_lwline_crossing_bugs);
1808  PG_ADD_TEST(suite,test_lwpoint_set_ordinate);
1809  PG_ADD_TEST(suite,test_lwpoint_get_ordinate);
1810  PG_ADD_TEST(suite,test_point_interpolate);
1811  PG_ADD_TEST(suite,test_lwline_interpolate_points);
1812  PG_ADD_TEST(suite, test_lwline_interpolate_point_3d);
1813  PG_ADD_TEST(suite,test_lwline_clip);
1814  PG_ADD_TEST(suite, test_lwpoly_clip);
1815  PG_ADD_TEST(suite, test_lwtriangle_clip);
1816  PG_ADD_TEST(suite,test_lwline_clip_big);
1817  PG_ADD_TEST(suite,test_lwmline_clip);
1818  PG_ADD_TEST(suite,test_geohash_point);
1819  PG_ADD_TEST(suite,test_geohash_precision);
1820  PG_ADD_TEST(suite,test_geohash);
1821  PG_ADD_TEST(suite,test_geohash_point_as_int);
1822  PG_ADD_TEST(suite, test_geohash_bbox);
1823  PG_ADD_TEST(suite,test_isclosed);
1824  PG_ADD_TEST(suite,test_lwgeom_simplify);
1825  PG_ADD_TEST(suite,test_lw_arc_center);
1826  PG_ADD_TEST(suite,test_point_density);
1827  PG_ADD_TEST(suite,test_kmeans);
1828  PG_ADD_TEST(suite,test_median_handles_3d_correctly);
1829  PG_ADD_TEST(suite,test_median_robustness);
1830  PG_ADD_TEST(suite,test_lwpoly_construct_circle);
1831  PG_ADD_TEST(suite,test_trim_bits);
1832  PG_ADD_TEST(suite,test_lwgeom_remove_repeated_points);
1833 }
test_kmeans
static void test_kmeans(void)
Definition: cu_algorithm.c:1709
test_geohash
static void test_geohash(void)
Definition: cu_algorithm.c:978
test_geohash_bbox
static void test_geohash_bbox(void)
Definition: cu_algorithm.c:1147
test_lw_segment_intersects
static void test_lw_segment_intersects(void)
Definition: cu_algorithm.c:120
pa22
POINTARRAY * pa22
Definition: cu_algorithm.c:27
test_lwline_clip_big
static void test_lwline_clip_big(void)
Definition: cu_algorithm.c:890
test_point_interpolate
static void test_point_interpolate(void)
Definition: cu_algorithm.c:436
l21
LWLINE * l21
Definition: cu_algorithm.c:28
init_cg_suite
static int init_cg_suite(void)
Definition: cu_algorithm.c:38
test_lwline_crossing_bugs
static void test_lwline_crossing_bugs(void)
Definition: cu_algorithm.c:385
test_geohash_point
static void test_geohash_point(void)
Definition: cu_algorithm.c:957
setpoint
#define setpoint(p, x1, y1)
test_isclosed
static void test_isclosed(void)
Definition: cu_algorithm.c:1021
test_lw_arc_center
static void test_lw_arc_center(void)
Definition: cu_algorithm.c:94
test_lwgeom_remove_repeated_points
static void test_lwgeom_remove_repeated_points(void)
Definition: cu_algorithm.c:1163
test_lwpoly_clip
static void test_lwpoly_clip(void)
Definition: cu_algorithm.c:708
test_point_density
static void test_point_density(void)
Definition: cu_algorithm.c:1573
test_lwpoly_construct_circle
static void test_lwpoly_construct_circle(void)
Definition: cu_algorithm.c:1669
test_lwline_interpolate_point_3d
static void test_lwline_interpolate_point_3d(void)
Definition: cu_algorithm.c:547
test_median_robustness
static void test_median_robustness(void)
Definition: cu_algorithm.c:1475
test_weighted_distance
static double test_weighted_distance(const POINT4D *curr, const POINT4D *points, uint32_t npoints)
Definition: cu_algorithm.c:1375
test_lwgeom_simplify
static void test_lwgeom_simplify(void)
Definition: cu_algorithm.c:1255
parse_result
LWGEOM_PARSER_RESULT parse_result
Definition: cu_algorithm.c:31
clean_cg_suite
static int clean_cg_suite(void)
Definition: cu_algorithm.c:52
test_geohash_precision
static void test_geohash_precision(void)
Definition: cu_algorithm.c:923
test_lwline_crossing_short_lines
static void test_lwline_crossing_short_lines(void)
Definition: cu_algorithm.c:262
test_lwmline_clip
static void test_lwmline_clip(void)
Definition: cu_algorithm.c:802
test_trim_bits
static void test_trim_bits(void)
Definition: cu_algorithm.c:1750
test_geohash_point_as_int
static void test_geohash_point_as_int(void)
Definition: cu_algorithm.c:1126
pa21
POINTARRAY * pa21
Definition: cu_algorithm.c:26
test_median_handles_3d_correctly
static void test_median_handles_3d_correctly(void)
Definition: cu_algorithm.c:1366
l22
LWLINE * l22
Definition: cu_algorithm.c:29
test_lwline_crossing_long_lines
static void test_lwline_crossing_long_lines(void)
Definition: cu_algorithm.c:325
do_median_dims_check
static void do_median_dims_check(char *wkt, int expected_dims)
Definition: cu_algorithm.c:1355
algorithms_suite_setup
void algorithms_suite_setup(void)
Definition: cu_algorithm.c:1800
test_lwtriangle_clip
static void test_lwtriangle_clip(void)
Definition: cu_algorithm.c:744
test_lwpoint_get_ordinate
static void test_lwpoint_get_ordinate(void)
Definition: cu_algorithm.c:419
test_lwline_interpolate_points
static void test_lwline_interpolate_points(void)
Definition: cu_algorithm.c:468
do_median_test
static void do_median_test(char *input, char *expected, int fail_if_not_converged, int iter_count)
Definition: cu_algorithm.c:1387
test_lwpoint_set_ordinate
static void test_lwpoint_set_ordinate(void)
Definition: cu_algorithm.c:399
test_lw_segment_side
static void test_lw_segment_side(void)
Definition: cu_algorithm.c:62
test_lwline_clip
static void test_lwline_clip(void)
Definition: cu_algorithm.c:599
precision
static uint8_t precision
Definition: cu_in_twkb.c:25
r
char * r
Definition: cu_in_wkt.c:24
result
char result[OUT_DOUBLE_BUFFER_SIZE]
Definition: cu_print.c:267
gbox_init
void gbox_init(GBOX *gbox)
Zero out all the entries in the GBOX.
Definition: gbox.c:40
cu_error_msg_reset
void cu_error_msg_reset()
Definition: liblwgeom/cunit/cu_tester.c:306
cu_error_msg
char cu_error_msg[MAX_CUNIT_ERROR_LENGTH+1]
Definition: liblwgeom/cunit/cu_tester.c:19
ASSERT_VARLENA_EQUAL
#define ASSERT_VARLENA_EQUAL(v, s)
Definition: liblwgeom/cunit/cu_tester.h:64
ASSERT_POINT4D_EQUAL
#define ASSERT_POINT4D_EQUAL(o, e, eps)
Definition: liblwgeom/cunit/cu_tester.h:112
ASSERT_INT_EQUAL
#define ASSERT_INT_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:38
PG_ADD_TEST
#define PG_ADD_TEST(suite, testfunc)
Definition: liblwgeom/cunit/cu_tester.h:21
ASSERT_STRING_EQUAL
#define ASSERT_STRING_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:62
ASSERT_POINT2D_EQUAL
#define ASSERT_POINT2D_EQUAL(o, e, eps)
Definition: liblwgeom/cunit/cu_tester.h:107
lwpoly_construct_circle
LWPOLY * lwpoly_construct_circle(int32_t srid, double x, double y, double radius, uint32_t segments_per_quarter, char exterior)
Definition: lwpoly.c:120
lwpoint_getPoint4d_p
int lwpoint_getPoint4d_p(const LWPOINT *point, POINT4D *out)
Definition: lwpoint.c:57
lwgeom_as_lwline
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:162
lwline_interpolate_point_3d
LWPOINT * lwline_interpolate_point_3d(const LWLINE *line, double distance)
Interpolate one point along a line in 3D.
Definition: lwline.c:605
getPoint4d
POINT4D getPoint4d(const POINTARRAY *pa, uint32_t n)
Definition: lwgeom_api.c:108
lwline_as_lwgeom
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:322
lwgeom_simplify
LWGEOM * lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
Simplification.
Definition: lwgeom.c:1847
LW_FALSE
#define LW_FALSE
Definition: liblwgeom.h:108
lwgeom_ndims
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:938
lwgeom_get_srid
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:910
lwpoint_make2d
LWPOINT * lwpoint_make2d(int32_t srid, double x, double y)
Definition: lwpoint.c:163
lwmpoint_free
void lwmpoint_free(LWMPOINT *mpt)
Definition: lwmpoint.c:72
lwgeom_as_lwmpoint
LWMPOINT * lwgeom_as_lwmpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:225
lwgeom_trim_bits_in_place
void lwgeom_trim_bits_in_place(LWGEOM *geom, int32_t prec_x, int32_t prec_y, int32_t prec_z, int32_t prec_m)
Trim the bits of an LWGEOM in place, to optimize it for compression.
Definition: lwgeom.c:2519
lwpoint_free
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:213
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1138
LW_PARSER_CHECK_NONE
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:2085
lwgeom_geohash
lwvarlena_t * lwgeom_geohash(const LWGEOM *lwgeom, int precision)
Calculate the GeoHash (http://geohash.org) string for a geometry.
Definition: lwalgorithm.c:863
LW_SUCCESS
#define LW_SUCCESS
Definition: liblwgeom.h:111
geohash_point_as_int
unsigned int geohash_point_as_int(POINT2D *pt)
Definition: lwalgorithm.c:664
lwpoly_as_lwgeom
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:312
lwline_crossing_direction
int lwline_crossing_direction(const LWLINE *l1, const LWLINE *l2)
Given two lines, characterize how (and if) they cross each other.
Definition: lwalgorithm.c:462
lwpoint_get_x
double lwpoint_get_x(const LWPOINT *point)
Definition: lwpoint.c:63
lwgeom_clone_deep
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:512
lwgeom_remove_repeated_points_in_place
int lwgeom_remove_repeated_points_in_place(LWGEOM *in, double tolerance)
Definition: lwgeom.c:1554
lwgeom_cluster_kmeans
int * lwgeom_cluster_kmeans(const LWGEOM **geoms, uint32_t ngeoms, uint32_t k)
Take a list of LWGEOMs and a number of clusters and return an integer array indicating which cluster ...
Definition: lwkmeans.c:206
ptarray_construct
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:51
lwgeom_has_z
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:917
POINTTYPE
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:116
lwgeom_to_ewkt
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:548
lwline_construct
LWLINE * lwline_construct(int32_t srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
lwgeom_area
double lwgeom_area(const LWGEOM *geom)
Definition: lwgeom.c:1862
ptarray_insert_point
int ptarray_insert_point(POINTARRAY *pa, const POINT4D *p, uint32_t where)
Insert a point into an existing POINTARRAY.
Definition: ptarray.c:85
lwgeom_count_vertices
uint32_t lwgeom_count_vertices(const LWGEOM *geom)
Count the total number of vertices in any LWGEOM.
Definition: lwgeom.c:1229
distance3d_pt_pt
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:1032
lwline_interpolate_points
POINTARRAY * lwline_interpolate_points(const LWLINE *line, double length_fraction, char repeat)
Interpolate one or more points along a line.
Definition: lwline.c:528
lwfree
void lwfree(void *mem)
Definition: lwutil.c:242
lwpoint_as_lwgeom
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:327
lwgeom_clip_to_ordinate_range
LWCOLLECTION * lwgeom_clip_to_ordinate_range(const LWGEOM *lwin, char ordinate, double from, double to, double offset)
Given a geometry clip based on the from/to range of one of its ordinates (x, y, z,...
Definition: lwlinearreferencing.c:829
LINE_MULTICROSS_END_RIGHT
@ LINE_MULTICROSS_END_RIGHT
Definition: liblwgeom.h:1653
LINE_MULTICROSS_END_SAME_FIRST_LEFT
@ LINE_MULTICROSS_END_SAME_FIRST_LEFT
Definition: liblwgeom.h:1654
LINE_MULTICROSS_END_LEFT
@ LINE_MULTICROSS_END_LEFT
Definition: liblwgeom.h:1652
LINE_CROSS_LEFT
@ LINE_CROSS_LEFT
Definition: liblwgeom.h:1650
LINE_CROSS_RIGHT
@ LINE_CROSS_RIGHT
Definition: liblwgeom.h:1651
LINE_NO_CROSS
@ LINE_NO_CROSS
Definition: liblwgeom.h:1649
ptarray_construct_empty
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:59
lwcollection_free
void lwcollection_free(LWCOLLECTION *col)
Definition: lwcollection.c:357
ptarray_free
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:327
lwgeom_get_bbox
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:726
lwgeom_from_wkt
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition: lwin_wkt.c:905
lwalloc
void * lwalloc(size_t size)
Definition: lwutil.c:227
lwpoly_free
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:175
lwgeom_median
LWPOINT * lwgeom_median(const LWGEOM *g, double tol, uint32_t maxiter, char fail_if_not_converged)
Definition: lwgeom_median.c:280
lwgeom_to_points
LWMPOINT * lwgeom_to_points(const LWGEOM *lwgeom, uint32_t npoints, int32_t seed)
Definition: liblwgeom/lwgeom_geos.c:1850
lwmline_free
void lwmline_free(LWMLINE *mline)
Definition: lwmline.c:112
LW_TRUE
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:107
SRID_UNKNOWN
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:229
lwgeom_has_m
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:924
ptarray_set_point4d
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:370
lwline_free
void lwline_free(LWLINE *line)
Definition: lwline.c:67
lwline_construct_empty
LWLINE * lwline_construct_empty(int32_t srid, char hasz, char hasm)
Definition: lwline.c:55
lwpoint_get_y
double lwpoint_get_y(const LWPOINT *point)
Definition: lwpoint.c:76
lwgeom_grid_in_place
void lwgeom_grid_in_place(LWGEOM *lwgeom, const gridspec *grid)
Definition: lwgeom.c:2143
lwpoint_set_ordinate
void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
Given a point, ordinate number and value, set that ordinate on the point.
Definition: lwlinearreferencing.c:283
lwgeom_geohash_precision
int lwgeom_geohash_precision(GBOX bbox, GBOX *bounds)
Definition: lwalgorithm.c:766
lwcircstring_is_closed
int lwcircstring_is_closed(const LWCIRCSTRING *curve)
Definition: lwcircstring.c:261
lwpoint_get_ordinate
double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
Given a POINT4D and an ordinate number, return the value of the ordinate.
Definition: lwlinearreferencing.c:255
geohash_point
lwvarlena_t * geohash_point(double longitude, double latitude, int precision)
Definition: lwalgorithm.c:598
lwline_is_closed
int lwline_is_closed(const LWLINE *line)
Definition: lwline.c:445
lwcompound_is_closed
int lwcompound_is_closed(const LWCOMPOUND *curve)
Definition: lwcompound.c:35
SEG_NO_INTERSECTION
@ SEG_NO_INTERSECTION
Definition: liblwgeom_internal.h:213
SEG_COLINEAR
@ SEG_COLINEAR
Definition: liblwgeom_internal.h:214
SEG_CROSS_RIGHT
@ SEG_CROSS_RIGHT
Definition: liblwgeom_internal.h:216
SEG_CROSS_LEFT
@ SEG_CROSS_LEFT
Definition: liblwgeom_internal.h:215
lw_segment_intersects
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:373
lwmpoint_extract_points_4d
POINT4D * lwmpoint_extract_points_4d(const LWMPOINT *g, uint32_t *npoints, int *input_empty)
Definition: lwgeom_median.c:172
lw_arc_center
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:229
lwpoint_is_empty
int lwpoint_is_empty(const LWPOINT *point)
FP_TOLERANCE
#define FP_TOLERANCE
Floating point comparators.
Definition: liblwgeom_internal.h:54
ptarray_has_z
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:37
lw_segment_side
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:65
point_interpolate
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...
Definition: lwlinearreferencing.c:316
decode_geohash_bbox
void decode_geohash_bbox(char *geohash, double *lat, double *lon, int precision)
Definition: lwalgorithm.c:721
ptarray_has_m
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:44
liblwgeom_internal.h
lwgeom_is_empty
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwinline.h:203
lwgeom_as_lwpoint
static LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwinline.h:131
distance
static double distance(double x1, double y1, double x2, double y2)
Definition: lwtree.c:1032
cu_tester.h
GBOX::ymax
double ymax
Definition: liblwgeom.h:371
GBOX::xmax
double xmax
Definition: liblwgeom.h:369
GBOX::ymin
double ymin
Definition: liblwgeom.h:370
GBOX::xmin
double xmin
Definition: liblwgeom.h:368
LWGEOM::type
uint8_t type
Definition: liblwgeom.h:476
LWLINE::points
POINTARRAY * points
Definition: liblwgeom.h:497
LWMPOINT::ngeoms
uint32_t ngeoms
Definition: liblwgeom.h:552
LWMPOINT::geoms
LWPOINT ** geoms
Definition: liblwgeom.h:547
POINT2D::y
double y
Definition: liblwgeom.h:404
POINT2D::x
double x
Definition: liblwgeom.h:404
POINT4D::m
double m
Definition: liblwgeom.h:428
POINT4D::x
double x
Definition: liblwgeom.h:428
POINT4D::z
double z
Definition: liblwgeom.h:428
POINT4D::y
double y
Definition: liblwgeom.h:428
POINTARRAY::npoints
uint32_t npoints
Definition: liblwgeom.h:441
gridspec_t::zsize
double zsize
Definition: liblwgeom.h:1373
gridspec_t::ysize
double ysize
Definition: liblwgeom.h:1372
gridspec_t::xsize
double xsize
Definition: liblwgeom.h:1371
gridspec_t::msize
double msize
Definition: liblwgeom.h:1374
gridspec_t
Snap-to-grid.
Definition: liblwgeom.h:1366
lwvarlena_t
Definition: liblwgeom.h:320
struct_lwgeom_parser_result
Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM.
Definition: liblwgeom.h:2093