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