PostGIS  2.5.7dev-r@@SVN_REVISION@@
cu_measures.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * Copyright (C) 2015 Sandro Santilli <strk@kbt.io>
7  * Copyright (C) 2009 Paul Ramsey <pramsey@cleverelephant.ca>
8  *
9  * This is free software; you can redistribute and/or modify it under
10  * the terms of the GNU General Public Licence. See the COPYING file.
11  *
12  **********************************************************************/
13 
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <string.h>
17 #include "CUnit/Basic.h"
18 
19 #include "liblwgeom_internal.h"
20 #include "cu_tester.h"
21 #include "measures.h"
22 #include "measures3d.h"
23 #include "lwtree.h"
24 
25 static LWGEOM* lwgeom_from_text(const char *str)
26 {
27  LWGEOM_PARSER_RESULT r;
28  if( LW_FAILURE == lwgeom_parse_wkt(&r, (char*)str, LW_PARSER_CHECK_NONE) )
29  return NULL;
30  return r.geom;
31 }
32 
33 #define DIST2DTEST(str1, str2, res) \
34  do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance2d_tolerance);\
35  do_test_mindistance_tolerance(str2, str1, res, __LINE__, lwgeom_mindistance2d_tolerance)
36 #define DIST3DTEST(str1, str2, res) \
37  do_test_mindistance_tolerance(str1, str2, res, __LINE__, lwgeom_mindistance3d_tolerance)
38 
39 static void
40 do_test_mindistance_tolerance(char *in1,
41  char *in2,
42  double expected_res,
43  int line,
44  double (*distancef)(const LWGEOM *, const LWGEOM *, double))
45 {
46  LWGEOM *lw1;
47  LWGEOM *lw2;
48  double distance;
49  char *msg1 = "test_mindistance2d_tolerance failed (got %g expected %g) at line %d\n";
50  char *msg2 = "\n\ndo_test_mindistance2d_tolerance: NULL lwgeom generated from WKT\n %s\n\n";
51 
52  lw1 = lwgeom_from_wkt(in1, LW_PARSER_CHECK_NONE);
53  lw2 = lwgeom_from_wkt(in2, LW_PARSER_CHECK_NONE);
54 
55  if ( ! lw1 )
56  {
57  printf(msg2, in1);
58  exit(1);
59  }
60  if ( ! lw2 )
61  {
62  printf(msg2, in2);
63  exit(1);
64  }
65 
66  distance = distancef(lw1, lw2, 0.0);
67  lwgeom_free(lw1);
68  lwgeom_free(lw2);
69 
70  if ( fabs(distance - expected_res) > 0.00001 )
71  {
72  printf(msg1, distance, expected_res, line);
73  CU_FAIL();
74  }
75  else
76  {
77  CU_PASS();
78  }
79 
80 }
81 
82 static void test_mindistance2d_tolerance(void)
83 {
84  /*
85  ** Simple case.
86  */
87  DIST2DTEST("POINT(0 0)", "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5);
88 
89  /*
90  ** Point vs Geometry Collection.
91  */
92  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
93 
94  /*
95  ** Point vs Geometry Collection Collection.
96  */
97  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
98 
99  /*
100  ** Point vs Geometry Collection Collection Collection.
101  */
102  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0);
103 
104  /*
105  ** Point vs Geometry Collection Collection Collection Multipoint.
106  */
107  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0);
108 
109  /*
110  ** Geometry Collection vs Geometry Collection
111  */
112  DIST2DTEST("GEOMETRYCOLLECTION(POINT(0 0))", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
113 
114  /*
115  ** Geometry Collection Collection vs Geometry Collection Collection
116  */
117  DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
118 
119  /*
120  ** Geometry Collection Collection Multipoint vs Geometry Collection Collection Multipoint
121  */
122  DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4)))", 5.0);
123 
124  /*
125  ** Linestring vs its start point
126  */
127  DIST2DTEST("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0);
128 
129  /*
130  ** Linestring vs its end point
131  */
132  DIST2DTEST("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0);
133 
134  /*
135  ** Linestring vs its start point (tricky number, see #1459)
136  */
137  DIST2DTEST("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0);
138 
139  /*
140  ** Linestring vs its end point (tricky number, see #1459)
141  */
142  DIST2DTEST("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0);
143 
144  /*
145  * Circular string and point
146  */
147  DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "POINT(0 0)", 1);
148  DIST2DTEST("CIRCULARSTRING(-3 0, -2 0, -1 0, 0 1, 1 0)", "POINT(0 0)", 1);
149 
150  /*
151  * Circular string and Circular string
152  */
153  DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "CIRCULARSTRING(0 0, 1 -1, 2 0)", 1);
154 
155  /*
156  * CurvePolygon and Point
157  */
158  static char *cs1 = "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))";
159  DIST2DTEST(cs1, "POINT(3 14)", 1);
160  DIST2DTEST(cs1, "POINT(3 8)", 0);
161  DIST2DTEST(cs1, "POINT(6 5)", 1);
162  DIST2DTEST(cs1, "POINT(6 4)", 0);
163 
164  /*
165  * CurvePolygon and Linestring
166  */
167  DIST2DTEST(cs1, "LINESTRING(0 0, 50 0)", 0.917484);
168  DIST2DTEST(cs1, "LINESTRING(6 0, 10 7)", 0);
169  DIST2DTEST(cs1, "LINESTRING(4 4, 4 8)", 0);
170  DIST2DTEST(cs1, "LINESTRING(4 7, 5 6, 6 7)", 0.585786);
171  DIST2DTEST(cs1, "LINESTRING(10 0, 10 2, 10 0)", 1.52913);
172 
173  /*
174  * CurvePolygon and Polygon
175  */
176  DIST2DTEST(cs1, "POLYGON((10 4, 10 8, 13 8, 13 4, 10 4))", 0.58415);
177  DIST2DTEST(cs1, "POLYGON((9 4, 9 8, 12 8, 12 4, 9 4))", 0);
178  DIST2DTEST(cs1, "POLYGON((1 4, 1 8, 4 8, 4 4, 1 4))", 0);
179 
180  /*
181  * CurvePolygon and CurvePolygon
182  */
183  DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(-1 4, 0 5, 1 4, 0 3, -1 4))", 0.0475666);
184  DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(1 4, 2 5, 3 4, 2 3, 1 4))", 0.0);
185 
186  /*
187  * MultiSurface and CurvePolygon
188  */
189  static char *cs2 = "MULTISURFACE(POLYGON((0 0,0 4,4 4,4 0,0 0)),CURVEPOLYGON(CIRCULARSTRING(8 2,10 4,12 2,10 0,8 2)))";
190  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))", 1);
191  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))", 0);
192  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))", 0);
193  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))", 0);
194  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))", 0.5);
195  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))", 0);
196  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))", 0);
197  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))", 2.60555);
198 
199  /*
200  * MultiCurve and Linestring
201  */
202  DIST2DTEST("LINESTRING(0.5 1,0.5 3)", "MULTICURVE(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3),(0 0, 0 5))", 0.5);
203 
204  /*
205  ** Ticket 4326
206  */
207  DIST2DTEST(
208  "CURVEPOLYGON(CIRCULARSTRING(7874821 8715927,8907663 8715927,8844683 7750316,7937800 7750316,7874821 8715927))",
209  "POINT(5433865 8243495)", 2271704.2698450615);
210 }
211 
212 static void
213 test_mindistance3d_tolerance(void)
214 {
215  /* 2D [Z=0] should work just the same */
216  DIST3DTEST("POINT(0 0 0)", "MULTIPOINT(0 1.5 0, 0 2 0, 0 2.5 0)", 1.5);
217  DIST3DTEST("POINT(0 0 0)", "MULTIPOINT(0 1.5 0, 0 2 0, 0 2.5 0)", 1.5);
218  DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(POINT(3 4 0))", 5.0);
219  DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0)))", 5.0);
220  DIST3DTEST("POINT(0 0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0))))", 5.0);
221  DIST3DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0);
222  DIST3DTEST("GEOMETRYCOLLECTION(POINT(0 0 0))", "GEOMETRYCOLLECTION(POINT(3 4 0))", 5.0);
223  DIST3DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0 0)))",
224  "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4 0)))",
225  5.0);
226  DIST3DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0 0)))",
227  "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4 0)))",
228  5.0);
229  DIST3DTEST("LINESTRING(-2 0 0, -0.2 0 0)", "POINT(-2 0 0)", 0);
230  DIST3DTEST("LINESTRING(-0.2 0 0, -2 0 0)", "POINT(-2 0 0)", 0);
231  DIST3DTEST("LINESTRING(-1e-8 0 0, -0.2 0 0)", "POINT(-1e-8 0 0)", 0);
232  DIST3DTEST("LINESTRING(-0.2 0 0, -1e-8 0 0)", "POINT(-1e-8 0 0)", 0);
233 
234  /* Tests around intersections */
235  DIST3DTEST("LINESTRING(1 0 0 , 2 0 0)", "POLYGON((1 1 0, 2 1 0, 2 2 0, 1 2 0, 1 1 0))", 1.0);
236  DIST3DTEST("LINESTRING(1 1 1 , 2 1 0)", "POLYGON((1 1 0, 2 1 0, 2 2 0, 1 2 0, 1 1 0))", 0.0);
237  DIST3DTEST("LINESTRING(1 1 1 , 2 1 1)", "POLYGON((1 1 0, 2 1 0, 2 2 0, 1 2 0, 1 1 0))", 1.0);
238  /* Line in polygon */
239  DIST3DTEST("LINESTRING(1 1 1 , 2 2 2)", "POLYGON((0 0 0, 2 2 2, 3 3 1, 0 0 0))", 0.0);
240 
241  /* Line has a point in the same plane as the polygon but isn't the closest*/
242  DIST3DTEST("LINESTRING(-10000 -10000 0, 0 0 1)", "POLYGON((0 0 0, 1 0 0, 1 1 0, 0 1 0, 0 0 0))", 1);
243 
244  /* This is an invalid polygon since it defines just a line */
245  DIST3DTEST("LINESTRING(1 1 1 , 2 2 2)", "POLYGON((0 0 0, 2 2 2, 3 3 3, 0 0 0))", 0);
246 }
247 
248 static int tree_pt(RECT_NODE *tree, double x, double y)
249 {
250  POINT2D pt;
251  pt.x = x; pt.y = y;
252  return rect_tree_contains_point(tree, &pt);
253 }
254 
255 static void test_rect_tree_contains_point(void)
256 {
257  LWGEOM *poly;
258  RECT_NODE* tree;
259 
260  /**********************************************************************
261  * curvepolygon
262  */
263  poly = lwgeom_from_wkt("CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0,1 5,0 10),(0 10,10 10,10 0, 0 0)),COMPOUNDCURVE(CIRCULARSTRING(3 7,5 8,7 7),(7 7,7 3,3 3, 3 7)))", LW_PARSER_CHECK_NONE);
264  tree = rect_tree_from_lwgeom(poly);
265  // char *wkt = rect_tree_to_wkt(tree);
266  // printf("%s\n", wkt);
267  // lwfree(wkt);
268  // return;
269 
270  /* in hole, within arc stroke */
271  CU_ASSERT_EQUAL(tree_pt(tree, 5, 7.5), 0);
272  /* inside */
273  CU_ASSERT_EQUAL(tree_pt(tree, 8, 9), 1);
274  /* outside */
275  CU_ASSERT_EQUAL(tree_pt(tree, -1, 5), 0);
276  /* outside */
277  CU_ASSERT_EQUAL(tree_pt(tree, -1, 7.5), 0);
278  /* outside, within arc stroke */
279  CU_ASSERT_EQUAL(tree_pt(tree, 0.2, 7.5), 0);
280  /* inside, within arc stroke */
281  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 0.5), 1);
282  /* inside, crossing arc stroke */
283  CU_ASSERT_EQUAL(tree_pt(tree, 2, 7.5), 1);
284  /* touching hole corner */
285  CU_ASSERT_EQUAL(tree_pt(tree, 7, 7), 1);
286 
287  rect_tree_free(tree);
288  lwgeom_free(poly);
289 
290 
291  /**********************************************************************
292  * polygon with hole and concavities
293  */
294  poly = lwgeom_from_wkt("POLYGON((0 0,0 10,10 10,10 0,9 0,9 9,8 6,8 0,2 0,2 9,1 6,1 0,0 0),(4 4,4 6,6 6,6 4,4 4))", LW_PARSER_CHECK_NONE);
295  tree = rect_tree_from_lwgeom(poly);
296 
297  /* inside, many grazings */
298  CU_ASSERT_EQUAL(tree_pt(tree, 3, 6), 1);
299  /* inside */
300  CU_ASSERT_EQUAL(tree_pt(tree, 3, 5.5), 1);
301  /* outside */
302  CU_ASSERT_EQUAL(tree_pt(tree, -3, 5.5), 0);
303  /* touching interior ring */
304  CU_ASSERT_EQUAL(tree_pt(tree, 4, 4), 1);
305  CU_ASSERT_EQUAL(tree_pt(tree, 6, 6), 1);
306  /* touching interior ring */
307  CU_ASSERT_EQUAL(tree_pt(tree, 4.5, 4), 1);
308  /* touching exterior ring */
309  CU_ASSERT_EQUAL(tree_pt(tree, 8, 0), 1);
310  CU_ASSERT_EQUAL(tree_pt(tree, 9, 0), 1);
311  CU_ASSERT_EQUAL(tree_pt(tree, 10, 1), 1);
312  CU_ASSERT_EQUAL(tree_pt(tree, 9.5, 1), 1);
313  CU_ASSERT_EQUAL(tree_pt(tree, 0, 10), 1);
314  /* touching grazing spike */
315  CU_ASSERT_EQUAL(tree_pt(tree, 1, 6), 1);
316  /* outide, many grazings */
317  CU_ASSERT_EQUAL(tree_pt(tree, -1, 6), 0);
318  /* within hole */
319  CU_ASSERT_EQUAL(tree_pt(tree, 5, 5), 0);
320  /* within */
321  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 4), 1);
322  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 6), 1);
323  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 9), 1);
324 
325  rect_tree_free(tree);
326  lwgeom_free(poly);
327 
328  /**********************************************************************
329  * square
330  */
331  poly = lwgeom_from_wkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", LW_PARSER_CHECK_NONE);
332  tree = rect_tree_from_lwgeom(poly);
333 
334  /* inside square */
335  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 0.5), 1);
336  /* outside square */
337  CU_ASSERT_EQUAL(tree_pt(tree, 1.5, 0.5), 0);
338  /* outside square grazing some edges */
339  CU_ASSERT_EQUAL(tree_pt(tree, -1, 1), 0);
340  /* inside square on corner */
341  CU_ASSERT_EQUAL(tree_pt(tree, 1, 1), 1);
342  /* inside square on top edge */
343  CU_ASSERT_EQUAL(tree_pt(tree, 0.5, 1), 1);
344  /* inside square on side edge */
345  CU_ASSERT_EQUAL(tree_pt(tree, 1, 0.5), 1);
346 
347  rect_tree_free(tree);
348  lwgeom_free(poly);
349 
350  /* ziggy zaggy horizontal saw tooth polygon */
351  poly = lwgeom_from_wkt("POLYGON((0 0, 1 3, 2 0, 3 3, 4 0, 4 5, 0 5, 0 0))", LW_PARSER_CHECK_NONE);
352  tree = rect_tree_from_lwgeom(poly);
353 
354  /* not in, left side */
355  CU_ASSERT_EQUAL(tree_pt(tree, -0.5, 0.5), 0);
356  /* not in, right side */
357  CU_ASSERT_EQUAL(tree_pt(tree, 3, 1), 0);
358  /* inside */
359  CU_ASSERT_EQUAL(tree_pt(tree, 2, 1), 1);
360  /* on left border */
361  CU_ASSERT_EQUAL(tree_pt(tree, 0, 1), 1);
362  /* on left border, grazing */
363  CU_ASSERT_EQUAL(tree_pt(tree, 0, 3), 1);
364  /* on right border */
365  CU_ASSERT_EQUAL(tree_pt(tree, 4, 0), 1);
366  /* on tooth concave */
367  CU_ASSERT_EQUAL(tree_pt(tree, 3, 3), 1);
368  /* on tooth convex */
369  CU_ASSERT_EQUAL(tree_pt(tree, 2, 0), 1);
370 
371  rect_tree_free(tree);
372  lwgeom_free(poly);
373 
374  /**********************************************************************
375  * ziggy zaggy vertical saw tooth polygon
376  */
377  poly = lwgeom_from_wkt("POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))", LW_PARSER_CHECK_NONE);
378  tree = rect_tree_from_lwgeom(poly);
379 
380  /* not in, left side */
381  CU_ASSERT_EQUAL(tree_pt(tree, -0.5, 3.5), 0);
382  /* not in, right side */
383  CU_ASSERT_EQUAL(tree_pt(tree, 6.0, 2.2), 0);
384  /* inside */
385  CU_ASSERT_EQUAL(tree_pt(tree, 3, 2), 1);
386  /* on bottom border */
387  CU_ASSERT_EQUAL(tree_pt(tree, 1, 0), 1);
388  /* on top border */
389  CU_ASSERT_EQUAL(tree_pt(tree, 3, 6), 1);
390  /* on tooth concave */
391  CU_ASSERT_EQUAL(tree_pt(tree, 3, 1), 1);
392  /* on tooth convex */
393  CU_ASSERT_EQUAL(tree_pt(tree, 0, 2), 1);
394  /* on tooth convex */
395  CU_ASSERT_EQUAL(tree_pt(tree, 0, 6), 1);
396 
397  rect_tree_free(tree);
398  lwgeom_free(poly);
399 
400 }
401 
402 static int tree_inter(const char *wkt1, const char *wkt2)
403 {
404  LWGEOM *g1 = lwgeom_from_wkt(wkt1, LW_PARSER_CHECK_NONE);
405  LWGEOM *g2 = lwgeom_from_wkt(wkt2, LW_PARSER_CHECK_NONE);
406  RECT_NODE *t1 = rect_tree_from_lwgeom(g1);
407  RECT_NODE *t2 = rect_tree_from_lwgeom(g2);
408  int result = rect_tree_intersects_tree(t1, t2);
409  rect_tree_free(t1);
410  rect_tree_free(t2);
411  lwgeom_free(g1);
412  lwgeom_free(g2);
413  return result;
414 }
415 
416 static void test_rect_tree_intersects_tree(void)
417 {
418  /* total overlap, A == B */
419  CU_ASSERT_EQUAL(tree_inter(
420  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
421  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))"),
422  LW_TRUE
423  );
424 
425  /* hiding between the tines of the comb */
426  CU_ASSERT_EQUAL(tree_inter(
427  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
428  "POLYGON((0.3 0.7, 0.3 0.8, 0.4 0.8, 0.4 0.7, 0.3 0.7))"),
429  LW_FALSE
430  );
431 
432  /* between the tines, but with a corner overlapping */
433  CU_ASSERT_EQUAL(tree_inter(
434  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
435  "POLYGON((0.3 0.7, 0.3 0.8, 0.4 0.8, 1.3 0.3, 0.3 0.7))"),
436  LW_TRUE
437  );
438 
439  /* Just touching the top left corner of the comb */
440  CU_ASSERT_EQUAL(tree_inter(
441  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
442  "POLYGON((-1 5, 0 5, 0 7, -1 7, -1 5))"),
443  LW_TRUE
444  );
445 
446  /* Contained, complex */
447  CU_ASSERT_EQUAL(tree_inter(
448  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
449  "GEOMETRYCOLLECTION(MULTILINESTRING((1 2, 3 2)),POINT(1 2))"),
450  LW_TRUE
451  );
452 
453  /* Touching, complex */
454  CU_ASSERT_EQUAL(tree_inter(
455  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
456  "GEOMETRYCOLLECTION(MULTILINESTRING((6 3, 8 4)),POINT(5 3))"),
457  LW_TRUE
458  );
459 
460  /* Not Touching, complex */
461  CU_ASSERT_EQUAL(tree_inter(
462  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
463  "GEOMETRYCOLLECTION(MULTILINESTRING((6 3, 8 4)),POINT(1 3.5))"),
464  LW_FALSE
465  );
466 
467  /* Crossing, complex */
468  CU_ASSERT_EQUAL(tree_inter(
469  "POLYGON((0 0, 3 1, 0 2, 3 3, 0 4, 3 5, 0 6, 5 6, 5 0, 0 0))",
470  "GEOMETRYCOLLECTION(MULTILINESTRING((1.5 4.1, 1.6 2)),POINT(1 3.5))"),
471  LW_TRUE
472  );
473 }
474 
475 
476 static double
477 test_rect_tree_distance_tree_case(const char *wkt1, const char *wkt2)
478 {
479  LWGEOM *lw1 = lwgeom_from_wkt(wkt1, LW_PARSER_CHECK_NONE);
480  LWGEOM *lw2 = lwgeom_from_wkt(wkt2, LW_PARSER_CHECK_NONE);
481  RECT_NODE *n1 = rect_tree_from_lwgeom(lw1);
482  RECT_NODE *n2 = rect_tree_from_lwgeom(lw2);
483 
484  // rect_tree_printf(n1, 0);
485  // rect_tree_printf(n2, 0);
486  //
487  // printf("%s\n", rect_tree_to_wkt(n1));
488  // printf("%s\n", rect_tree_to_wkt(n2));
489 
490  double dist = rect_tree_distance_tree(n1, n2, 0.0);
491  // printf("%g\n", dist);
492  rect_tree_free(n1);
493  rect_tree_free(n2);
494  lwgeom_free(lw1);
495  lwgeom_free(lw2);
496  return dist;
497 }
498 
499 #define TDT(w1, w2, d) CU_ASSERT_DOUBLE_EQUAL(test_rect_tree_distance_tree_case(w1, w2), d, 0.00001);
500 
501 static void
502 test_rect_tree_distance_tree(void)
503 {
504  const char *wkt;
505 
506  wkt = "MULTIPOLYGON(((-123.35702791281 48.4232302445918,-123.35689654493 48.4237265810249,-123.354053908057 48.4234039978588,-123.35417179975 48.4229151379279,-123.354369811539 48.4220987102936,-123.355779071731 48.4222571534228,-123.357238860904 48.4224209369449,-123.35702791281 48.4232302445918)))";
507  TDT(wkt, "MULTIPOLYGON(((-123.353452578038 48.4259519079838,-123.35072012771 48.4256699150083,-123.347337809991 48.4254740864963,-123.347469111645 48.4245757659326,-123.349409235923 48.4246224093429,-123.349966167324 48.4246562342604,-123.353650661317 48.4250703224683,-123.353452578038 48.4259519079838)))", 0.0017144228293396);
508 
509  wkt = "POINT(0 0)";
510  TDT(wkt, "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5);
511  TDT(wkt, "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
512  TDT(wkt, "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
513  TDT(wkt, "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0);
514  TDT(wkt, "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0);
515 
516  TDT("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0);
517  TDT("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0);
518  TDT("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0);
519  TDT("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0);
520 
521  wkt = "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))";
522  TDT(wkt, "POINT(3 14)", 1);
523  TDT(wkt, "POINT(3 8)", 0);
524  TDT(wkt, "POINT(6 5)", 1);
525  TDT(wkt, "POINT(6 4)", 0);
526 
527  wkt = "MULTISURFACE(POLYGON((0 0,0 4,4 4,4 0,0 0)),CURVEPOLYGON(CIRCULARSTRING(8 2,10 4,12 2,10 0,8 2)))";
528  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))", 2.60555);
529  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))", 1);
530  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))", 0);
531  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))", 0);
532  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))", 0);
533  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))", 0.5);
534  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))", 0);
535  TDT(wkt, "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))", 0);
536 
537  wkt = "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0,5 0,0 0)))";
538  TDT(wkt, "POINT(3 0)", 0.0);
539  TDT(wkt, "POINT(5 0)", 0.0);
540  TDT(wkt, "POINT(7 0)", 2.0);
541  TDT(wkt, "POINT(2.5 3.5)", 1.0);
542 
543  wkt = "POINT(0 0)";
544  TDT(wkt, "POINT(0 1)", 1.0);
545  TDT(wkt, "POINT(1 0)", 1.0);
546 
547  wkt = "LINESTRING(0 0,1 0)";
548  TDT(wkt, "LINESTRING(1 0,1 1)", 0.0);
549  TDT(wkt, "LINESTRING(0 1,1 1)", 1.0);
550 
551  wkt = "POLYGON((0 0,0 1,1 1,1 0,0 0))";
552  TDT(wkt, "POINT(2 2)", sqrt(2));
553  TDT(wkt, "POINT(0.5 0.5)", 0);
554  TDT(wkt, "POINT(1 1)", 0);
555 
556  wkt = "POLYGON((0 0,0 10,10 10,10 0,0 0), (4 4,4 6,6 6,6 4,4 4))";
557  TDT(wkt, "POINT(5 5)", 1);
558  TDT(wkt, "POLYGON((5 5,5 5.5,5.5 5.5,5.5 5, 5 5))", 0.5);
559 }
560 
561 
562 static void
563 test_lwgeom_segmentize2d(void)
564 {
565  LWGEOM *linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
566  LWGEOM *lineout = lwgeom_segmentize2d(linein, 5);
567  char *strout = lwgeom_to_ewkt(lineout);
568  ASSERT_STRING_EQUAL(strout, "LINESTRING(0 0,5 0,10 0)");
569  lwfree(strout);
570  lwgeom_free(linein);
571  lwgeom_free(lineout);
572 
573  /* test interruption */
574 
575  linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
576  lwgeom_request_interrupt();
577  lineout = lwgeom_segmentize2d(linein, 1e-100);
578  CU_ASSERT_EQUAL(lineout, NULL);
579  lwgeom_free(linein);
580 
581  linein = lwgeom_from_wkt("MULTILINESTRING((0 0,10 0),(20 0, 30 0))", LW_PARSER_CHECK_NONE);
582  lwgeom_request_interrupt();
583  lineout = lwgeom_segmentize2d(linein, 1e-100);
584  CU_ASSERT_EQUAL(lineout, NULL);
585  lwgeom_free(linein);
586 
587  linein = lwgeom_from_wkt(
588  "MULTIPOLYGON(((0 0,20 0,20 20,0 20,0 0),(2 2,2 4,4 4,4 2,2 2),(6 6,6 8,8 8,8 6,6 6)),((40 0,40 20,60 20,60 0,40 0),(42 2,42 4,44 4,44 2,42 2)))"
589  , LW_PARSER_CHECK_NONE);
590  lwgeom_request_interrupt();
591  lineout = lwgeom_segmentize2d(linein, 1e-100);
592  CU_ASSERT_EQUAL(lineout, NULL);
593  lwgeom_free(linein);
594 
595  linein = lwgeom_from_wkt(
596  "GEOMETRYCOLLECTION(MULTIPOLYGON(((0 0,20 0,20 20,0 20,0 0),(2 2,2 4,4 4,4 2,2 2),(6 6,6 8,8 8,8 6,6 6)),((40 0,40 20,60 20,60 0,40 0),(42 2,42 4,44 4,44 2,42 2))),MULTILINESTRING((0 0,10 0),(20 0, 30 0)),MULTIPOINT(0 0, 3 4))"
597  , LW_PARSER_CHECK_NONE);
598  CU_ASSERT_FATAL(linein != NULL);
599  lwgeom_request_interrupt();
600  lineout = lwgeom_segmentize2d(linein, 1e-100);
601  CU_ASSERT_EQUAL(lineout, NULL);
602  lwgeom_free(linein);
603 
604  linein = lwgeom_from_wkt("LINESTRING(20 0, 30 0)", LW_PARSER_CHECK_NONE);
605  /* NOT INTERRUPTED */
606  lineout = lwgeom_segmentize2d(linein, 5);
607  strout = lwgeom_to_ewkt(lineout);
608  ASSERT_STRING_EQUAL(strout, "LINESTRING(20 0,25 0,30 0)");
609  lwfree(strout);
610  lwgeom_free(linein);
611  lwgeom_free(lineout);
612 }
613 
614 static void
615 test_lwgeom_locate_along(void)
616 {
617  LWGEOM *geom = NULL;
618  LWGEOM *out = NULL;
619  double measure = 105.0;
620  char *str;
621 
622  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
623  geom = lwgeom_from_wkt("MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))", LW_PARSER_CHECK_NONE);
624  out = lwgeom_locate_along(geom, measure, 0.0);
625  str = lwgeom_to_wkt(out, WKT_ISO, 6, NULL);
626  lwgeom_free(geom);
627  lwgeom_free(out);
628  ASSERT_STRING_EQUAL(str, "MULTIPOINT M (55.226131 55.226131 105)");
629  lwfree(str);
630 
631  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
632  geom = lwgeom_from_wkt("MULTILINESTRING M ((1 2 3, 3 4 2, 9 4 3), (1 2 3, 5 4 5), (50 50 1, 60 60 200))", LW_PARSER_CHECK_NONE);
633  out = lwgeom_locate_along(geom, measure, 0.0);
634  str = lwgeom_to_wkt(out, WKT_ISO, 6, NULL);
635  lwgeom_free(geom);
636  lwgeom_free(out);
637  ASSERT_STRING_EQUAL(str, "MULTIPOINT M (55.226131 55.226131 105)");
638  lwfree(str);
639 }
640 
641 static void
642 test_lw_dist2d_pt_arc(void)
643 {
644  /* int lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const POINT2D* A3, DISTPTS* dl) */
645  DISTPTS dl;
646  POINT2D P, A1, A2, A3;
647  int rv;
648 
649 
650  /* Point within unit semicircle, 0.5 units from arc */
651  A1.x = -1; A1.y = 0;
652  A2.x = 0 ; A2.y = 1;
653  A3.x = 1 ; A3.y = 0;
654  P.x = 0 ; P.y = 0.5;
655 
656  lw_dist2d_distpts_init(&dl, DIST_MIN);
657  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
658  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
659  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
660 
661  /* Point outside unit semicircle, 0.5 units from arc */
662  P.x = 0 ; P.y = 1.5;
663  lw_dist2d_distpts_init(&dl, DIST_MIN);
664  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
665  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
666  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
667 
668  /* Point outside unit semicircle, sqrt(2) units from arc end point*/
669  P.x = 0 ; P.y = -1;
670  lw_dist2d_distpts_init(&dl, DIST_MIN);
671  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
672  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
673  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);
674 
675  /* Point outside unit semicircle, sqrt(2)-1 units from arc end point*/
676  P.x = 1 ; P.y = 1;
677  lw_dist2d_distpts_init(&dl, DIST_MIN);
678  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
679  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
680  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
681 
682  /* Point on unit semicircle midpoint */
683  P.x = 0 ; P.y = 1;
684  lw_dist2d_distpts_init(&dl, DIST_MIN);
685  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
686  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
687  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
688 
689  /* Point on unit semicircle endpoint */
690  P.x = 1 ; P.y = 0;
691  lw_dist2d_distpts_init(&dl, DIST_MIN);
692  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
693  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
694  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
695 
696  /* Point on semicircle center */
697  P.x = 0 ; P.y = 0;
698  lw_dist2d_distpts_init(&dl, DIST_MIN);
699  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
700  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
701  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
702 
703  /* Point inside closed circle */
704  P.x = 0 ; P.y = 0.5;
705  A2.x = 1; A2.y = 0;
706  A3 = A1;
707  lw_dist2d_distpts_init(&dl, DIST_MIN);
708  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
709  //printf("distance %g\n", dl.distance);
710  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
711  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
712 }
713 
714 static void
715 test_lw_dist2d_seg_arc(void)
716 {
717  /* int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl) */
718 
719  DISTPTS dl;
720  POINT2D A1, A2, B1, B2, B3;
721  int rv;
722 
723  /* Unit semicircle */
724  B1.x = -1; B1.y = 0;
725  B2.x = 0 ; B2.y = 1;
726  B3.x = 1 ; B3.y = 0;
727 
728  /* Edge above the unit semicircle */
729  lw_dist2d_distpts_init(&dl, DIST_MIN);
730  A1.x = -2; A1.y = 2;
731  A2.x = 2 ; A2.y = 2;
732  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
733  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
734  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
735 
736  /* Edge to the right of the unit semicircle */
737  lw_dist2d_distpts_init(&dl, DIST_MIN);
738  A1.x = 2; A1.y = -2;
739  A2.x = 2; A2.y = 2;
740  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
741  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
742  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
743 
744  /* Edge to the left of the unit semicircle */
745  lw_dist2d_distpts_init(&dl, DIST_MIN);
746  A1.x = -2; A1.y = -2;
747  A2.x = -2; A2.y = 2;
748  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
749  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
750  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
751 
752  /* Edge within the unit semicircle */
753  lw_dist2d_distpts_init(&dl, DIST_MIN);
754  A1.x = 0; A1.y = 0;
755  A2.x = 0; A2.y = 0.5;
756  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
757  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
758  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
759 
760  /* Edge grazing the unit semicircle */
761  lw_dist2d_distpts_init(&dl, DIST_MIN);
762  A1.x = -2; A1.y = 1;
763  A2.x = 2; A2.y = 1;
764  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
765  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
766  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0., 0.000001);
767 
768  /* Line grazing the unit semicircle, but edge not */
769  lw_dist2d_distpts_init(&dl, DIST_MIN);
770  A1.x = 1; A1.y = 1;
771  A2.x = 2; A2.y = 1;
772  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
773  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
774  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
775 
776  /* Edge intersecting the unit semicircle */
777  lw_dist2d_distpts_init(&dl, DIST_MIN);
778  A1.x = 0; A1.y = 0;
779  A2.x = 2; A2.y = 2;
780  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
781  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
782  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
783 
784  /* Line intersecting the unit semicircle, but edge not */
785  lw_dist2d_distpts_init(&dl, DIST_MIN);
786  A1.x = -1; A1.y = 1;
787  A2.x = -2; A2.y = 2;
788  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
789  //printf("distance %g\n", dl.distance);
790  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
791  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
792 }
793 
794 static void
795 test_lw_dist2d_arc_arc(void)
796 {
797  /* lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
798  const POINT2D *B1, const POINT2D *B2, const POINT2D *B3,
799  DISTPTS *dl) */
800  DISTPTS dl;
801  POINT2D A1, A2, A3, B1, B2, B3;
802  int rv;
803 
804  /* Ticket #4326 */
805  lw_dist2d_distpts_init(&dl, DIST_MIN);
806  A1.x = -1.0; A1.y = 4.0;
807  A2.x = 0.0; A2.y = 5.0;
808  A3.x = 1.0; A3.y = 4.0;
809  B1.x = 1.0; B1.y = 6.0;
810  B2.x = 6.0; B2.y = 1.0;
811  B3.x = 9.0; B3.y = 7.0;
812  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
813  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
814  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.0475666, 0.000001);
815 
816  /* Unit semicircle at 0,0 */
817  B1.x = -1; B1.y = 0;
818  B2.x = 0 ; B2.y = 1;
819  B3.x = 1 ; B3.y = 0;
820 
821  /* Arc above the unit semicircle */
822  lw_dist2d_distpts_init(&dl, DIST_MIN);
823  A1.x = -1; A1.y = 3;
824  A2.x = 0 ; A2.y = 2;
825  A3.x = 1 ; A3.y = 3;
826  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
827  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
828  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
829 
830  /* Arc grazes the unit semicircle */
831  lw_dist2d_distpts_init(&dl, DIST_MIN);
832  A1.x = -1; A1.y = 2;
833  A2.x = 0 ; A2.y = 1;
834  A3.x = 1 ; A3.y = 2;
835  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
836  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
837  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
838 
839  /* Circles intersect, but arcs do not */
840  lw_dist2d_distpts_init(&dl, DIST_MIN);
841  A1.x = -1; A1.y = 1;
842  A2.x = 0; A2.y = 2;
843  A3.x = 1; A3.y = 1;
844  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
845  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
846  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2)-1, 0.000001);
847 
848  /* Circles and arcs intersect */
849  lw_dist2d_distpts_init(&dl, DIST_MIN);
850  A1.x = -1; A1.y = 1;
851  A2.x = 0; A2.y = 0;
852  A3.x = 1; A3.y = 1;
853  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
854  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
855  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
856 
857  /* Concentric: and fully parallel */
858  lw_dist2d_distpts_init(&dl, DIST_MIN);
859  A1.x = -2.0; A1.y = 0.0;
860  A2.x = 0.0; A2.y = 2.0;
861  A3.x = 2.0; A3.y = 0.0;
862  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
863  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
864  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
865 
866  /* Concentric: with A fully included in B's range */
867  lw_dist2d_distpts_init(&dl, DIST_MIN);
868  A1.x = -0.5 / sqrt(2.0); A1.y = 0.5 / sqrt(2.0);
869  A2.x = 0.0; A2.y = 0.5;
870  A3.x = 0.5 / sqrt(2.0); A3.y = 0.5 / sqrt(2.0);
871  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
872  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
873  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
874 
875  /* Concentric: with A partially included in B's range */
876  lw_dist2d_distpts_init(&dl, DIST_MIN);
877  A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
878  A2.x = -0.5; A2.y = 0.0;
879  A3.x = -0.5 / sqrt(2.0); A3.y = 0.5 / sqrt(2.0);
880  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
881  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
882  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
883 
884  /* Concentric: with A and B without parallel segments */
885  lw_dist2d_distpts_init(&dl, DIST_MIN);
886  A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
887  A2.x = 0.0; A2.y = -0.5;
888  A3.x = 0.5 / sqrt(2.0); A3.y = -0.5 / sqrt(2.0);
889  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
890  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
891  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.736813, 0.000001);
892 
893  /* Concentric: Arcs are the same */
894  lw_dist2d_distpts_init(&dl, DIST_MIN);
895  A1.x = -1.0; A1.y = 0.0;
896  A2.x = 0.0; A2.y = 1.0;
897  A3.x = 1.0; A3.y = 0.0;
898  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
899  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
900  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.0, 0.000001);
901 
902  /* Check different orientations */
903  B1.x = -10.0; B1.y = 0.0;
904  B2.x = 0.0 ; B2.y = 10.0;
905  B3.x = 10.0 ; B3.y = 0.0;
906 
907  lw_dist2d_distpts_init(&dl, DIST_MIN);
908  A1.x = -22.0; A1.y = 0.0;
909  A2.x = -17.0; A2.y = -5.0;
910  A3.x = -12.0; A3.y = 0.0;
911  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
912  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
913  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 2.0, 0.000001);
914 
915  lw_dist2d_distpts_init(&dl, DIST_MIN);
916  A1.x = -19.0; A1.y = 0.0;
917  A2.x = -14.0; A2.y = -5.0;
918  A3.x = - 9.0; A3.y = 0.0;
919  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
920  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
921  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
922 
923  lw_dist2d_distpts_init(&dl, DIST_MIN);
924  A1.x = -9.0; A1.y = 0.0;
925  A2.x = -4.0; A2.y = -5.0;
926  A3.x = 1.0; A3.y = 0.0;
927  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
928  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
929  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
930 
931  lw_dist2d_distpts_init(&dl, DIST_MIN);
932  A1.x = -1.0; A1.y = 0.0;
933  A2.x = 4.0; A2.y = -5.0;
934  A3.x = 9.0; A3.y = 0.0;
935  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
936  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
937  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
938 
939  lw_dist2d_distpts_init(&dl, DIST_MIN);
940  A1.x = 1.0; A1.y = 0.0;
941  A2.x = 6.0; A2.y = -5.0;
942  A3.x = 11.0; A3.y = 0.0;
943  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
944  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
945  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
946 
947  lw_dist2d_distpts_init(&dl, DIST_MIN);
948  A1.x = 11.0; A1.y = 0.0;
949  A2.x = 16.0; A2.y = -5.0;
950  A3.x = 21.0; A3.y = 0.0;
951  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
952  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
953  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
954 
955 
956  lw_dist2d_distpts_init(&dl, DIST_MIN);
957  A1.x = -15.0; A1.y = -6.0;
958  A2.x = -10.0; A2.y = -1.0;
959  A3.x = - 5.0; A3.y = -6.0;
960  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
961  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
962  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
963 
964  lw_dist2d_distpts_init(&dl, DIST_MIN);
965  A1.x = -5.0; A1.y = 0.0;
966  A2.x = 0.0; A2.y = 5.0;
967  A3.x = 5.0; A3.y = 0.0;
968  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
969  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
970  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 5.0, 0.000001);
971 
972  lw_dist2d_distpts_init(&dl, DIST_MIN);
973  A1.x = -5.0; A1.y = 0.0;
974  A2.x = 0.0; A2.y = -5.0;
975  A3.x = 5.0; A3.y = 0.0;
976  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
977  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
978  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 5.0, 0.000001);
979 
980 }
981 
982 static void
983 test_lw_arc_length(void)
984 {
985 /* double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3) */
986 
987  POINT2D A1, A2, A3;
988  double d;
989 
990  /* Unit semicircle at 0,0 */
991  A1.x = -1; A1.y = 0;
992  A2.x = 0 ; A2.y = 1;
993  A3.x = 1 ; A3.y = 0;
994 
995  /* Arc above the unit semicircle */
996  d = lw_arc_length(&A1, &A2, &A3);
997  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
998  d = lw_arc_length(&A3, &A2, &A1);
999  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1000 
1001  /* Unit semicircle at 0,0 */
1002  A1.x = 0; A1.y = 1;
1003  A2.x = 1; A2.y = 0;
1004  A3.x = 0; A3.y = -1;
1005 
1006  /* Arc to right of the unit semicircle */
1007  d = lw_arc_length(&A1, &A2, &A3);
1008  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1009  d = lw_arc_length(&A3, &A2, &A1);
1010  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
1011 
1012  /* Unit 3/4 circle at 0,0 */
1013  A1.x = -1; A1.y = 0;
1014  A2.x = 1; A2.y = 0;
1015  A3.x = 0; A3.y = -1;
1016 
1017  /* Arc to right of the unit semicircle */
1018  d = lw_arc_length(&A1, &A2, &A3);
1019  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
1020  d = lw_arc_length(&A3, &A2, &A1);
1021  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
1022 }
1023 
1024 static void
1025 test_lw_dist2d_pt_ptarrayarc(void)
1026 {
1027  /* lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl) */
1028  DISTPTS dl;
1029  int rv;
1030  LWLINE *lwline;
1031  POINT2D P;
1032 
1033  /* Unit semi-circle above X axis */
1034  lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
1035 
1036  /* Point at origin */
1037  P.x = P.y = 0;
1038  lw_dist2d_distpts_init(&dl, DIST_MIN);
1039  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1040  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1041  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1042 
1043  /* Point above arc on Y axis */
1044  P.y = 2;
1045  lw_dist2d_distpts_init(&dl, DIST_MIN);
1046  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1047  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1048  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1049 
1050  /* Point 45 degrees off arc, 2 radii from center */
1051  P.y = P.x = 2 * cos(M_PI_4);
1052  lw_dist2d_distpts_init(&dl, DIST_MIN);
1053  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1054  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1055  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1056 
1057  /* Four unit semi-circles surrounding the 2x2 box around origin */
1058  lwline_free(lwline);
1059  lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 -1, -2 0, -1 1, 0 2, 1 1, 2 0, 1 -1, 0 -2, -1 -1)"));
1060 
1061  /* Point at origin */
1062  P.x = P.y = 0;
1063  lw_dist2d_distpts_init(&dl, DIST_MIN);
1064  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1065  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1066  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);
1067 
1068  /* Point on box edge */
1069  P.x = -1; P.y = 0;
1070  lw_dist2d_distpts_init(&dl, DIST_MIN);
1071  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1072  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1073  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1074 
1075  /* Point within a semicircle lobe */
1076  P.x = -1.5; P.y = 0;
1077  lw_dist2d_distpts_init(&dl, DIST_MIN);
1078  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1079  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1080  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
1081 
1082  /* Point outside a semicircle lobe */
1083  P.x = -2.5; P.y = 0;
1084  lw_dist2d_distpts_init(&dl, DIST_MIN);
1085  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1086  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1087  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
1088 
1089  /* Point outside a semicircle lobe */
1090  P.y = -2.5; P.x = 0;
1091  lw_dist2d_distpts_init(&dl, DIST_MIN);
1092  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1093  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1094  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
1095 
1096  /* Point outside a semicircle lobe */
1097  P.y = 2; P.x = 1;
1098  lw_dist2d_distpts_init(&dl, DIST_MIN);
1099  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
1100  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1101  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1.0, 0.000001);
1102 
1103  /* Clean up */
1104  lwline_free(lwline);
1105 }
1106 
1107 static void
1108 test_lw_dist2d_ptarray_ptarrayarc(void)
1109 {
1110  /* int lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl) */
1111  DISTPTS dl;
1112  int rv;
1113  LWLINE *lwline1;
1114  LWLINE *lwline2;
1115 
1116  /* Unit semi-circle above X axis */
1117  lwline1 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
1118 
1119  /* Line above top of semi-circle */
1120  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 2, -1 2, 1 2, 2 2)"));
1121  lw_dist2d_distpts_init(&dl, DIST_MIN);
1122  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
1123  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1124  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1125 
1126  /* Reversed arguments, should fail */
1127  lw_dist2d_distpts_init(&dl, DIST_MIN);
1128  cu_error_msg_reset();
1129  rv = lw_dist2d_ptarray_ptarrayarc(lwline1->points, lwline2->points, &dl);
1130  //printf("%s\n", cu_error_msg);
1131  CU_ASSERT_EQUAL( rv, LW_FAILURE );
1132  ASSERT_STRING_EQUAL(
1133  cu_error_msg,
1134  "lw_dist2d_ptarray_ptarrayarc called with non-arc input");
1135 
1136  lwline_free(lwline2);
1137 
1138  /* Line along side of semi-circle */
1139  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 -3, -2 -2, -2 2, -2 3)"));
1140  lw_dist2d_distpts_init(&dl, DIST_MIN);
1141  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
1142  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1143  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
1144 
1145  /* Four unit semi-circles surrounding the 2x2 box around origin */
1146  lwline_free(lwline1);
1147  lwline_free(lwline2);
1148  lwline1 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 -1, -2 0, -1 1, 0 2, 1 1, 2 0, 1 -1, 0 -2, -1 -1)"));
1149  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2.5 -3, -2.5 -2, -2.5 2, -2.5 3)"));
1150  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
1151  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
1152  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
1153 
1154  lwline_free(lwline2);
1155  lwline_free(lwline1);
1156 }
1157 
1158 static void
1159 test_lwgeom_tcpa(void)
1160 {
1161  LWGEOM *g1, *g2;
1162  double m, dist;
1163 
1164  /* Invalid input, lack of dimensions */
1165 
1166  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1167  g2 = lwgeom_from_wkt("LINESTRING (0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
1168  m = lwgeom_tcpa(g1, g2, NULL);
1169  lwgeom_free(g1);
1170  lwgeom_free(g2);
1171  ASSERT_DOUBLE_EQUAL(m, -1.0);
1172  ASSERT_STRING_EQUAL(
1173  cu_error_msg,
1174  "Both input geometries must have a measure dimension");
1175 
1176  /* Invalid input, not linestrings */
1177 
1178  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1179  g2 = lwgeom_from_wkt("POINT M (0 0 2)", LW_PARSER_CHECK_NONE);
1180  m = lwgeom_tcpa(g1, g2, NULL);
1181  lwgeom_free(g1);
1182  lwgeom_free(g2);
1183  ASSERT_DOUBLE_EQUAL(m, -1.0);
1184  ASSERT_STRING_EQUAL(cu_error_msg,
1185  "Both input geometries must be linestrings");
1186 
1187  /* Invalid input, too short linestring */
1188 
1189  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1190  g2 = lwgeom_from_wkt("LINESTRING M(2 0 1)", LW_PARSER_CHECK_NONE);
1191  dist = -77;
1192  m = lwgeom_tcpa(g1, g2, &dist);
1193  lwgeom_free(g1);
1194  lwgeom_free(g2);
1195  ASSERT_DOUBLE_EQUAL(dist, -77.0); /* not touched */
1196  ASSERT_DOUBLE_EQUAL(m, -1.0);
1197  ASSERT_STRING_EQUAL(
1198  cu_error_msg,
1199  "Both input lines must have at least 2 points" /* should be accepted
1200  ? */
1201 
1202  );
1203 
1204  /* Invalid input, empty linestring */
1205 
1206  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1207  g2 = lwgeom_from_wkt("LINESTRING M EMPTY", LW_PARSER_CHECK_NONE);
1208  m = lwgeom_tcpa(g1, g2, NULL);
1209  lwgeom_free(g1);
1210  lwgeom_free(g2);
1211  ASSERT_DOUBLE_EQUAL(m, -1.0);
1212  ASSERT_STRING_EQUAL(cu_error_msg,
1213  "Both input lines must have at least 2 points"
1214 
1215  );
1216 
1217  /* Timeranges do not overlap */
1218 
1219  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1220  g2 = lwgeom_from_wkt("LINESTRING M(0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
1221  m = lwgeom_tcpa(g1, g2, NULL);
1222  lwgeom_free(g1);
1223  lwgeom_free(g2);
1224  ASSERT_DOUBLE_EQUAL(m, -2.0); /* means timeranges do not overlap */
1225 
1226  /* One of the tracks is still, the other passes to that point */
1227 
1228  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1229  g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
1230  dist = -1;
1231  m = lwgeom_tcpa(g1, g2, &dist);
1232  ASSERT_DOUBLE_EQUAL(m, 1.0);
1233  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1234  CU_ASSERT( lwgeom_cpa_within(g1, g2, 0.0) == LW_TRUE );
1235  lwgeom_free(g1);
1236  lwgeom_free(g2);
1237 
1238  /* One of the tracks is still, the other passes at 10 meters from point */
1239 
1240  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
1241  g2 = lwgeom_from_wkt("LINESTRING M(-10 10 1, 10 10 5)", LW_PARSER_CHECK_NONE);
1242  dist = -1;
1243  m = lwgeom_tcpa(g1, g2, &dist);
1244  ASSERT_DOUBLE_EQUAL(m, 3.0);
1245  ASSERT_DOUBLE_EQUAL(dist, 10.0);
1246  CU_ASSERT( lwgeom_cpa_within(g1, g2, 11.0) == LW_TRUE );
1247  CU_ASSERT( lwgeom_cpa_within(g1, g2, 10.0) == LW_TRUE );
1248  CU_ASSERT( lwgeom_cpa_within(g1, g2, 9.0) == LW_FALSE );
1249  lwgeom_free(g1);
1250  lwgeom_free(g2);
1251 
1252  /* Equal tracks, 2d */
1253 
1254  g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1255  g2 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1256  dist = -1;
1257  m = lwgeom_tcpa(g1, g2, &dist);
1258  lwgeom_free(g1);
1259  lwgeom_free(g2);
1260  ASSERT_DOUBLE_EQUAL(m, 10.0);
1261  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1262 
1263  /* Reversed tracks, 2d */
1264 
1265  g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1266  g2 = lwgeom_from_wkt("LINESTRING M(10 0 10, 0 0 20)", LW_PARSER_CHECK_NONE);
1267  dist = -1;
1268  m = lwgeom_tcpa(g1, g2, &dist);
1269  lwgeom_free(g1);
1270  lwgeom_free(g2);
1271  ASSERT_DOUBLE_EQUAL(m, 15.0);
1272  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1273 
1274  /* Parallel tracks, same speed, 2d */
1275 
1276  g1 = lwgeom_from_wkt("LINESTRING M(2 0 10, 12 0 20)", LW_PARSER_CHECK_NONE);
1277  g2 = lwgeom_from_wkt("LINESTRING M(13 0 10, 23 0 20)", LW_PARSER_CHECK_NONE);
1278  dist = -1;
1279  m = lwgeom_tcpa(g1, g2, &dist);
1280  lwgeom_free(g1);
1281  lwgeom_free(g2);
1282  ASSERT_DOUBLE_EQUAL(m, 10.0);
1283  ASSERT_DOUBLE_EQUAL(dist, 11.0);
1284 
1285  /* Parallel tracks, different speed (g2 gets closer as time passes), 2d */
1286 
1287  g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1288  g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 9 0 20)", LW_PARSER_CHECK_NONE);
1289  dist = -1;
1290  m = lwgeom_tcpa(g1, g2, &dist);
1291  lwgeom_free(g1);
1292  lwgeom_free(g2);
1293  ASSERT_DOUBLE_EQUAL(m, 20.0);
1294  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1295 
1296  /* Parallel tracks, different speed (g2 left behind as time passes), 2d */
1297 
1298  g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1299  g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 6 0 20)", LW_PARSER_CHECK_NONE);
1300  dist = -1;
1301  m = lwgeom_tcpa(g1, g2, &dist);
1302  lwgeom_free(g1);
1303  lwgeom_free(g2);
1304  ASSERT_DOUBLE_EQUAL(m, 10.0);
1305  ASSERT_DOUBLE_EQUAL(dist, 2.0);
1306 
1307  /* Tracks, colliding, 2d */
1308 
1309  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
1310  g2 = lwgeom_from_wkt("LINESTRING M(5 -8 0, 5 8 10)", LW_PARSER_CHECK_NONE);
1311  dist = -1;
1312  m = lwgeom_tcpa(g1, g2, &dist);
1313  lwgeom_free(g1);
1314  lwgeom_free(g2);
1315  ASSERT_DOUBLE_EQUAL(m, 5.0);
1316  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1317 
1318  /* Tracks crossing, NOT colliding, 2d */
1319 
1320  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
1321  g2 = lwgeom_from_wkt("LINESTRING M(8 -5 0, 8 5 10)", LW_PARSER_CHECK_NONE);
1322  dist = -1;
1323  m = lwgeom_tcpa(g1, g2, &dist);
1324  lwgeom_free(g1);
1325  lwgeom_free(g2);
1326  ASSERT_DOUBLE_EQUAL(m, 6.5);
1327  ASSERT_DOUBLE_EQUAL(rint(dist*100), 212.0);
1328 
1329  /* Same origin, different direction, 2d */
1330 
1331  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1332  g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, -100 0 10)", LW_PARSER_CHECK_NONE);
1333  dist = -1;
1334  m = lwgeom_tcpa(g1, g2, &dist);
1335  lwgeom_free(g1);
1336  lwgeom_free(g2);
1337  ASSERT_DOUBLE_EQUAL(m, 1.0);
1338  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1339 
1340  /* Same ending, different direction, 2d */
1341 
1342  g1 = lwgeom_from_wkt("LINESTRING M(10 0 1, 0 0 10)", LW_PARSER_CHECK_NONE);
1343  g2 = lwgeom_from_wkt("LINESTRING M(0 -100 1, 0 0 10)", LW_PARSER_CHECK_NONE);
1344  dist = -1;
1345  m = lwgeom_tcpa(g1, g2, &dist);
1346  lwgeom_free(g1);
1347  lwgeom_free(g2);
1348  ASSERT_DOUBLE_EQUAL(m, 10.0);
1349  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1350 
1351  /* Converging tracks, 3d */
1352 
1353  g1 = lwgeom_from_wkt("LINESTRING ZM(0 0 0 10, 10 0 0 20)", LW_PARSER_CHECK_NONE);
1354  g2 = lwgeom_from_wkt("LINESTRING ZM(0 0 8 10, 10 0 5 20)", LW_PARSER_CHECK_NONE);
1355  dist = -1;
1356  m = lwgeom_tcpa(g1, g2, &dist);
1357  lwgeom_free(g1);
1358  lwgeom_free(g2);
1359  ASSERT_DOUBLE_EQUAL(m, 20.0);
1360  ASSERT_DOUBLE_EQUAL(dist, 5.0);
1361 
1362  /* G1 stops at t=1 until t=4 to let G2 pass by, then continues */
1363  /* G2 passes at 1 meter from G1 t=3 */
1364 
1365  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 0 1 1, 0 1 4, 0 10 13)", LW_PARSER_CHECK_NONE);
1366  g2 = lwgeom_from_wkt("LINESTRING M(-10 2 0, 0 2 3, 12 2 13)", LW_PARSER_CHECK_NONE);
1367  dist = -1;
1368  m = lwgeom_tcpa(g1, g2, &dist);
1369  lwgeom_free(g1);
1370  lwgeom_free(g2);
1371  ASSERT_DOUBLE_EQUAL(m, 3.0);
1372  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1373 
1374  /* Test for https://trac.osgeo.org/postgis/ticket/3136 */
1375 
1376  g1 = lwgeom_from_wkt("LINESTRING M (0 0 1432291464,2 0 1432291536) ", LW_PARSER_CHECK_NONE);
1377  g2 = lwgeom_from_wkt("LINESTRING M (0 0 1432291464,1 0 1432291466.25,2 0 1432291500)", LW_PARSER_CHECK_NONE);
1378  dist = -1;
1379  m = lwgeom_tcpa(g1, g2, &dist);
1380  lwgeom_free(g1);
1381  lwgeom_free(g2);
1382  ASSERT_DOUBLE_EQUAL(m, 1432291464.0);
1383  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1384 
1385  /* Tracks share a single point in time */
1386 
1387  g1 = lwgeom_from_wkt("LINESTRINGM(0 0 0, 1 0 2)", LW_PARSER_CHECK_NONE);
1388  g2 = lwgeom_from_wkt("LINESTRINGM(0 0 2, 1 0 3)", LW_PARSER_CHECK_NONE);
1389  dist = -1;
1390  m = lwgeom_tcpa(g1, g2, &dist);
1391  lwgeom_free(g1);
1392  lwgeom_free(g2);
1393  ASSERT_DOUBLE_EQUAL(m, 2.0);
1394  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1395 
1396 }
1397 
1398 static void
1399 test_lwgeom_is_trajectory(void)
1400 {
1401  LWGEOM *g;
1402  int ret;
1403 
1404  g = lwgeom_from_wkt("POINT M(0 0 1)", LW_PARSER_CHECK_NONE);
1405  ret = lwgeom_is_trajectory(g);
1406  lwgeom_free(g);
1407  ASSERT_INT_EQUAL(ret, LW_FALSE); /* not a linestring */
1408 
1409  g = lwgeom_from_wkt("LINESTRING Z(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1410  ret = lwgeom_is_trajectory(g);
1411  lwgeom_free(g);
1412  ASSERT_INT_EQUAL(ret, LW_FALSE); /* no measure */
1413 
1414  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1415  ret = lwgeom_is_trajectory(g);
1416  lwgeom_free(g);
1417  ASSERT_INT_EQUAL(ret, LW_FALSE); /* same measure in two points */
1418 
1419  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 0)", LW_PARSER_CHECK_NONE);
1420  ret = lwgeom_is_trajectory(g);
1421  lwgeom_free(g);
1422  ASSERT_INT_EQUAL(ret, LW_FALSE); /* backward measure */
1423 
1424  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 1 0 3, 2 2 2)", LW_PARSER_CHECK_NONE);
1425  ret = lwgeom_is_trajectory(g);
1426  lwgeom_free(g);
1427  ASSERT_INT_EQUAL(ret, LW_FALSE); /* backward measure */
1428 
1429  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 2)", LW_PARSER_CHECK_NONE);
1430  ret = lwgeom_is_trajectory(g);
1431  lwgeom_free(g);
1432  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (still) */
1433 
1434  g = lwgeom_from_wkt("LINESTRING M EMPTY", LW_PARSER_CHECK_NONE);
1435  ret = lwgeom_is_trajectory(g);
1436  lwgeom_free(g);
1437  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (empty) */
1438 
1439  g = lwgeom_from_wkt("LINESTRING M (0 0 1)", LW_PARSER_CHECK_NONE);
1440  ret = lwgeom_is_trajectory(g);
1441  lwgeom_free(g);
1442  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (corner case) */
1443 }
1444 
1445 /*
1446 ** Used by test harness to register the tests in this file.
1447 */
1448 void measures_suite_setup(void);
1449 void measures_suite_setup(void)
1450 {
1451  CU_pSuite suite = CU_add_suite("measures", NULL, NULL);
1452  PG_ADD_TEST(suite, test_mindistance2d_tolerance);
1453  PG_ADD_TEST(suite, test_mindistance3d_tolerance);
1454  PG_ADD_TEST(suite, test_rect_tree_contains_point);
1455  PG_ADD_TEST(suite, test_rect_tree_intersects_tree);
1456  PG_ADD_TEST(suite, test_lwgeom_segmentize2d);
1457  PG_ADD_TEST(suite, test_lwgeom_locate_along);
1458  PG_ADD_TEST(suite, test_lw_dist2d_pt_arc);
1459  PG_ADD_TEST(suite, test_lw_dist2d_seg_arc);
1460  PG_ADD_TEST(suite, test_lw_dist2d_arc_arc);
1461  PG_ADD_TEST(suite, test_lw_arc_length);
1462  PG_ADD_TEST(suite, test_lw_dist2d_pt_ptarrayarc);
1463  PG_ADD_TEST(suite, test_lw_dist2d_ptarray_ptarrayarc);
1464  PG_ADD_TEST(suite, test_lwgeom_tcpa);
1465  PG_ADD_TEST(suite, test_lwgeom_is_trajectory);
1466  PG_ADD_TEST(suite, test_rect_tree_distance_tree);
1467 }
r
char * r
Definition: cu_in_wkt.c:24
test_lw_dist2d_ptarray_ptarrayarc
static void test_lw_dist2d_ptarray_ptarrayarc(void)
Definition: cu_measures.c:1108
do_test_mindistance_tolerance
static void do_test_mindistance_tolerance(char *in1, char *in2, double expected_res, int line, double(*distancef)(const LWGEOM *, const LWGEOM *, double))
Definition: cu_measures.c:40
test_lw_arc_length
static void test_lw_arc_length(void)
Definition: cu_measures.c:983
test_rect_tree_intersects_tree
static void test_rect_tree_intersects_tree(void)
Definition: cu_measures.c:416
test_mindistance2d_tolerance
static void test_mindistance2d_tolerance(void)
Definition: cu_measures.c:82
test_lwgeom_is_trajectory
static void test_lwgeom_is_trajectory(void)
Definition: cu_measures.c:1399
measures_suite_setup
void measures_suite_setup(void)
Definition: cu_measures.c:1449
tree_inter
static int tree_inter(const char *wkt1, const char *wkt2)
Definition: cu_measures.c:402
test_lwgeom_tcpa
static void test_lwgeom_tcpa(void)
Definition: cu_measures.c:1159
test_lw_dist2d_seg_arc
static void test_lw_dist2d_seg_arc(void)
Definition: cu_measures.c:715
test_rect_tree_distance_tree
static void test_rect_tree_distance_tree(void)
Definition: cu_measures.c:502
TDT
#define TDT(w1, w2, d)
Definition: cu_measures.c:499
test_lw_dist2d_pt_ptarrayarc
static void test_lw_dist2d_pt_ptarrayarc(void)
Definition: cu_measures.c:1025
DIST2DTEST
#define DIST2DTEST(str1, str2, res)
Definition: cu_measures.c:33
DIST3DTEST
#define DIST3DTEST(str1, str2, res)
Definition: cu_measures.c:36
test_lw_dist2d_arc_arc
static void test_lw_dist2d_arc_arc(void)
Definition: cu_measures.c:795
test_rect_tree_contains_point
static void test_rect_tree_contains_point(void)
Definition: cu_measures.c:255
test_lwgeom_locate_along
static void test_lwgeom_locate_along(void)
Definition: cu_measures.c:615
test_lwgeom_segmentize2d
static void test_lwgeom_segmentize2d(void)
Definition: cu_measures.c:563
lwgeom_from_text
static LWGEOM * lwgeom_from_text(const char *str)
Definition: cu_measures.c:25
test_lw_dist2d_pt_arc
static void test_lw_dist2d_pt_arc(void)
Definition: cu_measures.c:642
test_rect_tree_distance_tree_case
static double test_rect_tree_distance_tree_case(const char *wkt1, const char *wkt2)
Definition: cu_measures.c:477
test_mindistance3d_tolerance
static void test_mindistance3d_tolerance(void)
Definition: cu_measures.c:213
tree_pt
static int tree_pt(RECT_NODE *tree, double x, double y)
Definition: cu_measures.c:248
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_DOUBLE_EQUAL
#define ASSERT_DOUBLE_EQUAL(o, e)
Definition: liblwgeom/cunit/cu_tester.h:31
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
lwgeom_as_lwline
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:170
LW_FALSE
#define LW_FALSE
Definition: liblwgeom.h:77
lwgeom_request_interrupt
void lwgeom_request_interrupt(void)
Request interruption of any running code.
Definition: lwgeom_api.c:728
lwgeom_locate_along
LWGEOM * lwgeom_locate_along(const LWGEOM *lwin, double m, double offset)
Determine the location(s) along a measured line where m occurs and return as a multipoint.
Definition: lwlinearreferencing.c:217
LW_FAILURE
#define LW_FAILURE
Definition: liblwgeom.h:79
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_segmentize2d
LWGEOM * lwgeom_segmentize2d(const LWGEOM *line, double dist)
Definition: lwgeom.c:762
LW_SUCCESS
#define LW_SUCCESS
Definition: liblwgeom.h:80
lwgeom_is_trajectory
int lwgeom_is_trajectory(const LWGEOM *geom)
Return LW_TRUE or LW_FALSE depending on whether or not a geometry is a linestring with measure value ...
Definition: lwgeom.c:2472
lwgeom_parse_wkt
int lwgeom_parse_wkt(LWGEOM_PARSER_RESULT *parser_result, char *wktstr, int parse_flags)
Parse a WKT geometry string into an LWGEOM structure.
Definition: lwin_wkt_parse.c:126
lwgeom_to_ewkt
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:556
lwgeom_tcpa
double lwgeom_tcpa(const LWGEOM *g1, const LWGEOM *g2, double *mindist)
Find the time of closest point of approach.
Definition: lwlinearreferencing.c:1077
lwfree
void lwfree(void *mem)
Definition: lwutil.c:244
WKT_ISO
#define WKT_ISO
Definition: liblwgeom.h:2075
lwgeom_to_wkt
char * lwgeom_to_wkt(const LWGEOM *geom, uint8_t variant, int precision, size_t *size_out)
WKT emitter function.
Definition: lwout_wkt.c:676
lwgeom_from_wkt
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition: lwin_wkt.c:904
LW_TRUE
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:76
lwline_free
void lwline_free(LWLINE *line)
Definition: lwline.c:76
lwgeom_cpa_within
int lwgeom_cpa_within(const LWGEOM *g1, const LWGEOM *g2, double maxdist)
Is the closest point of approach within a distance ?
Definition: lwlinearreferencing.c:1243
lw_arc_length
double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
Returns the length of a circular arc segment.
Definition: lwalgorithm.c:118
liblwgeom_internal.h
distance
Datum distance(PG_FUNCTION_ARGS)
Definition: lwgeom_backend_api.c:199
rect_tree_from_lwgeom
RECT_NODE * rect_tree_from_lwgeom(const LWGEOM *lwgeom)
Create a tree index on top an LWGEOM.
Definition: lwtree.c:861
rect_tree_contains_point
int rect_tree_contains_point(RECT_NODE *node, const POINT2D *pt)
Definition: lwtree.c:398
rect_tree_free
void rect_tree_free(RECT_NODE *node)
Recurse from top of node tree and free all children.
Definition: lwtree.c:69
rect_tree_distance_tree
double rect_tree_distance_tree(RECT_NODE *n1, RECT_NODE *n2, double threshold)
Return the distance between two RECT_NODE trees.
Definition: lwtree.c:1354
rect_tree_intersects_tree
int rect_tree_intersects_tree(RECT_NODE *n1, RECT_NODE *n2)
Test if two RECT_NODE trees intersect one another.
Definition: lwtree.c:990
lw_dist2d_pt_arc
int lw_dist2d_pt_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, DISTPTS *dl)
Definition: measures.c:1439
lw_dist2d_seg_arc
int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl)
Calculate the shortest distance between an arc and an edge.
Definition: measures.c:1292
lw_dist2d_distpts_init
void lw_dist2d_distpts_init(DISTPTS *dl, int mode)
Definition: measures.c:67
lw_dist2d_arc_arc
int lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl)
Definition: measures.c:1497
lw_dist2d_ptarray_ptarrayarc
int lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl)
Test each segment of pa against each arc of pb for distance.
Definition: measures.c:1187
lw_dist2d_pt_ptarrayarc
int lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl)
Search all the arcs of pointarray to see which one is closest to p1 Returns minimum distance between ...
Definition: measures.c:1083
DIST_MIN
#define DIST_MIN
Definition: measures.h:44
cu_tester.h
DISTPTS::distance
double distance
Definition: measures.h:51
DISTPTS
Structure used in distance-calculations.
Definition: measures.h:50
LWLINE::points
POINTARRAY * points
Definition: liblwgeom.h:425
POINT2D::y
double y
Definition: liblwgeom.h:331
POINT2D::x
double x
Definition: liblwgeom.h:331
rect_node
Definition: lwtree.h:66
struct_lwgeom_parser_result
Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM.
Definition: liblwgeom.h:2013