PostGIS  2.5.7dev-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 test_lwline_clip(void)
547 {
548  LWCOLLECTION *c;
549  LWLINE *line = NULL;
550  LWLINE *l51 = NULL;
551  char *ewkt;
552 
553  /* Vertical line with vertices at y integers */
554  l51 = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
555 
556  /* Clip in the middle, mid-range. */
557  c = lwline_clip_to_ordinate_range(l51, 'Y', 1.5, 2.5);
558  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
559  //printf("c = %s\n", ewkt);
560  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
561  lwfree(ewkt);
562  lwcollection_free(c);
563 
564  /* Clip off the top. */
565  c = lwline_clip_to_ordinate_range(l51, 'Y', 3.5, 5.5);
566  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
567  //printf("c = %s\n", ewkt);
568  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 3.5,0 4))");
569  lwfree(ewkt);
570  lwcollection_free(c);
571 
572  /* Clip off the bottom. */
573  c = lwline_clip_to_ordinate_range(l51, 'Y', -1.5, 2.5);
574  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
575  //printf("c = %s\n", ewkt);
576  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 2.5))" );
577  lwfree(ewkt);
578  lwcollection_free(c);
579 
580  /* Range holds entire object. */
581  c = lwline_clip_to_ordinate_range(l51, 'Y', -1.5, 5.5);
582  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
583  //printf("c = %s\n", ewkt);
584  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))" );
585  lwfree(ewkt);
586  lwcollection_free(c);
587 
588  /* Clip on vertices. */
589  c = lwline_clip_to_ordinate_range(l51, 'Y', 1.0, 2.0);
590  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
591  //printf("c = %s\n", ewkt);
592  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1,0 2))" );
593  lwfree(ewkt);
594  lwcollection_free(c);
595 
596  /* Clip on vertices off the bottom. */
597  c = lwline_clip_to_ordinate_range(l51, 'Y', -1.0, 2.0);
598  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
599  //printf("c = %s\n", ewkt);
600  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2))" );
601  lwfree(ewkt);
602  lwcollection_free(c);
603 
604  /* Clip on top. */
605  c = lwline_clip_to_ordinate_range(l51, 'Y', -1.0, 0.0);
606  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
607  //printf("c = %s\n", ewkt);
608  CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(0 0))" );
609  lwfree(ewkt);
610  lwcollection_free(c);
611 
612  /* ST_LocateBetweenElevations(ST_GeomFromEWKT('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)'), 1, 2)) */
613  line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
614  c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 2.0);
615  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
616  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))" );
617  lwfree(ewkt);
618  lwcollection_free(c);
619  lwline_free(line);
620 
621  /* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 2)) */
622  line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
623  c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 2.0);
624  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
625  //printf("a = %s\n", ewkt);
626  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))" );
627  lwfree(ewkt);
628  lwcollection_free(c);
629  lwline_free(line);
630 
631  /* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 1)) */
632  line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
633  c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 1.0);
634  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
635  //printf("b = %s\n", ewkt);
636  CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))" );
637  lwfree(ewkt);
638  lwcollection_free(c);
639  lwline_free(line);
640 
641  /* ST_LocateBetweenElevations('LINESTRING(1 1 1, 1 2 2)', 1,1) */
642  line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1 1, 1 2 2)", LW_PARSER_CHECK_NONE);
643  c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 1.0);
644  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
645  //printf("c = %s\n", ewkt);
646  CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))" );
647  lwfree(ewkt);
648  lwcollection_free(c);
649  lwline_free(line);
650 
651  lwline_free(l51);
652 
653 }
654 
655 static void test_lwmline_clip(void)
656 {
657  LWCOLLECTION *c;
658  char *ewkt;
659  LWMLINE *mline = NULL;
660  LWLINE *line = NULL;
661 
662  /*
663  ** Set up the input line. Trivial one-member case.
664  */
665  mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
666 
667  /* Clip in the middle, mid-range. */
668  c = lwmline_clip_to_ordinate_range(mline, 'Y', 1.5, 2.5);
669  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
670  //printf("c = %s\n", ewkt);
671  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
672  lwfree(ewkt);
673  lwcollection_free(c);
674 
675  lwmline_free(mline);
676 
677  /*
678  ** Set up the input line. Two-member case.
679  */
680  mline = (LWMLINE*)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);
681 
682  /* Clip off the top. */
683  c = lwmline_clip_to_ordinate_range(mline, 'Y', 3.5, 5.5);
684  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
685  //printf("c = %s\n", ewkt);
686  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((1 3.5,1 4),(0 3.5,0 4))");
687  lwfree(ewkt);
688  lwcollection_free(c);
689 
690  lwmline_free(mline);
691 
692  /*
693  ** Set up staggered input line to create multi-type output.
694  */
695  mline = (LWMLINE*)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);
696 
697  /* Clip from 0 upwards.. */
698  c = lwmline_clip_to_ordinate_range(mline, 'Y', 0.0, 2.5);
699  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
700  //printf("c = %s\n", ewkt);
701  CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 0),LINESTRING(0 0,0 1,0 2,0 2.5))");
702  lwfree(ewkt);
703  lwcollection_free(c);
704 
705  lwmline_free(mline);
706 
707  /*
708  ** Set up input line from MAC
709  */
710  line = (LWLINE*)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)", LW_PARSER_CHECK_NONE);
711 
712  /* Clip from 3 to 3.5 */
713  c = lwline_clip_to_ordinate_range(line, 'Z', 3.0, 3.5);
714  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
715  //printf("c = %s\n", ewkt);
716  CU_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))");
717  lwfree(ewkt);
718  lwcollection_free(c);
719 
720  /* Clip from 2 to 3.5 */
721  c = lwline_clip_to_ordinate_range(line, 'Z', 2.0, 3.5);
722  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
723  //printf("c = %s\n", ewkt);
724  CU_ASSERT_STRING_EQUAL(ewkt, "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))");
725  lwfree(ewkt);
726  lwcollection_free(c);
727 
728  /* Clip from 3 to 4 */
729  c = lwline_clip_to_ordinate_range(line, 'Z', 3.0, 4.0);
730  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
731  //printf("c = %s\n", ewkt);
732  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,4 4 4 4,3 3 3 5))");
733  lwfree(ewkt);
734  lwcollection_free(c);
735 
736  /* Clip from 2 to 3 */
737  c = lwline_clip_to_ordinate_range(line, 'Z', 2.0, 3.0);
738  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
739  //printf("c = %s\n", ewkt);
740  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3),(3 3 3 5,2 2 2 6))");
741  lwfree(ewkt);
742  lwcollection_free(c);
743 
744 
745  lwline_free(line);
746 
747 }
748 
749 
750 
751 static void test_lwline_clip_big(void)
752 {
753  POINTARRAY *pa = ptarray_construct(1, 0, 3);
754  LWLINE *line = lwline_construct(SRID_UNKNOWN, NULL, pa);
755  LWCOLLECTION *c;
756  char *ewkt;
757  POINT4D p;
758 
759  p.x = 0.0;
760  p.y = 0.0;
761  p.z = 0.0;
762  ptarray_set_point4d(pa, 0, &p);
763 
764  p.x = 1.0;
765  p.y = 1.0;
766  p.z = 1.0;
767  ptarray_set_point4d(pa, 1, &p);
768 
769  p.x = 2.0;
770  p.y = 2.0;
771  p.z = 2.0;
772  ptarray_set_point4d(pa, 2, &p);
773 
774  c = lwline_clip_to_ordinate_range(line, 'Z', 0.5, 1.5);
775  ewkt = lwgeom_to_ewkt((LWGEOM*)c);
776  //printf("c = %s\n", ewkt);
777  CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0.5 0.5 0.5,1 1 1,1.5 1.5 1.5))" );
778 
779  lwfree(ewkt);
780  lwcollection_free(c);
781  lwline_free(line);
782 }
783 
784 static void test_geohash_precision(void)
785 {
786  GBOX bbox;
787  GBOX bounds;
788  int precision = 0;
789  gbox_init(&bbox);
790  gbox_init(&bounds);
791 
792  bbox.xmin = 23.0;
793  bbox.xmax = 23.0;
794  bbox.ymin = 25.2;
795  bbox.ymax = 25.2;
796  precision = lwgeom_geohash_precision(bbox, &bounds);
797  //printf("\nprecision %d\n",precision);
798  CU_ASSERT_EQUAL(precision, 20);
799 
800  bbox.xmin = 23.0;
801  bbox.ymin = 23.0;
802  bbox.xmax = 23.1;
803  bbox.ymax = 23.1;
804  precision = lwgeom_geohash_precision(bbox, &bounds);
805  //printf("precision %d\n",precision);
806  CU_ASSERT_EQUAL(precision, 3);
807 
808  bbox.xmin = 23.0;
809  bbox.ymin = 23.0;
810  bbox.xmax = 23.0001;
811  bbox.ymax = 23.0001;
812  precision = lwgeom_geohash_precision(bbox, &bounds);
813  //printf("precision %d\n",precision);
814  CU_ASSERT_EQUAL(precision, 7);
815 
816 }
817 
818 static void test_geohash_point(void)
819 {
820  char *geohash;
821 
822  geohash = geohash_point(0, 0, 16);
823  //printf("\ngeohash %s\n",geohash);
824  CU_ASSERT_STRING_EQUAL(geohash, "s000000000000000");
825  lwfree(geohash);
826 
827  geohash = geohash_point(90, 0, 16);
828  //printf("\ngeohash %s\n",geohash);
829  CU_ASSERT_STRING_EQUAL(geohash, "w000000000000000");
830  lwfree(geohash);
831 
832  geohash = geohash_point(20.012345, -20.012345, 15);
833  //printf("\ngeohash %s\n",geohash);
834  CU_ASSERT_STRING_EQUAL(geohash, "kkqnpkue9ktbpe5");
835  lwfree(geohash);
836 
837 }
838 
839 static void test_geohash(void)
840 {
841  LWPOINT *lwpoint = NULL;
842  LWLINE *lwline = NULL;
843  LWMLINE *lwmline = NULL;
844  char *geohash = NULL;
845 
846  lwpoint = (LWPOINT*)lwgeom_from_wkt("POINT(23.0 25.2)", LW_PARSER_CHECK_NONE);
847  geohash = lwgeom_geohash((LWGEOM*)lwpoint,0);
848  //printf("\ngeohash %s\n",geohash);
849  CU_ASSERT_STRING_EQUAL(geohash, "ss2r77s0du7p2ewb8hmx");
850  lwpoint_free(lwpoint);
851  lwfree(geohash);
852 
853  lwpoint = (LWPOINT*)lwgeom_from_wkt("POINT(23.0 25.2 2.0)", LW_PARSER_CHECK_NONE);
854  geohash = lwgeom_geohash((LWGEOM*)lwpoint,0);
855  //printf("geohash %s\n",geohash);
856  CU_ASSERT_STRING_EQUAL(geohash, "ss2r77s0du7p2ewb8hmx");
857  lwpoint_free(lwpoint);
858  lwfree(geohash);
859 
860  lwline = (LWLINE*)lwgeom_from_wkt("LINESTRING(23.0 23.0,23.1 23.1)", LW_PARSER_CHECK_NONE);
861  geohash = lwgeom_geohash((LWGEOM*)lwline,0);
862  //printf("geohash %s\n",geohash);
863  CU_ASSERT_STRING_EQUAL(geohash, "ss0");
864  lwline_free(lwline);
865  lwfree(geohash);
866 
867  lwline = (LWLINE*)lwgeom_from_wkt("LINESTRING(23.0 23.0,23.001 23.001)", LW_PARSER_CHECK_NONE);
868  geohash = lwgeom_geohash((LWGEOM*)lwline,0);
869  //printf("geohash %s\n",geohash);
870  CU_ASSERT_STRING_EQUAL(geohash, "ss06g7h");
871  lwline_free(lwline);
872  lwfree(geohash);
873 
874  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);
875  geohash = lwgeom_geohash((LWGEOM*)lwmline,0);
876  //printf("geohash %s\n",geohash);
877  CU_ASSERT_STRING_EQUAL(geohash, "ss0");
878  lwmline_free(lwmline);
879  lwfree(geohash);
880 }
881 
882 static void test_isclosed(void)
883 {
884  LWGEOM *geom;
885 
886  /* LINESTRING */
887 
888  /* Not Closed on 2D */
889  geom = lwgeom_from_wkt("LINESTRING(1 2,3 4)", LW_PARSER_CHECK_NONE);
890  CU_ASSERT(!lwline_is_closed((LWLINE *) geom));
891  lwgeom_free(geom);
892 
893  /* Closed on 2D */
894  geom = lwgeom_from_wkt("LINESTRING(1 2,3 4,1 2)", LW_PARSER_CHECK_NONE);
895  CU_ASSERT(lwline_is_closed((LWLINE *) geom));
896  lwgeom_free(geom);
897 
898  /* Not closed on 3D */
899  geom = lwgeom_from_wkt("LINESTRING(1 2 3,4 5 6)", LW_PARSER_CHECK_NONE);
900  CU_ASSERT(!lwline_is_closed((LWLINE *) geom));
901  lwgeom_free(geom);
902 
903  /* Closed on 3D */
904  geom = lwgeom_from_wkt("LINESTRING(1 2 3,4 5 6,1 2 3)", LW_PARSER_CHECK_NONE);
905  CU_ASSERT(lwline_is_closed((LWLINE *) geom));
906  lwgeom_free(geom);
907 
908  /* Closed on 4D, even if M is not the same */
909  geom = lwgeom_from_wkt("LINESTRING(1 2 3 4,5 6 7 8,1 2 3 0)", LW_PARSER_CHECK_NONE);
910  CU_ASSERT(lwline_is_closed((LWLINE *) geom));
911  lwgeom_free(geom);
912 
913 
914  /* CIRCULARSTRING */
915 
916  /* Not Closed on 2D */
917  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2,3 4,5 6)", LW_PARSER_CHECK_NONE);
918  CU_ASSERT(!lwcircstring_is_closed((LWCIRCSTRING *) geom));
919  lwgeom_free(geom);
920 
921  /* Closed on 2D */
922  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2,3 4,1 2)", LW_PARSER_CHECK_NONE);
923  CU_ASSERT(lwcircstring_is_closed((LWCIRCSTRING *) geom));
924  lwgeom_free(geom);
925 
926  /* Not closed on 3D */
927  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2 3,4 5 6,7 8 9)", LW_PARSER_CHECK_NONE);
928  CU_ASSERT(!lwcircstring_is_closed((LWCIRCSTRING *) geom));
929  lwgeom_free(geom);
930 
931  /* Closed on 3D */
932  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2 3,4 5 6,1 2 3)", LW_PARSER_CHECK_NONE);
933  CU_ASSERT(lwcircstring_is_closed((LWCIRCSTRING *) geom));
934  lwgeom_free(geom);
935 
936  /* Closed on 4D, even if M is not the same */
937  geom = lwgeom_from_wkt("CIRCULARSTRING(1 2 3 4,5 6 7 8,1 2 3 0)", LW_PARSER_CHECK_NONE);
938  CU_ASSERT(lwcircstring_is_closed((LWCIRCSTRING *) geom));
939  lwgeom_free(geom);
940 
941 
942  /* COMPOUNDCURVE */
943 
944  /* Not Closed on 2D */
945  geom = lwgeom_from_wkt("COMPOUNDCURVE(CIRCULARSTRING(1 2,3 4,1 2),(1 2,7 8,5 6))", LW_PARSER_CHECK_NONE);
946  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
947  lwgeom_free(geom);
948 
949  geom = lwgeom_from_wkt("COMPOUNDCURVE((1 2,3 4,1 2),CIRCULARSTRING(1 2,7 8,5 6))", LW_PARSER_CHECK_NONE);
950  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
951  lwgeom_free(geom);
952 
953  /* Closed on 2D */
954  geom = lwgeom_from_wkt("COMPOUNDCURVE(CIRCULARSTRING(1 2,3 4,5 6), (5 6,7 8,1 2))", LW_PARSER_CHECK_NONE);
955  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
956  lwgeom_free(geom);
957 
958  geom = lwgeom_from_wkt("COMPOUNDCURVE((1 2,3 4,5 6),CIRCULARSTRING(5 6,7 8,1 2))", LW_PARSER_CHECK_NONE);
959  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
960  lwgeom_free(geom);
961 
962  /* Not Closed on 3D */
963  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);
964  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
965  lwgeom_free(geom);
966 
967  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);
968  CU_ASSERT(!lwcompound_is_closed((LWCOMPOUND *) geom));
969  lwgeom_free(geom);
970 
971  /* Closed on 3D */
972  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);
973  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
974  lwgeom_free(geom);
975 
976  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);
977  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
978  lwgeom_free(geom);
979 
980  /* Closed on 4D, even if M is not the same */
981  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);
982  CU_ASSERT(lwcompound_is_closed((LWCOMPOUND *) geom));
983  lwgeom_free(geom);
984 }
985 
986 
987 static void test_geohash_point_as_int(void)
988 {
989  unsigned int gh;
990  POINT2D p;
991  unsigned long long rs;
992 
993  p.x = 50; p.y = 35;
994  gh = geohash_point_as_int(&p);
995  rs = 3440103613;
996  CU_ASSERT_EQUAL(gh, rs);
997  p.x = 140; p.y = 45;
998  gh = geohash_point_as_int(&p);
999  rs = 3982480893;
1000  CU_ASSERT_EQUAL(gh, rs);
1001  p.x = 140; p.y = 55;
1002  gh = geohash_point_as_int(&p);
1003  rs = 4166944232;
1004  CU_ASSERT_EQUAL(gh, rs);
1005 }
1006 
1007 static void test_lwgeom_remove_repeated_points(void)
1008 {
1009  LWGEOM *g;
1010  char *ewkt;
1011 
1012  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);
1013  lwgeom_remove_repeated_points_in_place(g, 1);
1014  ewkt = lwgeom_to_ewkt(g);
1015  CU_ASSERT_STRING_EQUAL(ewkt, "MULTIPOINT(0 0,10 0,10 10,0 10,5 5)");
1016  lwgeom_free(g);
1017  lwfree(ewkt);
1018 
1019  g = lwgeom_from_wkt("LINESTRING(1612830.15445 4841287.12672,1612830.15824 4841287.12674,1612829.98813 4841274.56198)", LW_PARSER_CHECK_NONE);
1020  lwgeom_remove_repeated_points_in_place(g, 0.01);
1021  ewkt = lwgeom_to_ewkt(g);
1022  CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(1612830.15445 4841287.12672,1612829.98813 4841274.56198)");
1023  lwgeom_free(g);
1024  lwfree(ewkt);
1025 
1026  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);
1027  lwgeom_remove_repeated_points_in_place(g, 1);
1028  ewkt = lwgeom_to_ewkt(g);
1029  CU_ASSERT_STRING_EQUAL(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)");
1030  lwgeom_free(g);
1031  lwfree(ewkt);
1032 
1033  g = lwgeom_from_wkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", LW_PARSER_CHECK_NONE);
1034  lwgeom_remove_repeated_points_in_place(g, 10000);
1035  ewkt = lwgeom_to_ewkt(g);
1036  CU_ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 1,1 0,0 0))");
1037  lwgeom_free(g);
1038  lwfree(ewkt);
1039 }
1040 
1041 static void
1042 test_geohash_bbox(void)
1043 {
1044  double lat[2], lon[2];
1045 
1046  /* SELECT ST_GeoHash(ST_SetSRID(ST_MakePoint(-126,48),4326)) */
1047  decode_geohash_bbox("c0w3hf1s70w3hf1s70w3", lat, lon, 100);
1048  CU_ASSERT_DOUBLE_EQUAL(lat[0], 48, 1e-11);
1049  CU_ASSERT_DOUBLE_EQUAL(lat[1], 48, 1e-11);
1050  CU_ASSERT_DOUBLE_EQUAL(lon[0], -126, 1e-11);
1051  CU_ASSERT_DOUBLE_EQUAL(lon[1], -126, 1e-11);
1052 
1053  cu_error_msg_reset();
1054  decode_geohash_bbox("@@@@@@", lat, lon, 100);
1055  ASSERT_STRING_EQUAL(cu_error_msg, "decode_geohash_bbox: Invalid character '@'");
1056 }
1057 
1058 static void test_lwgeom_simplify(void)
1059 {
1060  LWGEOM *l;
1061  LWGEOM *g;
1062  char *ewkt;
1063 
1064  /* Simplify but only so far... */
1065  g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
1066  l = lwgeom_simplify(g, 10, LW_TRUE);
1067  ewkt = lwgeom_to_ewkt(l);
1068  CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)");
1069  lwgeom_free(g);
1070  lwgeom_free(l);
1071  lwfree(ewkt);
1072 
1073  /* Simplify but only so far... */
1074  g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
1075  l = lwgeom_simplify(g, 10, LW_TRUE);
1076  ewkt = lwgeom_to_ewkt(l);
1077  CU_ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))");
1078  lwgeom_free(g);
1079  lwgeom_free(l);
1080  lwfree(ewkt);
1081 
1082  /* Simplify and collapse */
1083  g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
1084  l = lwgeom_simplify(g, 10, LW_FALSE);
1085  CU_ASSERT_EQUAL(l, NULL);
1086  lwgeom_free(g);
1087  lwgeom_free(l);
1088 
1089  /* Simplify and collapse */
1090  g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
1091  l = lwgeom_simplify(g, 10, LW_FALSE);
1092  CU_ASSERT_EQUAL(l, NULL);
1093  lwgeom_free(g);
1094  lwgeom_free(l);
1095 
1096  /* Not simplifiable */
1097  g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE);
1098  l = lwgeom_simplify(g, 1.0, LW_FALSE);
1099  ewkt = lwgeom_to_ewkt(l);
1100  CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)");
1101  lwgeom_free(g);
1102  lwgeom_free(l);
1103  lwfree(ewkt);
1104 
1105  /* Simplifiable */
1106  g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE);
1107  l = lwgeom_simplify(g, 1.0, LW_FALSE);
1108  ewkt = lwgeom_to_ewkt(l);
1109  CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)");
1110  lwgeom_free(g);
1111  lwgeom_free(l);
1112  lwfree(ewkt);
1113 }
1114 
1115 
1116 static void do_median_dims_check(char* wkt, int expected_dims)
1117 {
1118  LWGEOM* g = lwgeom_from_wkt(wkt, LW_PARSER_CHECK_NONE);
1119  LWPOINT* result = lwgeom_median(g, 1e-8, 100, LW_FALSE);
1120 
1121  CU_ASSERT_EQUAL(expected_dims, lwgeom_ndims((LWGEOM*) result));
1122 
1123  lwgeom_free(g);
1124  lwpoint_free(result);
1125 }
1126 
1127 static void test_median_handles_3d_correctly(void)
1128 {
1129  do_median_dims_check("MULTIPOINT ((1 3), (4 7), (2 9), (0 4), (2 2))", 2);
1130  do_median_dims_check("MULTIPOINT Z ((1 3 4), (4 7 8), (2 9 1), (0 4 4), (2 2 3))", 3);
1131  do_median_dims_check("MULTIPOINT M ((1 3 4), (4 7 8), (2 9 1), (0 4 4), (2 2 3))", 2);
1132  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);
1133 }
1134 
1135 static double
1136 test_weighted_distance(const POINT4D* curr, const POINT4D* points, uint32_t npoints)
1137 {
1138  double distance = 0.0;
1139  uint32_t i;
1140  for (i = 0; i < npoints; i++)
1141  {
1142  distance += distance3d_pt_pt((POINT3D*)curr, (POINT3D*)&points[i]) * points[i].m;
1143  }
1144 
1145  return distance;
1146 }
1147 
1148 static void do_median_test(char* input, char* expected, int fail_if_not_converged, int iter_count)
1149 {
1150  cu_error_msg_reset();
1151  LWGEOM* g = lwgeom_from_wkt(input, LW_PARSER_CHECK_NONE);
1152  LWPOINT* expected_result = NULL;
1153  POINT4D actual_pt;
1154  POINT4D expected_pt;
1155  const double tolerance = FP_TOLERANCE / 10.0;
1156 
1157  LWPOINT* result = lwgeom_median(g, tolerance, iter_count, fail_if_not_converged);
1158  int passed = LW_FALSE;
1159 
1160  if (expected != NULL)
1161  {
1162  expected_result = lwgeom_as_lwpoint(lwgeom_from_wkt(expected, LW_PARSER_CHECK_NONE));
1163  lwpoint_getPoint4d_p(expected_result, &expected_pt);
1164  }
1165  if (result != NULL)
1166  {
1167  lwpoint_getPoint4d_p(result, &actual_pt);
1168  }
1169 
1170  if (result != NULL && expected != NULL) /* got something, expecting something */
1171  {
1172  passed = LW_TRUE;
1173  passed = passed && lwgeom_is_empty((LWGEOM*) expected_result) == lwgeom_is_empty((LWGEOM*) result);
1174  passed = passed && (lwgeom_has_z((LWGEOM*) expected_result) == lwgeom_has_z((LWGEOM*) result));
1175 
1176  if (passed && !lwgeom_is_empty((LWGEOM*) result))
1177  {
1178  if (g->type == POINTTYPE)
1179  {
1180  passed &= fabs(actual_pt.x - expected_pt.x) < tolerance;
1181  passed &= fabs(actual_pt.y - expected_pt.y) < tolerance;
1182  passed &= (!lwgeom_has_z((LWGEOM*) expected_result) || fabs(actual_pt.z - expected_pt.z) < tolerance);
1183  passed &= (!lwgeom_has_m((LWGEOM*) expected_result) || fabs(actual_pt.m - expected_pt.m) < tolerance);
1184  }
1185  else
1186  {
1187  /* Check that the difference between the obtained geometric
1188  median and the expected point is within tolerance */
1189  uint32_t npoints = 1;
1190  int input_empty = LW_TRUE;
1191  POINT4D* points = lwmpoint_extract_points_4d(lwgeom_as_lwmpoint(g), &npoints, &input_empty);
1192  double distance_expected = test_weighted_distance(&expected_pt, points, npoints);
1193  double distance_result = test_weighted_distance(&actual_pt, points, npoints);
1194 
1195  passed = distance_result <= (1.0 + tolerance) * distance_expected;
1196  if (!passed)
1197  {
1198  printf("Diff: Got %.10f Expected %.10f\n", distance_result, distance_expected);
1199  }
1200  lwfree(points);
1201  }
1202  }
1203 
1204  if (!passed)
1205  {
1206  printf("median_test input %s (parsed %s) expected %s got %s\n",
1207  input, lwgeom_to_ewkt(g),
1208  lwgeom_to_ewkt((LWGEOM*) expected_result),
1209  lwgeom_to_ewkt((LWGEOM*) result));
1210  }
1211 
1212  }
1213  else if (result == NULL && expected == NULL) /* got nothing, expecting nothing */
1214  {
1215  passed = LW_TRUE;
1216  }
1217  else if (result != NULL && expected == NULL) /* got something, expecting nothing */
1218  {
1219  passed = LW_FALSE;
1220  printf("median_test input %s (parsed %s) expected NULL got %s\n", input, lwgeom_to_ewkt(g), lwgeom_to_ewkt((LWGEOM*) result));
1221  }
1222  else if (result == NULL && expected != NULL) /* got nothing, expecting something */
1223  {
1224  passed = LW_FALSE;
1225  printf("%s", cu_error_msg);
1226  printf("median_test input %s (parsed %s) expected %s got NULL\n", input, lwgeom_to_ewkt(g), lwgeom_to_ewkt((LWGEOM*) expected_result));
1227  }
1228 
1229  CU_ASSERT_TRUE(passed);
1230 
1231  lwgeom_free(g);
1232  lwpoint_free(expected_result);
1233  lwpoint_free(result);
1234 }
1235 
1236 static void test_median_robustness(void)
1237 {
1238  /* A simple implementation of Weiszfeld's algorithm will fail if the median is equal
1239  * to any one of the inputs, during any iteration of the algorithm.
1240  *
1241  * Because the algorithm uses the centroid as a starting point, this situation will
1242  * occur in the test case below.
1243  */
1244  do_median_test("MULTIPOINT ((0 -1), (0 0), (0 1))", "POINT (0 0)", LW_TRUE, 1000);
1245 
1246  /* Same as above but 3D, and shifter */
1247  do_median_test("MULTIPOINT ((1 -1 3), (1 0 2), (2 1 1))", "POINT (1 0 2)", LW_TRUE, 1000);
1248 
1249  /* Starting point is duplicated */
1250  do_median_test("MULTIPOINT ((0 -1), (0 0), (0 0), (0 1))", "POINT (0 0)", LW_TRUE, 1000);
1251 
1252  /* Cube */
1253  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))",
1254  "POINT (15 15 15)", LW_TRUE, 1000);
1255 
1256  /* Some edge cases */
1257  do_median_test("POINT (7 6)", "POINT (7 6)", LW_TRUE, 1000);
1258  do_median_test("POINT (7 6 2)", "POINT (7 6 2)", LW_TRUE, 1000);
1259  do_median_test("MULTIPOINT ((7 6 2), EMPTY)", "POINT (7 6 2)", LW_TRUE, 1000);
1260 
1261  /* Empty input */
1262  do_median_test("MULTIPOINT EMPTY", "POINT EMPTY", LW_FALSE, 1000);
1263  do_median_test("MULTIPOINT (EMPTY)", "POINT EMPTY", LW_FALSE, 1000);
1264  do_median_test("MULTIPOINT EMPTY", "POINT EMPTY", LW_TRUE, 1000);
1265  do_median_test("MULTIPOINT (EMPTY)", "POINT EMPTY", LW_TRUE, 1000);
1266  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);
1267 
1268  /* Weighted input */
1269  do_median_test("MULTIPOINT ZM (1 -1 3 1, 1 0 2 7, 2 1 1 1)", "POINT (1 0 2)", LW_TRUE, 1000);
1270  do_median_test("MULTIPOINT ZM (-1 1 -3 1, -1 0 -2 7, -2 -1 -1 1)", "POINT (-1 0 -2)", LW_TRUE, 1000);
1271  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);
1272 
1273  /* Point that is replaced by two half-weighted */
1274  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);
1275  /* Point is doubled and then erased by negative weight */
1276  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);
1277  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);
1278  /* Weightless input won't converge */
1279  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);
1280  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);
1281  /* Negative weight won't converge */
1282  do_median_test("MULTIPOINT ZM ((0 -1 0 -1), (0 0 0 -1), (0 1 0 -1))", NULL, LW_FALSE, 1000);
1283  do_median_test("MULTIPOINT ZM ((0 -1 0 -1), (0 0 0 -1), (0 1 0 -1))", NULL, LW_TRUE, 1000);
1284 
1285  /* Bind convergence too tightly */
1286  do_median_test("MULTIPOINT ((0 0), (1 1), (0 1), (2 2))", "POINT(0.75 1.0)", LW_FALSE, 0);
1287  do_median_test("MULTIPOINT ((0 0), (1 1), (0 1), (2 2))", NULL, LW_TRUE, 1);
1288  /* Unsupported geometry type */
1289  do_median_test("POLYGON((1 0,0 1,1 2,2 1,1 0))", NULL, LW_TRUE, 1000);
1290  do_median_test("POLYGON((1 0,0 1,1 2,2 1,1 0))", NULL, LW_FALSE, 1000);
1291 
1292  /* Median point is included */
1293  do_median_test("MULTIPOINT ZM ("
1294  "(1480 0 200 100),"
1295  "(620 0 200 100),"
1296  "(1000 0 -200 100),"
1297  "(1000 0 -590 100),"
1298  "(1025 0 65 100),"
1299  "(1025 0 -65 100)"
1300  ")",
1301  "POINT (1025 0 -65)", LW_TRUE, 10000);
1302 
1303 #if 0
1304  /* Leads to invalid result (0 0 0) with 80bit (fmulp + faddp) precision. ok with 64 bit float ops */
1305  do_median_test("MULTIPOINT ZM ("
1306  "(0 0 20000 0.5),"
1307  "(0 0 59000 0.5),"
1308  "(0 -3000 -3472.22222222222262644208967685699462890625 1),"
1309  "(0 3000 3472.22222222222262644208967685699462890625 1),"
1310  "(0 0 -1644.736842105263121993630193173885345458984375 1),"
1311  "(0 0 1644.736842105263121993630193173885345458984375 1),"
1312  "(0 48000 -20000 1.3),"
1313  "(0 -48000 -20000 1.3)"
1314  ")",
1315  "POINT (0 0 0)", LW_TRUE, 10000);
1316 #endif
1317 
1318 #if 0
1319  /* Leads to invalid result (0 0 0) with 64bit (vfmadd231sd) precision. Ok with 80 bit float ops */
1320  do_median_test("MULTIPOINT ZM ("
1321  "(0 0 20000 0.5),"
1322  "(0 0 59000 0.5),"
1323  "(0 -3000 -3472.22222222222262644208967685699462890625 1),"
1324  "(0 3000 3472.22222222222262644208967685699462890625 1),"
1325  "(0 -0.00000000000028047739569477638384522295466033823196 -1644.736842105263121993630193173885345458984375 1),"
1326  "(0 0.00000000000028047739569477638384522295466033823196 1644.736842105263121993630193173885345458984375 1),"
1327  "(0 48000 -20000 1.3),"
1328  "(0 -48000 -20000 1.3)"
1329  ")",
1330  "POINT (0 0 0)", LW_TRUE, 10000);
1331 #endif
1332 }
1333 
1334 static void test_point_density(void)
1335 {
1336  LWGEOM *geom;
1337  LWMPOINT *mpt;
1338  // char *ewkt;
1339 
1340  /* POLYGON */
1341  geom = lwgeom_from_wkt("POLYGON((1 0,0 1,1 2,2 1,1 0))", LW_PARSER_CHECK_NONE);
1342  mpt = lwgeom_to_points(geom, 100);
1343  CU_ASSERT_EQUAL(mpt->ngeoms,100);
1344  // ewkt = lwgeom_to_ewkt((LWGEOM*)mpt);
1345  // printf("%s\n", ewkt);
1346  // lwfree(ewkt);
1347  lwmpoint_free(mpt);
1348 
1349  mpt = lwgeom_to_points(geom, 1);
1350  CU_ASSERT_EQUAL(mpt->ngeoms,1);
1351  lwmpoint_free(mpt);
1352 
1353  mpt = lwgeom_to_points(geom, 0);
1354  CU_ASSERT_EQUAL(mpt, NULL);
1355  lwmpoint_free(mpt);
1356 
1357  lwgeom_free(geom);
1358 
1359  /* MULTIPOLYGON */
1360  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);
1361 
1362  mpt = lwgeom_to_points(geom, 1000);
1363  CU_ASSERT_EQUAL(mpt->ngeoms,1000);
1364  lwmpoint_free(mpt);
1365 
1366  mpt = lwgeom_to_points(geom, 1);
1367  CU_ASSERT_EQUAL(mpt->ngeoms,1);
1368  lwmpoint_free(mpt);
1369 
1370  lwgeom_free(geom);
1371 }
1372 
1373 static void test_lwpoly_construct_circle(void)
1374 {
1375  LWPOLY* p;
1376  const GBOX* g;
1377  const int srid = 4326;
1378  const uint32_t segments_per_quad = 17;
1379  const int x = 10;
1380  const int y = 20;
1381  const int r = 5;
1382 
1383  /* With normal arguments you should get something circle-y */
1384  p = lwpoly_construct_circle(srid, x, y, r, segments_per_quad, LW_TRUE);
1385 
1386  ASSERT_INT_EQUAL(lwgeom_count_vertices(lwpoly_as_lwgeom(p)), segments_per_quad * 4 + 1);
1387  ASSERT_INT_EQUAL(lwgeom_get_srid(lwpoly_as_lwgeom(p)), srid)
1388 
1389  g = lwgeom_get_bbox(lwpoly_as_lwgeom(p));
1390  CU_ASSERT_DOUBLE_EQUAL(g->xmin, x-r, 0.1);
1391  CU_ASSERT_DOUBLE_EQUAL(g->xmax, x+r, 0.1);
1392  CU_ASSERT_DOUBLE_EQUAL(g->ymin, y-r, 0.1);
1393  CU_ASSERT_DOUBLE_EQUAL(g->ymax, y+r, 0.1);
1394 
1395  CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(lwpoly_as_lwgeom(p)), M_PI*5*5, 0.1);
1396 
1397  lwpoly_free(p);
1398 
1399  /* No segments? No circle. */
1400  p = lwpoly_construct_circle(srid, x, y, r, 0, LW_TRUE);
1401  CU_ASSERT_TRUE(p == NULL);
1402 
1403  /* Negative radius? No circle. */
1404  p = lwpoly_construct_circle(srid, x, y, -1e-3, segments_per_quad, LW_TRUE);
1405  CU_ASSERT_TRUE(p == NULL);
1406 
1407  /* Zero radius? Invalid circle */
1408  p = lwpoly_construct_circle(srid, x, y, 0, segments_per_quad, LW_TRUE);
1409  CU_ASSERT_TRUE(p != NULL);
1410  lwpoly_free(p);
1411 }
1412 
1413 static void test_kmeans(void)
1414 {
1415  static int cluster_size = 100;
1416  static int num_clusters = 4;
1417  static double spread = 1.5;
1418  int N = cluster_size * num_clusters;
1419  LWGEOM **geoms;
1420  int i, j, k=0;
1421  int *r;
1422 
1423  geoms = lwalloc(sizeof(LWGEOM*) * N);
1424 
1425  for (j = 0; j < num_clusters; j++) {
1426  for (i = 0; i < cluster_size; i++)
1427  {
1428  double u1 = 1.0 * rand() / RAND_MAX;
1429  double u2 = 1.0 * rand() / RAND_MAX;
1430  double z1 = spread * j + sqrt(-2*log2(u1))*cos(2*M_PI*u2);
1431  double z2 = spread * j + sqrt(-2*log2(u1))*sin(2*M_PI*u2);
1432 
1433  LWPOINT *lwp = lwpoint_make2d(SRID_UNKNOWN, z1, z2);
1434  geoms[k++] = lwpoint_as_lwgeom(lwp);
1435  }
1436  }
1437 
1438  r = lwgeom_cluster_2d_kmeans((const LWGEOM **)geoms, N, num_clusters);
1439 
1440  // for (i = 0; i < k; i++)
1441  // {
1442  // printf("[%d] %s\n", r[i], lwgeom_to_ewkt(geoms[i]));
1443  // }
1444 
1445  /* Clean up */
1446  lwfree(r);
1447  for (i = 0; i < k; i++)
1448  lwgeom_free(geoms[i]);
1449  lwfree(geoms);
1450 
1451  return;
1452 }
1453 
1454 static void test_trim_bits(void)
1455 {
1456  POINTARRAY *pta = ptarray_construct_empty(LW_TRUE, LW_TRUE, 2);
1457  POINT4D pt;
1458  LWLINE *line;
1459  int precision;
1460  uint32_t i;
1461 
1462  pt.x = 1.2345678901234;
1463  pt.y = 2.3456789012345;
1464  pt.z = 3.4567890123456;
1465  pt.m = 4.5678901234567;
1466 
1467  ptarray_insert_point(pta, &pt, 0);
1468 
1469  pt.x *= 3;
1470  pt.y *= 3;
1471  pt.y *= 3;
1472  pt.z *= 3;
1473 
1474  ptarray_insert_point(pta, &pt, 1);
1475 
1476  line = lwline_construct(0, NULL, pta);
1477 
1478  for (precision = -15; precision <= 15; precision++)
1479  {
1480  LWLINE *line2 = (LWLINE*) lwgeom_clone_deep((LWGEOM*) line);
1481  lwgeom_trim_bits_in_place((LWGEOM*) line2, precision, precision, precision, precision);
1482 
1483  for (i = 0; i < line->points->npoints; i++)
1484  {
1485  POINT4D pt1 = getPoint4d(line->points, i);
1486  POINT4D pt2 = getPoint4d(line2->points, i);
1487 
1488  CU_ASSERT_DOUBLE_EQUAL(pt1.x, pt2.x, pow(10, -1*precision));
1489  CU_ASSERT_DOUBLE_EQUAL(pt1.y, pt2.y, pow(10, -1*precision));
1490  CU_ASSERT_DOUBLE_EQUAL(pt1.z, pt2.z, pow(10, -1*precision));
1491  CU_ASSERT_DOUBLE_EQUAL(pt1.m, pt2.m, pow(10, -1*precision));
1492  }
1493 
1494  lwline_free(line2);
1495  }
1496 
1497  lwline_free(line);
1498 }
1499 
1500 /*
1501 ** Used by test harness to register the tests in this file.
1502 */
1503 void algorithms_suite_setup(void);
1504 void algorithms_suite_setup(void)
1505 {
1506  CU_pSuite suite = CU_add_suite("algorithm", init_cg_suite, clean_cg_suite);
1507  PG_ADD_TEST(suite,test_lw_segment_side);
1508  PG_ADD_TEST(suite,test_lw_segment_intersects);
1509  PG_ADD_TEST(suite,test_lwline_crossing_short_lines);
1510  PG_ADD_TEST(suite,test_lwline_crossing_long_lines);
1511  PG_ADD_TEST(suite,test_lwline_crossing_bugs);
1512  PG_ADD_TEST(suite,test_lwpoint_set_ordinate);
1513  PG_ADD_TEST(suite,test_lwpoint_get_ordinate);
1514  PG_ADD_TEST(suite,test_point_interpolate);
1515  PG_ADD_TEST(suite,test_lwline_interpolate_points);
1516  PG_ADD_TEST(suite,test_lwline_clip);
1517  PG_ADD_TEST(suite,test_lwline_clip_big);
1518  PG_ADD_TEST(suite,test_lwmline_clip);
1519  PG_ADD_TEST(suite,test_geohash_point);
1520  PG_ADD_TEST(suite,test_geohash_precision);
1521  PG_ADD_TEST(suite,test_geohash);
1522  PG_ADD_TEST(suite,test_geohash_point_as_int);
1523  PG_ADD_TEST(suite, test_geohash_bbox);
1524  PG_ADD_TEST(suite,test_isclosed);
1525  PG_ADD_TEST(suite,test_lwgeom_simplify);
1526  PG_ADD_TEST(suite,test_lw_arc_center);
1527  PG_ADD_TEST(suite,test_point_density);
1528  PG_ADD_TEST(suite,test_kmeans);
1529  PG_ADD_TEST(suite,test_median_handles_3d_correctly);
1530  PG_ADD_TEST(suite,test_median_robustness);
1531  PG_ADD_TEST(suite,test_lwpoly_construct_circle);
1532  PG_ADD_TEST(suite,test_trim_bits);
1533  PG_ADD_TEST(suite,test_lwgeom_remove_repeated_points);
1534 }
test_kmeans
static void test_kmeans(void)
Definition: cu_algorithm.c:1413
test_geohash
static void test_geohash(void)
Definition: cu_algorithm.c:839
test_geohash_bbox
static void test_geohash_bbox(void)
Definition: cu_algorithm.c:1042
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:751
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:818
setpoint
#define setpoint(p, x1, y1)
test_isclosed
static void test_isclosed(void)
Definition: cu_algorithm.c:882
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:1007
test_point_density
static void test_point_density(void)
Definition: cu_algorithm.c:1334
test_lwpoly_construct_circle
static void test_lwpoly_construct_circle(void)
Definition: cu_algorithm.c:1373
test_median_robustness
static void test_median_robustness(void)
Definition: cu_algorithm.c:1236
test_weighted_distance
static double test_weighted_distance(const POINT4D *curr, const POINT4D *points, uint32_t npoints)
Definition: cu_algorithm.c:1136
test_lwgeom_simplify
static void test_lwgeom_simplify(void)
Definition: cu_algorithm.c:1058
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:784
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:655
test_trim_bits
static void test_trim_bits(void)
Definition: cu_algorithm.c:1454
test_geohash_point_as_int
static void test_geohash_point_as_int(void)
Definition: cu_algorithm.c:987
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:1127
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:1116
algorithms_suite_setup
void algorithms_suite_setup(void)
Definition: cu_algorithm.c:1504
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:1148
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:546
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: g_box.c:47
cu_error_msg_reset
void cu_error_msg_reset()
Definition: liblwgeom/cunit/cu_tester.c:301
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
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:170
lwpoint_make2d
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
getPoint4d
POINT4D getPoint4d(const POINTARRAY *pa, uint32_t n)
Definition: lwgeom_api.c:106
lwline_as_lwgeom
LWGEOM * lwline_as_lwgeom(const LWLINE *obj)
Definition: lwgeom.c:330
lwgeom_simplify
LWGEOM * lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
Definition: lwgeom.c:1857
LW_FALSE
#define LW_FALSE
Definition: liblwgeom.h:77
lwgeom_ndims
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:944
lwgeom_get_srid
int32_t lwgeom_get_srid(const LWGEOM *geom)
Return SRID number.
Definition: lwgeom.c:916
lwmpoint_free
void lwmpoint_free(LWMPOINT *mpt)
Definition: lwmpoint.c:72
lwgeom_as_lwmpoint
LWMPOINT * lwgeom_as_lwmpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:233
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:2518
lwpoint_free
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:213
lwgeom_free
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1144
LW_PARSER_CHECK_NONE
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:2005
lwgeom_remove_repeated_points_in_place
void lwgeom_remove_repeated_points_in_place(LWGEOM *in, double tolerance)
Definition: lwgeom.c:1603
LW_SUCCESS
#define LW_SUCCESS
Definition: liblwgeom.h:80
geohash_point_as_int
unsigned int geohash_point_as_int(POINT2D *pt)
Definition: lwalgorithm.c:657
lwpoly_as_lwgeom
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:320
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:461
lwgeom_clone_deep
LWGEOM * lwgeom_clone_deep(const LWGEOM *lwgeom)
Deep clone an LWGEOM, everything is copied.
Definition: lwgeom.c:520
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:62
lwpoly_construct_circle
LWPOLY * lwpoly_construct_circle(int srid, double x, double y, double radius, uint32_t segments_per_quarter, char exterior)
Definition: lwpoly.c:120
lwgeom_has_z
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:930
POINTTYPE
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:85
lwgeom_to_ewkt
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:556
lwgeom_as_lwpoint
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:161
lwgeom_area
double lwgeom_area(const LWGEOM *geom)
Definition: lwgeom.c:1872
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:96
lwgeom_count_vertices
uint32_t lwgeom_count_vertices(const LWGEOM *geom)
Count the total number of vertices in any LWGEOM.
Definition: lwgeom.c:1235
lwgeom_geohash
char * lwgeom_geohash(const LWGEOM *lwgeom, int precision)
Calculate the GeoHash (http://geohash.org) string for a geometry.
Definition: lwalgorithm.c:856
distance3d_pt_pt
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:815
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:542
lwgeom_to_points
LWMPOINT * lwgeom_to_points(const LWGEOM *lwgeom, uint32_t npoints)
Definition: liblwgeom/lwgeom_geos.c:1705
lwfree
void lwfree(void *mem)
Definition: lwutil.c:244
lwpoint_as_lwgeom
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:335
LINE_MULTICROSS_END_RIGHT
@ LINE_MULTICROSS_END_RIGHT
Definition: liblwgeom.h:1522
LINE_MULTICROSS_END_SAME_FIRST_LEFT
@ LINE_MULTICROSS_END_SAME_FIRST_LEFT
Definition: liblwgeom.h:1523
LINE_MULTICROSS_END_LEFT
@ LINE_MULTICROSS_END_LEFT
Definition: liblwgeom.h:1521
LINE_CROSS_LEFT
@ LINE_CROSS_LEFT
Definition: liblwgeom.h:1519
LINE_CROSS_RIGHT
@ LINE_CROSS_RIGHT
Definition: liblwgeom.h:1520
LINE_NO_CROSS
@ LINE_NO_CROSS
Definition: liblwgeom.h:1518
ptarray_construct_empty
POINTARRAY * ptarray_construct_empty(char hasz, char hasm, uint32_t maxpoints)
Create a new POINTARRAY with no points.
Definition: ptarray.c:70
lwcollection_free
void lwcollection_free(LWCOLLECTION *col)
Definition: lwcollection.c:356
lwgeom_is_empty
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwgeom.c:1393
ptarray_free
void ptarray_free(POINTARRAY *pa)
Definition: ptarray.c:328
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:734
lwgeom_from_wkt
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition: lwin_wkt.c:904
lwalloc
void * lwalloc(size_t size)
Definition: lwutil.c:229
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
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:76
SRID_UNKNOWN
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:188
lwgeom_has_m
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:937
lwline_construct
LWLINE * lwline_construct(int srid, GBOX *bbox, POINTARRAY *points)
Definition: lwline.c:42
ptarray_set_point4d
void ptarray_set_point4d(POINTARRAY *pa, uint32_t n, const POINT4D *p4d)
Definition: lwgeom_api.c:435
lwline_free
void lwline_free(LWLINE *line)
Definition: lwline.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:240
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:760
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:251
lwline_clip_to_ordinate_range
LWCOLLECTION * lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to)
Clip a line based on the from/to range of one of its ordinates.
Definition: lwlinearreferencing.c:547
lwline_is_closed
int lwline_is_closed(const LWLINE *line)
Definition: lwline.c:454
lwcompound_is_closed
int lwcompound_is_closed(const LWCOMPOUND *curve)
Definition: lwcompound.c:35
SEG_NO_INTERSECTION
@ SEG_NO_INTERSECTION
Definition: liblwgeom_internal.h:212
SEG_COLINEAR
@ SEG_COLINEAR
Definition: liblwgeom_internal.h:213
SEG_CROSS_RIGHT
@ SEG_CROSS_RIGHT
Definition: liblwgeom_internal.h:215
SEG_CROSS_LEFT
@ SEG_CROSS_LEFT
Definition: liblwgeom_internal.h:214
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:372
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:228
FP_TOLERANCE
#define FP_TOLERANCE
Floating point comparators.
Definition: liblwgeom_internal.h:53
ptarray_has_z
int ptarray_has_z(const POINTARRAY *pa)
Definition: ptarray.c:36
lw_segment_side
int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
lw_segment_side()
Definition: lwalgorithm.c:64
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:321
decode_geohash_bbox
void decode_geohash_bbox(char *geohash, double *lat, double *lon, int precision)
Definition: lwalgorithm.c:714
ptarray_has_m
int ptarray_has_m(const POINTARRAY *pa)
Definition: ptarray.c:43
lwmline_clip_to_ordinate_range
LWCOLLECTION * lwmline_clip_to_ordinate_range(const LWMLINE *mline, char ordinate, double from, double to)
Clip a multi-line based on the from/to range of one of its ordinates.
Definition: lwlinearreferencing.c:468
geohash_point
char * geohash_point(double longitude, double latitude, int precision)
Definition: lwalgorithm.c:591
liblwgeom_internal.h
distance
Datum distance(PG_FUNCTION_ARGS)
Definition: lwgeom_backend_api.c:199
cu_tester.h
GBOX::ymax
double ymax
Definition: liblwgeom.h:298
GBOX::xmax
double xmax
Definition: liblwgeom.h:296
GBOX::ymin
double ymin
Definition: liblwgeom.h:297
GBOX::xmin
double xmin
Definition: liblwgeom.h:295
LWGEOM::type
uint8_t type
Definition: liblwgeom.h:399
LWLINE::points
POINTARRAY * points
Definition: liblwgeom.h:425
LWMPOINT::ngeoms
uint32_t ngeoms
Definition: liblwgeom.h:471
POINT2D::y
double y
Definition: liblwgeom.h:331
POINT2D::x
double x
Definition: liblwgeom.h:331
POINT4D::m
double m
Definition: liblwgeom.h:355
POINT4D::x
double x
Definition: liblwgeom.h:355
POINT4D::z
double z
Definition: liblwgeom.h:355
POINT4D::y
double y
Definition: liblwgeom.h:355
POINTARRAY::npoints
uint32_t npoints
Definition: liblwgeom.h:374
struct_lwgeom_parser_result
Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM.
Definition: liblwgeom.h:2013
uint32_t
unsigned int uint32_t
Definition: uthash.h:78