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