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