PostGIS  3.3.9dev-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, *gn;
1166  char *ewkt;
1167  char *ewkt_exp;
1168  int modified = LW_FALSE;
1169 
1170  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);
1171  modified = lwgeom_remove_repeated_points_in_place(g, 1);
1172  ASSERT_INT_EQUAL(modified, LW_TRUE);
1173  gn = lwgeom_normalize(g);
1174  ewkt = lwgeom_to_ewkt(gn);
1175  ASSERT_STRING_EQUAL(ewkt, "MULTIPOINT(10 10,10 0,5 5,0 10,0 0)");
1176  lwgeom_free(g);
1177  lwgeom_free(gn);
1178  lwfree(ewkt);
1179 
1180  g = lwgeom_from_wkt("LINESTRING(1612830.15445 4841287.12672,1612830.15824 4841287.12674,1612829.98813 4841274.56198)", LW_PARSER_CHECK_NONE);
1181  modified = lwgeom_remove_repeated_points_in_place(g, 0.01);
1182  ASSERT_INT_EQUAL(modified, LW_TRUE);
1183  ewkt = lwgeom_to_ewkt(g);
1184  ASSERT_STRING_EQUAL(ewkt, "LINESTRING(1612830.15445 4841287.12672,1612829.98813 4841274.56198)");
1185  lwgeom_free(g);
1186  lwfree(ewkt);
1187 
1188  /* remove points */
1189  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);
1190  modified = lwgeom_remove_repeated_points_in_place(g, 1);
1191  ASSERT_INT_EQUAL(modified, LW_TRUE);
1192  gn = lwgeom_normalize(g);
1193  ewkt = lwgeom_to_ewkt(gn);
1194  lwgeom_free(g);
1195  lwgeom_free(gn);
1196  /* build expected and normalize */
1197  g = lwgeom_from_wkt("MULTIPOINT(0 0,0 10,5 5,5 8,8 5,8 8,10 0,10 10,50 50,50 60,55 55,55 58,58 55,58 58,60 50,60 60)", LW_PARSER_CHECK_NONE);
1198  gn = lwgeom_normalize(g);
1199  ewkt_exp = lwgeom_to_ewkt(gn);
1200  ASSERT_STRING_EQUAL(ewkt, ewkt_exp);
1201  lwgeom_free(g);
1202  lwgeom_free(gn);
1203  lwfree(ewkt);
1204  lwfree(ewkt_exp);
1205 
1206  g = lwgeom_from_wkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", LW_PARSER_CHECK_NONE);
1207  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1208  ASSERT_INT_EQUAL(modified, LW_TRUE);
1209  ewkt = lwgeom_to_ewkt(g);
1210  ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 1,1 0,0 0))");
1211  lwgeom_free(g);
1212  lwfree(ewkt);
1213 
1214  // Test the return value (modified or not)
1215  g = lwgeom_from_wkt("POINT(0 0)", LW_PARSER_CHECK_NONE);
1216  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1217  ASSERT_INT_EQUAL(modified, LW_FALSE);
1218  lwgeom_free(g);
1219 
1220  g = lwgeom_from_wkt("TRIANGLE((0 0, 5 0, 3 3, 0 0))", LW_PARSER_CHECK_NONE);
1221  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1222  ASSERT_INT_EQUAL(modified, LW_FALSE);
1223  lwgeom_free(g);
1224 
1225  g = lwgeom_from_wkt("POLYGON((0 0,0 1,1 1,1 0,0 0))", LW_PARSER_CHECK_NONE);
1226  modified = lwgeom_remove_repeated_points_in_place(g, 0.1);
1227  ASSERT_INT_EQUAL(modified, LW_FALSE);
1228  ewkt = lwgeom_to_ewkt(g);
1229  ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,0 1,1 1,1 0,0 0))");
1230  lwgeom_free(g);
1231  lwfree(ewkt);
1232 
1233  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))",
1234  LW_PARSER_CHECK_NONE);
1235  modified = lwgeom_remove_repeated_points_in_place(g, 0.1);
1236  ASSERT_INT_EQUAL(modified, LW_TRUE);
1237  ewkt = lwgeom_to_ewkt(g);
1238  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))");
1239  lwgeom_free(g);
1240  lwfree(ewkt);
1241 
1242  g = lwgeom_from_wkt("GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0,1 0,1 1,0 1,0 0)))", LW_PARSER_CHECK_NONE);
1243  modified = lwgeom_remove_repeated_points_in_place(g, 0.1);
1244  ASSERT_INT_EQUAL(modified, LW_FALSE);
1245  ewkt = lwgeom_to_ewkt(g);
1246  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0,1 0,1 1,0 1,0 0)))");
1247  lwgeom_free(g);
1248  lwfree(ewkt);
1249 
1250  g = lwgeom_from_wkt("GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0, 0 1, 1 1, 1 0, 0 0)))", LW_PARSER_CHECK_NONE);
1251  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1252  ASSERT_INT_EQUAL(modified, LW_TRUE);
1253  ewkt = lwgeom_to_ewkt(g);
1254  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0,1 1,1 0,0 0)))");
1255  lwgeom_free(g);
1256  lwfree(ewkt);
1257 
1258  g = lwgeom_from_wkt("GEOMETRYCOLLECTION(POLYGON((0 0, 0 1, 1 1, 1 0, 0 0)),POINT(2 0))", LW_PARSER_CHECK_NONE);
1259  modified = lwgeom_remove_repeated_points_in_place(g, 10000);
1260  ASSERT_INT_EQUAL(modified, LW_TRUE);
1261  ewkt = lwgeom_to_ewkt(g);
1262  ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POLYGON((0 0,1 1,1 0,0 0)),POINT(2 0))");
1263  lwgeom_free(g);
1264  lwfree(ewkt);
1265 }
1266 
1267 static void test_lwgeom_simplify(void)
1268 {
1269  LWGEOM *l;
1270  LWGEOM *g;
1271  char *ewkt;
1272 
1273  /* Simplify but only so far... */
1274  g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
1275  l = lwgeom_simplify(g, 10, LW_TRUE);
1276  ewkt = lwgeom_to_ewkt(l);
1277  ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)");
1278  lwgeom_free(g);
1279  lwgeom_free(l);
1280  lwfree(ewkt);
1281 
1282  /* Simplify but only so far... */
1283  g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
1284  l = lwgeom_simplify(g, 10, LW_TRUE);
1285  ewkt = lwgeom_to_ewkt(l);
1286  ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))");
1287  lwgeom_free(g);
1288  lwgeom_free(l);
1289  lwfree(ewkt);
1290 
1291  /* Simplify and collapse */
1292  g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
1293  l = lwgeom_simplify(g, 10, LW_FALSE);
1294  CU_ASSERT_EQUAL(l, NULL);
1295  lwgeom_free(g);
1296  lwgeom_free(l);
1297 
1298  /* Simplify and collapse */
1299  g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
1300  l = lwgeom_simplify(g, 10, LW_FALSE);
1301  CU_ASSERT_EQUAL(l, NULL);
1302  lwgeom_free(g);
1303  lwgeom_free(l);
1304 
1305  /* Not simplifiable */
1306  g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE);
1307  l = lwgeom_simplify(g, 1.0, LW_FALSE);
1308  ewkt = lwgeom_to_ewkt(l);
1309  ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)");
1310  lwgeom_free(g);
1311  lwgeom_free(l);
1312  lwfree(ewkt);
1313 
1314  /* Simplifiable */
1315  g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE);
1316  l = lwgeom_simplify(g, 1.0, LW_FALSE);
1317  ewkt = lwgeom_to_ewkt(l);
1318  ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)");
1319  lwgeom_free(g);
1320  lwgeom_free(l);
1321  lwfree(ewkt);
1322 
1323  /* POLYGON with multiple inner rings*/
1324  g = lwgeom_from_wkt(
1325  "POLYGON("
1326  "(0 0, 100 0, 100 100, 0 100, 0 0),"
1327  "(1 1, 1 5, 5 5, 5 1, 1 1),"
1328  "(20 20, 20 40, 40 40, 40 20, 20 20)"
1329  ")",
1330  LW_PARSER_CHECK_NONE);
1331  l = lwgeom_simplify(g, 10, LW_FALSE);
1332  ewkt = lwgeom_to_ewkt(l);
1333  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))");
1334  lwgeom_free(g);
1335  lwgeom_free(l);
1336  lwfree(ewkt);
1337 
1338  /* Reorder inner rings: Same result */
1339  g = lwgeom_from_wkt(
1340  "POLYGON("
1341  "(0 0, 100 0, 100 100, 0 100, 0 0),"
1342  "(20 20, 20 40, 40 40, 40 20, 20 20),"
1343  "(1 1, 1 5, 5 5, 5 1, 1 1)"
1344  ")",
1345  LW_PARSER_CHECK_NONE);
1346  l = lwgeom_simplify(g, 10, LW_FALSE);
1347  ewkt = lwgeom_to_ewkt(l);
1348  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))");
1349  lwgeom_free(g);
1350  lwgeom_free(l);
1351  lwfree(ewkt);
1352 
1353  g = lwgeom_from_wkt(
1354  "POLYGON("
1355  "(0 0, 100 0, 100 100, 0 100, 0 0),"
1356  "(20 20, 20 40, 40 40, 40 20, 20 20),"
1357  "(1 1, 1 5, 5 5, 5 1, 1 1)"
1358  ")",
1359  LW_PARSER_CHECK_NONE);
1360  l = lwgeom_simplify(g, 100, LW_FALSE);
1361  CU_ASSERT_EQUAL(l, NULL);
1362  lwgeom_free(g);
1363  lwgeom_free(l);
1364 }
1365 
1366 
1367 static void do_median_dims_check(char* wkt, int expected_dims)
1368 {
1369  LWGEOM* g = lwgeom_from_wkt(wkt, LW_PARSER_CHECK_NONE);
1370  LWPOINT* result = lwgeom_median(g, 1e-8, 100, LW_FALSE);
1371 
1372  CU_ASSERT_EQUAL(expected_dims, lwgeom_ndims((LWGEOM*) result));
1373 
1374  lwgeom_free(g);
1375  lwpoint_free(result);
1376 }
1377 
1378 static void test_median_handles_3d_correctly(void)
1379 {
1380  do_median_dims_check("MULTIPOINT ((1 3), (4 7), (2 9), (0 4), (2 2))", 2);
1381  do_median_dims_check("MULTIPOINT Z ((1 3 4), (4 7 8), (2 9 1), (0 4 4), (2 2 3))", 3);
1382  do_median_dims_check("MULTIPOINT M ((1 3 4), (4 7 8), (2 9 1), (0 4 4), (2 2 3))", 2);
1383  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);
1384 }
1385 
1386 static double
1387 test_weighted_distance(const POINT4D* curr, const POINT4D* points, uint32_t npoints)
1388 {
1389  double distance = 0.0;
1390  uint32_t i;
1391  for (i = 0; i < npoints; i++)
1392  {
1393  distance += distance3d_pt_pt((POINT3D*)curr, (POINT3D*)&points[i]) * points[i].m;
1394  }
1395 
1396  return distance;
1397 }
1398 
1399 static void do_median_test(char* input, char* expected, int fail_if_not_converged, int iter_count)
1400 {
1401  cu_error_msg_reset();
1402  LWGEOM* g = lwgeom_from_wkt(input, LW_PARSER_CHECK_NONE);
1403  LWPOINT* expected_result = NULL;
1404  POINT4D actual_pt;
1405  POINT4D expected_pt;
1406  const double tolerance = FP_TOLERANCE / 10.0;
1407 
1408  LWPOINT* result = lwgeom_median(g, tolerance, iter_count, fail_if_not_converged);
1409  int passed = LW_FALSE;
1410 
1411  if (expected != NULL)
1412  {
1413  expected_result = lwgeom_as_lwpoint(lwgeom_from_wkt(expected, LW_PARSER_CHECK_NONE));
1414  lwpoint_getPoint4d_p(expected_result, &expected_pt);
1415  }
1416  if (result != NULL)
1417  {
1418  lwpoint_getPoint4d_p(result, &actual_pt);
1419  }
1420 
1421  if (result != NULL && expected != NULL) /* got something, expecting something */
1422  {
1423  passed = LW_TRUE;
1424  passed = passed && lwgeom_is_empty((LWGEOM*) expected_result) == lwgeom_is_empty((LWGEOM*) result);
1425  passed = passed && (lwgeom_has_z((LWGEOM*) expected_result) == lwgeom_has_z((LWGEOM*) result));
1426 
1427  if (passed && !lwgeom_is_empty((LWGEOM*) result))
1428  {
1429  if (g->type == POINTTYPE)
1430  {
1431  passed &= fabs(actual_pt.x - expected_pt.x) < tolerance;
1432  passed &= fabs(actual_pt.y - expected_pt.y) < tolerance;
1433  passed &= (!lwgeom_has_z((LWGEOM*) expected_result) || fabs(actual_pt.z - expected_pt.z) < tolerance);
1434  passed &= (!lwgeom_has_m((LWGEOM*) expected_result) || fabs(actual_pt.m - expected_pt.m) < tolerance);
1435  }
1436  else
1437  {
1438  /* Check that the difference between the obtained geometric
1439  median and the expected point is within tolerance */
1440  uint32_t npoints = 1;
1441  int input_empty = LW_TRUE;
1442  POINT4D* points = lwmpoint_extract_points_4d(lwgeom_as_lwmpoint(g), &npoints, &input_empty);
1443  double distance_expected = test_weighted_distance(&expected_pt, points, npoints);
1444  double distance_result = test_weighted_distance(&actual_pt, points, npoints);
1445 
1446  passed = distance_result <= (1.0 + tolerance) * distance_expected;
1447  if (!passed)
1448  {
1449  printf("Diff: Got %.10f Expected %.10f\n", distance_result, distance_expected);
1450  }
1451  lwfree(points);
1452  }
1453  }
1454 
1455  if (!passed)
1456  {
1457  printf("median_test input %s (parsed %s) expected %s got %s\n",
1458  input, lwgeom_to_ewkt(g),
1459  lwgeom_to_ewkt((LWGEOM*) expected_result),
1460  lwgeom_to_ewkt((LWGEOM*) result));
1461  }
1462 
1463  }
1464  else if (result == NULL && expected == NULL) /* got nothing, expecting nothing */
1465  {
1466  passed = LW_TRUE;
1467  }
1468  else if (result != NULL && expected == NULL) /* got something, expecting nothing */
1469  {
1470  passed = LW_FALSE;
1471  printf("median_test input %s (parsed %s) expected NULL got %s\n", input, lwgeom_to_ewkt(g), lwgeom_to_ewkt((LWGEOM*) result));
1472  }
1473  else if (result == NULL && expected != NULL) /* got nothing, expecting something */
1474  {
1475  passed = LW_FALSE;
1476  printf("%s", cu_error_msg);
1477  printf("median_test input %s (parsed %s) expected %s got NULL\n", input, lwgeom_to_ewkt(g), lwgeom_to_ewkt((LWGEOM*) expected_result));
1478  }
1479 
1480  CU_ASSERT_TRUE(passed);
1481 
1482  lwgeom_free(g);
1483  lwpoint_free(expected_result);
1484  lwpoint_free(result);
1485 }
1486 
1487 static void test_median_robustness(void)
1488 {
1489  /* A simple implementation of Weiszfeld's algorithm will fail if the median is equal
1490  * to any one of the inputs, during any iteration of the algorithm.
1491  *
1492  * Because the algorithm uses the centroid as a starting point, this situation will
1493  * occur in the test case below.
1494  */
1495  do_median_test("MULTIPOINT ((0 -1), (0 0), (0 1))", "POINT (0 0)", LW_TRUE, 1000);
1496 
1497  /* Same as above but 3D, and shifter */
1498  do_median_test("MULTIPOINT ((1 -1 3), (1 0 2), (2 1 1))", "POINT (1 0 2)", LW_TRUE, 1000);
1499 
1500  /* Starting point is duplicated */
1501  do_median_test("MULTIPOINT ((0 -1), (0 0), (0 0), (0 1))", "POINT (0 0)", LW_TRUE, 1000);
1502 
1503  /* Cube */
1504  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))",
1505  "POINT (15 15 15)", LW_TRUE, 1000);
1506 
1507  /* Some edge cases */
1508  do_median_test("POINT (7 6)", "POINT (7 6)", LW_TRUE, 1000);
1509  do_median_test("POINT (7 6 2)", "POINT (7 6 2)", LW_TRUE, 1000);
1510  do_median_test("MULTIPOINT ((7 6 2), EMPTY)", "POINT (7 6 2)", LW_TRUE, 1000);
1511 
1512  /* Empty input */
1513  do_median_test("MULTIPOINT EMPTY", "POINT EMPTY", LW_FALSE, 1000);
1514  do_median_test("MULTIPOINT (EMPTY)", "POINT EMPTY", LW_FALSE, 1000);
1515  do_median_test("MULTIPOINT EMPTY", "POINT EMPTY", LW_TRUE, 1000);
1516  do_median_test("MULTIPOINT (EMPTY)", "POINT EMPTY", LW_TRUE, 1000);
1517  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);
1518 
1519  /* Weighted input */
1520  do_median_test("MULTIPOINT ZM (1 -1 3 1, 1 0 2 7, 2 1 1 1)", "POINT (1 0 2)", LW_TRUE, 1000);
1521  do_median_test("MULTIPOINT ZM (-1 1 -3 1, -1 0 -2 7, -2 -1 -1 1)", "POINT (-1 0 -2)", LW_TRUE, 1000);
1522  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);
1523 
1524  /* Point that is replaced by two half-weighted */
1525  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);
1526  /* Point is doubled and then erased by negative weight */
1527  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);
1528  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);
1529  /* Weightless input won't converge */
1530  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);
1531  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);
1532  /* Negative weight won't converge */
1533  do_median_test("MULTIPOINT ZM ((0 -1 0 -1), (0 0 0 -1), (0 1 0 -1))", NULL, LW_FALSE, 1000);
1534  do_median_test("MULTIPOINT ZM ((0 -1 0 -1), (0 0 0 -1), (0 1 0 -1))", NULL, LW_TRUE, 1000);
1535 
1536  /* Bind convergence too tightly */
1537  do_median_test("MULTIPOINT ((0 0), (1 1), (0 1), (2 2))", "POINT(0.75 1.0)", LW_FALSE, 0);
1538  do_median_test("MULTIPOINT ((0 0), (1 1), (0 1), (2 2))", NULL, LW_TRUE, 1);
1539  /* Unsupported geometry type */
1540  do_median_test("POLYGON((1 0,0 1,1 2,2 1,1 0))", NULL, LW_TRUE, 1000);
1541  do_median_test("POLYGON((1 0,0 1,1 2,2 1,1 0))", NULL, LW_FALSE, 1000);
1542 
1543  /* Median point is included */
1544  do_median_test("MULTIPOINT ZM ("
1545  "(1480 0 200 100),"
1546  "(620 0 200 100),"
1547  "(1000 0 -200 100),"
1548  "(1000 0 -590 100),"
1549  "(1025 0 65 100),"
1550  "(1025 0 -65 100)"
1551  ")",
1552  "POINT (1025 0 -65)", LW_TRUE, 10000);
1553 
1554 #if 0
1555  /* Leads to invalid result (0 0 0) with 80bit (fmulp + faddp) precision. ok with 64 bit float ops */
1556  do_median_test("MULTIPOINT ZM ("
1557  "(0 0 20000 0.5),"
1558  "(0 0 59000 0.5),"
1559  "(0 -3000 -3472.22222222222262644208967685699462890625 1),"
1560  "(0 3000 3472.22222222222262644208967685699462890625 1),"
1561  "(0 0 -1644.736842105263121993630193173885345458984375 1),"
1562  "(0 0 1644.736842105263121993630193173885345458984375 1),"
1563  "(0 48000 -20000 1.3),"
1564  "(0 -48000 -20000 1.3)"
1565  ")",
1566  "POINT (0 0 0)", LW_TRUE, 10000);
1567 #endif
1568 
1569 #if 0
1570  /* Leads to invalid result (0 0 0) with 64bit (vfmadd231sd) precision. Ok with 80 bit float ops */
1571  do_median_test("MULTIPOINT ZM ("
1572  "(0 0 20000 0.5),"
1573  "(0 0 59000 0.5),"
1574  "(0 -3000 -3472.22222222222262644208967685699462890625 1),"
1575  "(0 3000 3472.22222222222262644208967685699462890625 1),"
1576  "(0 -0.00000000000028047739569477638384522295466033823196 -1644.736842105263121993630193173885345458984375 1),"
1577  "(0 0.00000000000028047739569477638384522295466033823196 1644.736842105263121993630193173885345458984375 1),"
1578  "(0 48000 -20000 1.3),"
1579  "(0 -48000 -20000 1.3)"
1580  ")",
1581  "POINT (0 0 0)", LW_TRUE, 10000);
1582 #endif
1583 }
1584 
1585 static void test_point_density(void)
1586 {
1587  LWGEOM *geom;
1588  LWMPOINT *mpt;
1589  LWMPOINT *mpt2;
1590  LWPOINT *pt;
1591  LWPOINT *pt2;
1592  int eq, i;
1593  // char *ewkt;
1594 
1595  /* POLYGON */
1596  geom = lwgeom_from_wkt("POLYGON((0 0,1 0,1 1,0 1,0 0))", LW_PARSER_CHECK_NONE);
1597  mpt = lwgeom_to_points(geom, 100, 0); /* Set a zero seed to base it on Unix time and process ID */
1598  CU_ASSERT_EQUAL(mpt->ngeoms,100);
1599 
1600  /* Run a second time with a zero seed to get a different multipoint sequence */
1601  mpt2 = lwgeom_to_points(geom, 100, 0);
1602  eq = 0;
1603  for (i = 0; i < 100; i++)
1604  {
1605  pt = (LWPOINT*)mpt->geoms[i];
1606  pt2 = (LWPOINT*)mpt2->geoms[i];
1607  if (lwpoint_get_x(pt) == lwpoint_get_x(pt2) && lwpoint_get_y(pt) == lwpoint_get_y(pt2))
1608  eq++;
1609  }
1610  CU_ASSERT_EQUAL(eq, 0);
1611  lwmpoint_free(mpt);
1612  lwmpoint_free(mpt2);
1613  pt = NULL;
1614  pt2 = NULL;
1615 
1616  /* Set seed to get a deterministic sequence */
1617  mpt = lwgeom_to_points(geom, 1000, 12345);
1618 
1619  /* Check to find a different multipoint sequence with different seed */
1620  mpt2 = lwgeom_to_points(geom, 1000, 54321);
1621  eq = 0;
1622  for (i = 0; i < 1000; i++)
1623  {
1624  pt = (LWPOINT*)mpt->geoms[i];
1625  pt2 = (LWPOINT*)mpt2->geoms[i];
1626  if (lwpoint_get_x(pt) == lwpoint_get_x(pt2) && lwpoint_get_y(pt) == lwpoint_get_y(pt2))
1627  eq++;
1628  }
1629  CU_ASSERT_EQUAL(eq, 0);
1630  lwmpoint_free(mpt2);
1631  pt = NULL;
1632  pt2 = NULL;
1633 
1634  /* Check to find an identical multipoint sequence with same seed */
1635  mpt2 = lwgeom_to_points(geom, 1000, 12345);
1636  eq = 0;
1637  for (i = 0; i < 1000; i++)
1638  {
1639  pt = (LWPOINT*)mpt->geoms[i];
1640  pt2 = (LWPOINT*)mpt2->geoms[i];
1641  if (lwpoint_get_x(pt) == lwpoint_get_x(pt2) && lwpoint_get_y(pt) == lwpoint_get_y(pt2))
1642  eq++;
1643  }
1644  CU_ASSERT_EQUAL(eq, 1000);
1645  lwmpoint_free(mpt2);
1646  pt = NULL;
1647  pt2 = NULL;
1648 
1649 
1650  /* Check if the 1000th point is the expected value.
1651  * Note that if the RNG is not portable, this test may fail. */
1652  pt = (LWPOINT*)mpt->geoms[999];
1653  // ewkt = lwgeom_to_ewkt((LWGEOM*)pt);
1654  // printf("%s\n", ewkt);
1655  // lwfree(ewkt);
1656  CU_ASSERT_DOUBLE_EQUAL(lwpoint_get_x(pt), 0.801167838758, 1e-11);
1657  CU_ASSERT_DOUBLE_EQUAL(lwpoint_get_y(pt), 0.345281131175, 1e-11);
1658  lwmpoint_free(mpt);
1659  pt = NULL;
1660 
1661  mpt = lwgeom_to_points(geom, 0, 0);
1662  CU_ASSERT_EQUAL(mpt, NULL);
1663  lwmpoint_free(mpt);
1664 
1665  lwgeom_free(geom);
1666 
1667  /* MULTIPOLYGON */
1668  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);
1669 
1670  mpt = lwgeom_to_points(geom, 1000, 0);
1671  CU_ASSERT_EQUAL(mpt->ngeoms,1000);
1672  lwmpoint_free(mpt);
1673 
1674  mpt = lwgeom_to_points(geom, 1, 0);
1675  CU_ASSERT_EQUAL(mpt->ngeoms,1);
1676  lwmpoint_free(mpt);
1677 
1678  lwgeom_free(geom);
1679 }
1680 
1681 static void test_lwpoly_construct_circle(void)
1682 {
1683  LWPOLY* p;
1684  const GBOX* g;
1685  const int32_t srid = 4326;
1686  const uint32_t segments_per_quad = 17;
1687  const int x = 10;
1688  const int y = 20;
1689  const int r = 5;
1690 
1691  /* With normal arguments you should get something circle-y */
1692  p = lwpoly_construct_circle(srid, x, y, r, segments_per_quad, LW_TRUE);
1693 
1694  ASSERT_INT_EQUAL(lwgeom_count_vertices(lwpoly_as_lwgeom(p)), segments_per_quad * 4 + 1);
1695  ASSERT_INT_EQUAL(lwgeom_get_srid(lwpoly_as_lwgeom(p)), srid)
1696 
1697  g = lwgeom_get_bbox(lwpoly_as_lwgeom(p));
1698  CU_ASSERT_DOUBLE_EQUAL(g->xmin, x-r, 0.1);
1699  CU_ASSERT_DOUBLE_EQUAL(g->xmax, x+r, 0.1);
1700  CU_ASSERT_DOUBLE_EQUAL(g->ymin, y-r, 0.1);
1701  CU_ASSERT_DOUBLE_EQUAL(g->ymax, y+r, 0.1);
1702 
1703  CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(lwpoly_as_lwgeom(p)), M_PI*5*5, 0.1);
1704 
1705  lwpoly_free(p);
1706 
1707  /* No segments? No circle. */
1708  p = lwpoly_construct_circle(srid, x, y, r, 0, LW_TRUE);
1709  CU_ASSERT_TRUE(p == NULL);
1710 
1711  /* Negative radius? No circle. */
1712  p = lwpoly_construct_circle(srid, x, y, -1e-3, segments_per_quad, LW_TRUE);
1713  CU_ASSERT_TRUE(p == NULL);
1714 
1715  /* Zero radius? Invalid circle */
1716  p = lwpoly_construct_circle(srid, x, y, 0, segments_per_quad, LW_TRUE);
1717  CU_ASSERT_TRUE(p != NULL);
1718  lwpoly_free(p);
1719 }
1720 
1721 static void test_kmeans(void)
1722 {
1723  static int cluster_size = 100;
1724  static int num_clusters = 4;
1725  static double spread = 1.5;
1726  int N = cluster_size * num_clusters;
1727  LWGEOM **geoms;
1728  int i, j, k=0;
1729  int *r;
1730 
1731  geoms = lwalloc(sizeof(LWGEOM*) * N);
1732 
1733  for (j = 0; j < num_clusters; j++) {
1734  for (i = 0; i < cluster_size; i++)
1735  {
1736  double u1 = 1.0 * rand() / RAND_MAX;
1737  double u2 = 1.0 * rand() / RAND_MAX;
1738  double z1 = spread * j + sqrt(-2*log2(u1))*cos(2*M_PI*u2);
1739  double z2 = spread * j + sqrt(-2*log2(u1))*sin(2*M_PI*u2);
1740 
1741  LWPOINT *lwp = lwpoint_make2d(SRID_UNKNOWN, z1, z2);
1742  geoms[k++] = lwpoint_as_lwgeom(lwp);
1743  }
1744  }
1745 
1746  r = lwgeom_cluster_kmeans((const LWGEOM **)geoms, N, num_clusters, 0.0);
1747 
1748  // for (i = 0; i < k; i++)
1749  // {
1750  // printf("[%d] %s\n", r[i], lwgeom_to_ewkt(geoms[i]));
1751  // }
1752 
1753  /* Clean up */
1754  lwfree(r);
1755  for (i = 0; i < k; i++)
1756  lwgeom_free(geoms[i]);
1757  lwfree(geoms);
1758 
1759  return;
1760 }
1761 
1762 static void test_trim_bits(void)
1763 {
1764  POINTARRAY *pta = ptarray_construct_empty(LW_TRUE, LW_TRUE, 2);
1765  POINT4D pt;
1766  LWLINE *line;
1767  int precision;
1768  uint32_t i;
1769 
1770  pt.x = 1.2345678901234;
1771  pt.y = 2.3456789012345;
1772  pt.z = 3.4567890123456;
1773  pt.m = 4.5678901234567;
1774 
1775  ptarray_insert_point(pta, &pt, 0);
1776 
1777  pt.x *= 3;
1778  pt.y *= 3;
1779  pt.y *= 3;
1780  pt.z *= 3;
1781 
1782  ptarray_insert_point(pta, &pt, 1);
1783 
1784  line = lwline_construct(0, NULL, pta);
1785 
1786  for (precision = -15; precision <= 15; precision++)
1787  {
1788  LWLINE *line2 = (LWLINE*) lwgeom_clone_deep((LWGEOM*) line);
1789  lwgeom_trim_bits_in_place((LWGEOM*) line2, precision, precision, precision, precision);
1790 
1791  for (i = 0; i < line->points->npoints; i++)
1792  {
1793  POINT4D pt1 = getPoint4d(line->points, i);
1794  POINT4D pt2 = getPoint4d(line2->points, i);
1795 
1796  CU_ASSERT_DOUBLE_EQUAL(pt1.x, pt2.x, pow(10, -1*precision));
1797  CU_ASSERT_DOUBLE_EQUAL(pt1.y, pt2.y, pow(10, -1*precision));
1798  CU_ASSERT_DOUBLE_EQUAL(pt1.z, pt2.z, pow(10, -1*precision));
1799  CU_ASSERT_DOUBLE_EQUAL(pt1.m, pt2.m, pow(10, -1*precision));
1800  }
1801 
1802  lwline_free(line2);
1803  }
1804 
1805  lwline_free(line);
1806 }
1807 
1808 /*
1809 ** Used by test harness to register the tests in this file.
1810 */
1811 void algorithms_suite_setup(void);
1812 void algorithms_suite_setup(void)
1813 {
1814  CU_pSuite suite = CU_add_suite("algorithm", init_cg_suite, clean_cg_suite);
1815  PG_ADD_TEST(suite,test_lw_segment_side);
1816  PG_ADD_TEST(suite,test_lw_segment_intersects);
1817  PG_ADD_TEST(suite,test_lwline_crossing_short_lines);
1818  PG_ADD_TEST(suite,test_lwline_crossing_long_lines);
1819  PG_ADD_TEST(suite,test_lwline_crossing_bugs);
1820  PG_ADD_TEST(suite,test_lwpoint_set_ordinate);
1821  PG_ADD_TEST(suite,test_lwpoint_get_ordinate);
1822  PG_ADD_TEST(suite,test_point_interpolate);
1823  PG_ADD_TEST(suite,test_lwline_interpolate_points);
1824  PG_ADD_TEST(suite, test_lwline_interpolate_point_3d);
1825  PG_ADD_TEST(suite,test_lwline_clip);
1826  PG_ADD_TEST(suite, test_lwpoly_clip);
1827  PG_ADD_TEST(suite, test_lwtriangle_clip);
1828  PG_ADD_TEST(suite,test_lwline_clip_big);
1829  PG_ADD_TEST(suite,test_lwmline_clip);
1830  PG_ADD_TEST(suite,test_geohash_point);
1831  PG_ADD_TEST(suite,test_geohash_precision);
1832  PG_ADD_TEST(suite,test_geohash);
1833  PG_ADD_TEST(suite,test_geohash_point_as_int);
1834  PG_ADD_TEST(suite, test_geohash_bbox);
1835  PG_ADD_TEST(suite,test_isclosed);
1836  PG_ADD_TEST(suite,test_lwgeom_simplify);
1837  PG_ADD_TEST(suite,test_lw_arc_center);
1838  PG_ADD_TEST(suite,test_point_density);
1839  PG_ADD_TEST(suite,test_kmeans);
1840  PG_ADD_TEST(suite,test_median_handles_3d_correctly);
1841  PG_ADD_TEST(suite,test_median_robustness);
1842  PG_ADD_TEST(suite,test_lwpoly_construct_circle);
1843  PG_ADD_TEST(suite,test_trim_bits);
1844  PG_ADD_TEST(suite,test_lwgeom_remove_repeated_points);
1845 }
test_kmeans
static void test_kmeans(void)
Definition: cu_algorithm.c:1721
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:1585
test_lwpoly_construct_circle
static void test_lwpoly_construct_circle(void)
Definition: cu_algorithm.c:1681
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:1487
test_weighted_distance
static double test_weighted_distance(const POINT4D *curr, const POINT4D *points, uint32_t npoints)
Definition: cu_algorithm.c:1387
test_lwgeom_simplify
static void test_lwgeom_simplify(void)
Definition: cu_algorithm.c:1267
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:1762
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:1378
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:1367
algorithms_suite_setup
void algorithms_suite_setup(void)
Definition: cu_algorithm.c:1812
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:1399
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:305
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:96
ASSERT_POINT4D_EQUAL
#define ASSERT_POINT4D_EQUAL(o, e, eps)
Definition: liblwgeom/cunit/cu_tester.h:144
ASSERT_INT_EQUAL
#define ASSERT_INT_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:44
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:94
ASSERT_POINT2D_EQUAL
#define ASSERT_POINT2D_EQUAL(o, e, eps)
Definition: liblwgeom/cunit/cu_tester.h:139
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:179
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:339
lwgeom_simplify
LWGEOM * lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
Simplification.
Definition: lwgeom.c:1893
LW_FALSE
#define LW_FALSE
Definition: liblwgeom.h:109
lwgeom_ndims
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:955
lwgeom_get_srid
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:927
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:242
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:2565
lwpoint_free
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:213
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1155
LW_PARSER_CHECK_NONE
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:2096
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:112
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:329
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:529
lwgeom_remove_repeated_points_in_place
int lwgeom_remove_repeated_points_in_place(LWGEOM *in, double tolerance)
Definition: lwgeom.c:1604
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:934
POINTTYPE
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:117
lwgeom_to_ewkt
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:565
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:1908
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:1246
lwgeom_cluster_kmeans
int * lwgeom_cluster_kmeans(const LWGEOM **geoms, uint32_t n, uint32_t k, double max_radius)
Take a list of LWGEOMs and a number of clusters and return an integer array indicating which cluster ...
Definition: lwkmeans.c:321
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:344
lwgeom_normalize
LWGEOM * lwgeom_normalize(const LWGEOM *geom)
Definition: liblwgeom/lwgeom_geos.c:676
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:835
LINE_MULTICROSS_END_RIGHT
@ LINE_MULTICROSS_END_RIGHT
Definition: liblwgeom.h:1670
LINE_MULTICROSS_END_SAME_FIRST_LEFT
@ LINE_MULTICROSS_END_SAME_FIRST_LEFT
Definition: liblwgeom.h:1671
LINE_MULTICROSS_END_LEFT
@ LINE_MULTICROSS_END_LEFT
Definition: liblwgeom.h:1669
LINE_CROSS_LEFT
@ LINE_CROSS_LEFT
Definition: liblwgeom.h:1667
LINE_CROSS_RIGHT
@ LINE_CROSS_RIGHT
Definition: liblwgeom.h:1668
LINE_NO_CROSS
@ LINE_NO_CROSS
Definition: liblwgeom.h:1666
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:743
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:1893
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:108
SRID_UNKNOWN
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:230
lwgeom_has_m
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:941
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:2189
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:217
SEG_COLINEAR
@ SEG_COLINEAR
Definition: liblwgeom_internal.h:218
SEG_CROSS_RIGHT
@ SEG_CROSS_RIGHT
Definition: liblwgeom_internal.h:220
SEG_CROSS_LEFT
@ SEG_CROSS_LEFT
Definition: liblwgeom_internal.h:219
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:372
GBOX::xmax
double xmax
Definition: liblwgeom.h:370
GBOX::ymin
double ymin
Definition: liblwgeom.h:371
GBOX::xmin
double xmin
Definition: liblwgeom.h:369
LWGEOM::type
uint8_t type
Definition: liblwgeom.h:477
LWLINE::points
POINTARRAY * points
Definition: liblwgeom.h:498
LWMPOINT::ngeoms
uint32_t ngeoms
Definition: liblwgeom.h:553
LWMPOINT::geoms
LWPOINT ** geoms
Definition: liblwgeom.h:548
POINT2D::y
double y
Definition: liblwgeom.h:405
POINT2D::x
double x
Definition: liblwgeom.h:405
POINT4D::m
double m
Definition: liblwgeom.h:429
POINT4D::x
double x
Definition: liblwgeom.h:429
POINT4D::z
double z
Definition: liblwgeom.h:429
POINT4D::y
double y
Definition: liblwgeom.h:429
POINTARRAY::npoints
uint32_t npoints
Definition: liblwgeom.h:442
gridspec_t::zsize
double zsize
Definition: liblwgeom.h:1390
gridspec_t::ysize
double ysize
Definition: liblwgeom.h:1389
gridspec_t::xsize
double xsize
Definition: liblwgeom.h:1388
gridspec_t::msize
double msize
Definition: liblwgeom.h:1391
gridspec_t
Snap-to-grid.
Definition: liblwgeom.h:1383
lwvarlena_t
Definition: liblwgeom.h:321
struct_lwgeom_parser_result
Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM.
Definition: liblwgeom.h:2104