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