PostGIS  2.2.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@keybit.net>
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 "lwtree.h"
23 
24 static LWGEOM* lwgeom_from_text(const char *str)
25 {
27  if( LW_FAILURE == lwgeom_parse_wkt(&r, (char*)str, LW_PARSER_CHECK_NONE) )
28  return NULL;
29  return r.geom;
30 }
31 
32 #define DIST2DTEST(str1, str2, res) do_test_mindistance2d_tolerance(str1, str2, res, __LINE__)
33 
34 static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expected_res, int line)
35 {
36  LWGEOM *lw1;
37  LWGEOM *lw2;
38  double distance;
39  char *msg1 = "test_mindistance2d_tolerance failed (got %g expected %g) at line %d\n";
40  char *msg2 = "\n\ndo_test_mindistance2d_tolerance: NULL lwgeom generated from WKT\n %s\n\n";
41 
44 
45  if ( ! lw1 )
46  {
47  printf(msg2, in1);
48  exit(1);
49  }
50  if ( ! lw2 )
51  {
52  printf(msg2, in2);
53  exit(1);
54  }
55 
56  distance = lwgeom_mindistance2d_tolerance(lw1, lw2, 0.0);
57  lwgeom_free(lw1);
58  lwgeom_free(lw2);
59 
60  if ( fabs(distance - expected_res) > 0.00001 )
61  {
62  printf(msg1, distance, expected_res, line);
63  CU_FAIL();
64  }
65  else
66  {
67  CU_PASS();
68  }
69 
70 }
71 
73 {
74  /*
75  ** Simple case.
76  */
77  DIST2DTEST("POINT(0 0)", "MULTIPOINT(0 1.5,0 2,0 2.5)", 1.5);
78 
79  /*
80  ** Point vs Geometry Collection.
81  */
82  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
83 
84  /*
85  ** Point vs Geometry Collection Collection.
86  */
87  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
88 
89  /*
90  ** Point vs Geometry Collection Collection Collection.
91  */
92  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4))))", 5.0);
93 
94  /*
95  ** Point vs Geometry Collection Collection Collection Multipoint.
96  */
97  DIST2DTEST("POINT(0 0)", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4))))", 5.0);
98 
99  /*
100  ** Geometry Collection vs Geometry Collection
101  */
102  DIST2DTEST("GEOMETRYCOLLECTION(POINT(0 0))", "GEOMETRYCOLLECTION(POINT(3 4))", 5.0);
103 
104  /*
105  ** Geometry Collection Collection vs Geometry Collection Collection
106  */
107  DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(POINT(3 4)))", 5.0);
108 
109  /*
110  ** Geometry Collection Collection Multipoint vs Geometry Collection Collection Multipoint
111  */
112  DIST2DTEST("GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(0 0)))", "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION(MULTIPOINT(3 4)))", 5.0);
113 
114  /*
115  ** Linestring vs its start point
116  */
117  DIST2DTEST("LINESTRING(-2 0, -0.2 0)", "POINT(-2 0)", 0);
118 
119  /*
120  ** Linestring vs its end point
121  */
122  DIST2DTEST("LINESTRING(-0.2 0, -2 0)", "POINT(-2 0)", 0);
123 
124  /*
125  ** Linestring vs its start point (tricky number, see #1459)
126  */
127  DIST2DTEST("LINESTRING(-1e-8 0, -0.2 0)", "POINT(-1e-8 0)", 0);
128 
129  /*
130  ** Linestring vs its end point (tricky number, see #1459)
131  */
132  DIST2DTEST("LINESTRING(-0.2 0, -1e-8 0)", "POINT(-1e-8 0)", 0);
133 
134  /*
135  * Circular string and point
136  */
137  DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "POINT(0 0)", 1);
138  DIST2DTEST("CIRCULARSTRING(-3 0, -2 0, -1 0, 0 1, 1 0)", "POINT(0 0)", 1);
139 
140  /*
141  * Circular string and Circular string
142  */
143  DIST2DTEST("CIRCULARSTRING(-1 0, 0 1, 1 0)", "CIRCULARSTRING(0 0, 1 -1, 2 0)", 1);
144 
145  /*
146  * CurvePolygon and Point
147  */
148  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)))";
149  DIST2DTEST(cs1, "POINT(3 14)", 1);
150  DIST2DTEST(cs1, "POINT(3 8)", 0);
151  DIST2DTEST(cs1, "POINT(6 5)", 1);
152  DIST2DTEST(cs1, "POINT(6 4)", 0);
153 
154  /*
155  * CurvePolygon and Linestring
156  */
157  DIST2DTEST(cs1, "LINESTRING(0 0, 50 0)", 0.917484);
158  DIST2DTEST(cs1, "LINESTRING(6 0, 10 7)", 0);
159  DIST2DTEST(cs1, "LINESTRING(4 4, 4 8)", 0);
160  DIST2DTEST(cs1, "LINESTRING(4 7, 5 6, 6 7)", 0.585786);
161  DIST2DTEST(cs1, "LINESTRING(10 0, 10 2, 10 0)", 1.52913);
162 
163  /*
164  * CurvePolygon and Polygon
165  */
166  DIST2DTEST(cs1, "POLYGON((10 4, 10 8, 13 8, 13 4, 10 4))", 0.58415);
167  DIST2DTEST(cs1, "POLYGON((9 4, 9 8, 12 8, 12 4, 9 4))", 0);
168  DIST2DTEST(cs1, "POLYGON((1 4, 1 8, 4 8, 4 4, 1 4))", 0);
169 
170  /*
171  * CurvePolygon and CurvePolygon
172  */
173  DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(-1 4, 0 5, 1 4, 0 3, -1 4))", 0.0475666);
174  DIST2DTEST(cs1, "CURVEPOLYGON(CIRCULARSTRING(1 4, 2 5, 3 4, 2 3, 1 4))", 0.0);
175 
176  /*
177  * MultiSurface and CurvePolygon
178  */
179  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)))";
180  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))", 1);
181  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))", 0);
182  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))", 0);
183  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))", 0);
184  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))", 0.5);
185  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))", 0);
186  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))", 0);
187  DIST2DTEST(cs2, "CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))", 2.60555);
188 
189  /*
190  * MultiCurve and Linestring
191  */
192  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);
193 
194 }
195 
197 {
198  LWPOLY *poly;
199  POINT2D p;
200  RECT_NODE* tree;
201  int result;
202  int boundary = 0;
203 
204  /* square */
205  poly = (LWPOLY*)lwgeom_from_wkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", LW_PARSER_CHECK_NONE);
206  tree = rect_tree_new(poly->rings[0]);
207 
208  /* inside square */
209  boundary = 0;
210  p.x = 0.5;
211  p.y = 0.5;
212  result = rect_tree_contains_point(tree, &p, &boundary);
213  CU_ASSERT_NOT_EQUAL(result, 0);
214 
215  /* outside square */
216  boundary = 0;
217  p.x = 1.5;
218  p.y = 0.5;
219  result = rect_tree_contains_point(tree, &p, &boundary);
220  CU_ASSERT_EQUAL(result, 0);
221 
222  rect_tree_free(tree);
223  lwpoly_free(poly);
224 
225  /* ziggy zaggy horizontal saw tooth polygon */
226  poly = (LWPOLY*)lwgeom_from_wkt("POLYGON((0 0, 1 3, 2 0, 3 3, 4 0, 4 5, 0 5, 0 0))", LW_PARSER_CHECK_NONE);
227  tree = rect_tree_new(poly->rings[0]);
228 
229  /* not in, left side */
230  boundary = 0;
231  p.x = -0.5;
232  p.y = 0.5;
233  result = rect_tree_contains_point(tree, &p, &boundary);
234  CU_ASSERT_EQUAL(result, 0);
235 
236  /* not in, right side */
237  boundary = 0;
238  p.x = 3.0;
239  p.y = 1.0;
240  result = rect_tree_contains_point(tree, &p, &boundary);
241  CU_ASSERT_EQUAL(result, 0);
242 
243  /* inside */
244  boundary = 0;
245  p.x = 2.0;
246  p.y = 1.0;
247  result = rect_tree_contains_point(tree, &p, &boundary);
248  CU_ASSERT_NOT_EQUAL(result, 0);
249 
250  /* on left border */
251  boundary = 0;
252  p.x = 0.0;
253  p.y = 1.0;
254  result = rect_tree_contains_point(tree, &p, &boundary);
255  CU_ASSERT_EQUAL(boundary, 1);
256 
257  /* on right border */
258  boundary = 0;
259  p.x = 4.0;
260  p.y = 0.0;
261  result = rect_tree_contains_point(tree, &p, &boundary);
262  CU_ASSERT_EQUAL(boundary, 1);
263 
264  /* on tooth concave */
265  boundary = 0;
266  p.x = 3.0;
267  p.y = 3.0;
268  result = rect_tree_contains_point(tree, &p, &boundary);
269  CU_ASSERT_EQUAL(boundary, 1);
270 
271  /* on tooth convex */
272  boundary = 0;
273  p.x = 2.0;
274  p.y = 0.0;
275  result = rect_tree_contains_point(tree, &p, &boundary);
276  CU_ASSERT_EQUAL(boundary, 1);
277 
278  rect_tree_free(tree);
279  lwpoly_free(poly);
280 
281  /* ziggy zaggy vertical saw tooth polygon */
282  poly = (LWPOLY*)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);
283  tree = rect_tree_new(poly->rings[0]);
284 
285  /* not in, left side */
286  boundary = 0;
287  p.x = -0.5;
288  p.y = 3.5;
289  result = rect_tree_contains_point(tree, &p, &boundary);
290  CU_ASSERT_EQUAL(result, 0);
291 
292  /* not in, right side */
293  boundary = 0;
294  p.x = 6.0;
295  p.y = 2.2;
296  result = rect_tree_contains_point(tree, &p, &boundary);
297  CU_ASSERT_EQUAL(result, 0);
298 
299  /* inside */
300  boundary = 0;
301  p.x = 3.0;
302  p.y = 2.0;
303  result = rect_tree_contains_point(tree, &p, &boundary);
304  CU_ASSERT_NOT_EQUAL(result, 0);
305 
306  /* on bottom border */
307  boundary = 0;
308  p.x = 1.0;
309  p.y = 0.0;
310  result = rect_tree_contains_point(tree, &p, &boundary);
311  CU_ASSERT_EQUAL(boundary, 1);
312 
313  /* on top border */
314  boundary = 0;
315  p.x = 3.0;
316  p.y = 6.0;
317  result = rect_tree_contains_point(tree, &p, &boundary);
318  CU_ASSERT_EQUAL(boundary, 1);
319 
320  /* on tooth concave */
321  boundary = 0;
322  p.x = 3.0;
323  p.y = 1.0;
324  result = rect_tree_contains_point(tree, &p, &boundary);
325  CU_ASSERT_EQUAL(boundary, 1);
326 
327  /* on tooth convex */
328  boundary = 0;
329  p.x = 0.0;
330  p.y = 2.0;
331  result = rect_tree_contains_point(tree, &p, &boundary);
332  CU_ASSERT_EQUAL(boundary, 1);
333 
334  /* on tooth convex */
335  boundary = 0;
336  p.x = 0.0;
337  p.y = 6.0;
338  result = rect_tree_contains_point(tree, &p, &boundary);
339  CU_ASSERT_EQUAL(boundary, 1);
340 
341  rect_tree_free(tree);
342  lwpoly_free(poly);
343 
344 }
345 
347 {
348  LWPOLY *poly1, *poly2;
349  RECT_NODE *tree1, *tree2;
350  int result;
351 
352  /* total overlap, A == B */
353  poly1 = (LWPOLY*)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);
354  poly2 = (LWPOLY*)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);
355  tree1 = rect_tree_new(poly1->rings[0]);
356  tree2 = rect_tree_new(poly2->rings[0]);
357  result = rect_tree_intersects_tree(tree1, tree2);
358  CU_ASSERT_EQUAL(result, LW_TRUE);
359  lwpoly_free(poly1);
360  lwpoly_free(poly2);
361  rect_tree_free(tree1);
362  rect_tree_free(tree2);
363 
364  /* hiding between the tines of the comb */
365  poly1 = (LWPOLY*)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);
366  poly2 = (LWPOLY*)lwgeom_from_wkt("POLYGON((0.3 0.7, 0.3 0.8, 0.4 0.8, 0.4 0.7, 0.3 0.7))", LW_PARSER_CHECK_NONE);
367  tree1 = rect_tree_new(poly1->rings[0]);
368  tree2 = rect_tree_new(poly2->rings[0]);
369  result = rect_tree_intersects_tree(tree1, tree2);
370  CU_ASSERT_EQUAL(result, LW_FALSE);
371  lwpoly_free(poly1);
372  lwpoly_free(poly2);
373  rect_tree_free(tree1);
374  rect_tree_free(tree2);
375 
376  /* between the tines, but with a corner overlapping */
377  poly1 = (LWPOLY*)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  poly2 = (LWPOLY*)lwgeom_from_wkt("POLYGON((0.3 0.7, 0.3 0.8, 0.4 0.8, 1.3 0.3, 0.3 0.7))", LW_PARSER_CHECK_NONE);
379  tree1 = rect_tree_new(poly1->rings[0]);
380  tree2 = rect_tree_new(poly2->rings[0]);
381  result = rect_tree_intersects_tree(tree1, tree2);
382  CU_ASSERT_EQUAL(result, LW_TRUE);
383  lwpoly_free(poly1);
384  lwpoly_free(poly2);
385  rect_tree_free(tree1);
386  rect_tree_free(tree2);
387 
388  /* Just touching the top left corner of the comb */
389  poly1 = (LWPOLY*)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);
390  poly2 = (LWPOLY*)lwgeom_from_wkt("POLYGON((-1 5, 0 5, 0 7, -1 7, -1 5))", LW_PARSER_CHECK_NONE);
391  tree1 = rect_tree_new(poly1->rings[0]);
392  tree2 = rect_tree_new(poly2->rings[0]);
393  result = rect_tree_intersects_tree(tree1, tree2);
394  CU_ASSERT_EQUAL(result, LW_TRUE);
395  lwpoly_free(poly1);
396  lwpoly_free(poly2);
397  rect_tree_free(tree1);
398  rect_tree_free(tree2);
399 
400 
401 }
402 
403 static void
405 {
406  LWGEOM *linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
407  LWGEOM *lineout = lwgeom_segmentize2d(linein, 5);
408  char *strout = lwgeom_to_ewkt(lineout);
409  CU_ASSERT_STRING_EQUAL(strout, "LINESTRING(0 0,5 0,10 0)");
410  lwfree(strout);
411  lwgeom_free(linein);
412  lwgeom_free(lineout);
413 
414  /* test interruption */
415 
416  linein = lwgeom_from_wkt("LINESTRING(0 0,10 0)", LW_PARSER_CHECK_NONE);
418  lineout = lwgeom_segmentize2d(linein, 1e-100);
419  CU_ASSERT_EQUAL(lineout, NULL);
420  lwgeom_free(linein);
421 
422  linein = lwgeom_from_wkt("MULTILINESTRING((0 0,10 0),(20 0, 30 0))", LW_PARSER_CHECK_NONE);
424  lineout = lwgeom_segmentize2d(linein, 1e-100);
425  CU_ASSERT_EQUAL(lineout, NULL);
426  lwgeom_free(linein);
427 
428  linein = lwgeom_from_wkt(
429  "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)))"
432  lineout = lwgeom_segmentize2d(linein, 1e-100);
433  CU_ASSERT_EQUAL(lineout, NULL);
434  lwgeom_free(linein);
435 
436  linein = lwgeom_from_wkt(
437  "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))"
439  CU_ASSERT_FATAL(linein != NULL);
441  lineout = lwgeom_segmentize2d(linein, 1e-100);
442  CU_ASSERT_EQUAL(lineout, NULL);
443  lwgeom_free(linein);
444 
445  linein = lwgeom_from_wkt("LINESTRING(20 0, 30 0)", LW_PARSER_CHECK_NONE);
446  /* NOT INTERRUPTED */
447  lineout = lwgeom_segmentize2d(linein, 5);
448  strout = lwgeom_to_ewkt(lineout);
449  CU_ASSERT_STRING_EQUAL(strout, "LINESTRING(20 0,25 0,30 0)");
450  lwfree(strout);
451  lwgeom_free(linein);
452  lwgeom_free(lineout);
453 }
454 
455 static void
457 {
458  LWGEOM *geom = NULL;
459  LWGEOM *out = NULL;
460  double measure = 105.0;
461  char *str;
462 
463  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
464  geom = lwgeom_from_wkt("MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))", LW_PARSER_CHECK_NONE);
465  out = lwgeom_locate_along(geom, measure, 0.0);
466  str = lwgeom_to_wkt(out, WKT_ISO, 8, NULL);
467  lwgeom_free(geom);
468  lwgeom_free(out);
469  CU_ASSERT_STRING_EQUAL("MULTIPOINT M (55.226131 55.226131 105)", str);
470  lwfree(str);
471 
472  /* ST_Locatealong(ST_GeomFromText('MULTILINESTRING M ((1 2 3, 5 4 5), (50 50 1, 60 60 200))'), 105) */
473  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);
474  out = lwgeom_locate_along(geom, measure, 0.0);
475  str = lwgeom_to_wkt(out, WKT_ISO, 8, NULL);
476  lwgeom_free(geom);
477  lwgeom_free(out);
478  CU_ASSERT_STRING_EQUAL("MULTIPOINT M (55.226131 55.226131 105)", str);
479  lwfree(str);
480 }
481 
482 static void
484 {
485  /* int lw_dist2d_pt_arc(const POINT2D* P, const POINT2D* A1, const POINT2D* A2, const POINT2D* A3, DISTPTS* dl) */
486  DISTPTS dl;
487  POINT2D P, A1, A2, A3;
488  int rv;
489 
490 
491  /* Point within unit semicircle, 0.5 units from arc */
492  A1.x = -1; A1.y = 0;
493  A2.x = 0 ; A2.y = 1;
494  A3.x = 1 ; A3.y = 0;
495  P.x = 0 ; P.y = 0.5;
496 
498  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
499  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
500  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
501 
502  /* Point outside unit semicircle, 0.5 units from arc */
503  P.x = 0 ; P.y = 1.5;
505  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
506  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
507  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
508 
509  /* Point outside unit semicircle, sqrt(2) units from arc end point*/
510  P.x = 0 ; P.y = -1;
512  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
513  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
514  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);
515 
516  /* Point outside unit semicircle, sqrt(2)-1 units from arc end point*/
517  P.x = 1 ; P.y = 1;
519  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
520  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
521  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
522 
523  /* Point on unit semicircle midpoint */
524  P.x = 0 ; P.y = 1;
526  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
527  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
528  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
529 
530  /* Point on unit semicircle endpoint */
531  P.x = 1 ; P.y = 0;
533  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
534  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
535  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
536 
537  /* Point on semicircle center */
538  P.x = 0 ; P.y = 0;
540  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
541  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
542  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
543 
544  /* Point inside closed circle */
545  P.x = 0 ; P.y = 0.5;
546  A2.x = 1; A2.y = 0;
547  A3 = A1;
549  rv = lw_dist2d_pt_arc(&P, &A1, &A2, &A3, &dl);
550  //printf("distance %g\n", dl.distance);
551  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
552  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
553 }
554 
555 static void
557 {
558  /* int lw_dist2d_seg_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *B1, const POINT2D *B2, const POINT2D *B3, DISTPTS *dl) */
559 
560  DISTPTS dl;
561  POINT2D A1, A2, B1, B2, B3;
562  int rv;
563 
564  /* Unit semicircle */
565  B1.x = -1; B1.y = 0;
566  B2.x = 0 ; B2.y = 1;
567  B3.x = 1 ; B3.y = 0;
568 
569  /* Edge above the unit semicircle */
571  A1.x = -2; A1.y = 2;
572  A2.x = 2 ; A2.y = 2;
573  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
574  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
575  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
576 
577  /* Edge to the right of the unit semicircle */
579  A1.x = 2; A1.y = -2;
580  A2.x = 2; A2.y = 2;
581  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
582  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
583  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
584 
585  /* Edge to the left of the unit semicircle */
587  A1.x = -2; A1.y = -2;
588  A2.x = -2; A2.y = 2;
589  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
590  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
591  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
592 
593  /* Edge within the unit semicircle */
595  A1.x = 0; A1.y = 0;
596  A2.x = 0; A2.y = 0.5;
597  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
598  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
599  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
600 
601  /* Edge grazing the unit semicircle */
603  A1.x = -2; A1.y = 1;
604  A2.x = 2; A2.y = 1;
605  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
606  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
607  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0., 0.000001);
608 
609  /* Line grazing the unit semicircle, but edge not */
611  A1.x = 1; A1.y = 1;
612  A2.x = 2; A2.y = 1;
613  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
614  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
615  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
616 
617  /* Edge intersecting the unit semicircle */
619  A1.x = 0; A1.y = 0;
620  A2.x = 2; A2.y = 2;
621  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
622  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
623  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
624 
625  /* Line intersecting the unit semicircle, but edge not */
627  A1.x = -1; A1.y = 1;
628  A2.x = -2; A2.y = 2;
629  rv = lw_dist2d_seg_arc(&A1, &A2, &B1, &B2, &B3, &dl);
630  //printf("distance %g\n", dl.distance);
631  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
632  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1, 0.000001);
633 }
634 
635 static void
637 {
638  /* lw_dist2d_arc_arc(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3,
639  const POINT2D *B1, const POINT2D *B2, const POINT2D *B3,
640  DISTPTS *dl) */
641  DISTPTS dl;
642  POINT2D A1, A2, A3, B1, B2, B3;
643  int rv;
644 
645  /* Unit semicircle at 0,0 */
646  B1.x = -1; B1.y = 0;
647  B2.x = 0 ; B2.y = 1;
648  B3.x = 1 ; B3.y = 0;
649 
650  /* Arc above the unit semicircle */
652  A1.x = -1; A1.y = 3;
653  A2.x = 0 ; A2.y = 2;
654  A3.x = 1 ; A3.y = 3;
655  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
656  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
657  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
658 
659  /* Arc grazes the unit semicircle */
661  A1.x = -1; A1.y = 2;
662  A2.x = 0 ; A2.y = 1;
663  A3.x = 1 ; A3.y = 2;
664  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
665  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
666  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
667 
668  /* Circles intersect, but arcs do not */
670  A1.x = -1; A1.y = 1;
671  A2.x = 0; A2.y = 2;
672  A3.x = 1; A3.y = 1;
673  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
674  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
675  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2)-1, 0.000001);
676 
677  /* Circles and arcs intersect */
679  A1.x = -1; A1.y = 1;
680  A2.x = 0; A2.y = 0;
681  A3.x = 1; A3.y = 1;
682  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
683  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
684  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0, 0.000001);
685 
686  /* Concentric: and fully parallel */
688  A1.x = -2.0; A1.y = 0.0;
689  A2.x = 0.0; A2.y = 2.0;
690  A3.x = 2.0; A3.y = 0.0;
691  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
692  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
693  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
694 
695  /* Concentric: with A fully included in B's range */
697  A1.x = -0.5 / sqrt(2.0); A1.y = 0.5 / sqrt(2.0);
698  A2.x = 0.0; A2.y = 0.5;
699  A3.x = 0.5 / sqrt(2.0); A3.y = 0.5 / sqrt(2.0);
700  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
701  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
702  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
703 
704  /* Concentric: with A partially included in B's range */
706  A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
707  A2.x = -0.5; A2.y = 0.0;
708  A3.x = -0.5 / sqrt(2.0); A3.y = 0.5 / sqrt(2.0);
709  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
710  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
711  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
712 
713  /* Concentric: with A and B without parallel segments */
715  A1.x = -0.5 / sqrt(2.0); A1.y = -0.5 / sqrt(2.0);
716  A2.x = 0.0; A2.y = -0.5;
717  A3.x = 0.5 / sqrt(2.0); A3.y = -0.5 / sqrt(2.0);
718  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
719  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
720  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.736813, 0.000001);
721 
722  /* Concentric: Arcs are the same */
724  A1.x = -1.0; A1.y = 0.0;
725  A2.x = 0.0; A2.y = 1.0;
726  A3.x = 1.0; A3.y = 0.0;
727  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
728  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
729  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.0, 0.000001);
730 
731  /* Check different orientations */
732  B1.x = -10.0; B1.y = 0.0;
733  B2.x = 0.0 ; B2.y = 10.0;
734  B3.x = 10.0 ; B3.y = 0.0;
735 
737  A1.x = -22.0; A1.y = 0.0;
738  A2.x = -17.0; A2.y = -5.0;
739  A3.x = -12.0; A3.y = 0.0;
740  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
741  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
742  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 2.0, 0.000001);
743 
745  A1.x = -19.0; A1.y = 0.0;
746  A2.x = -14.0; A2.y = -5.0;
747  A3.x = - 9.0; A3.y = 0.0;
748  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
749  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
750  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
751 
753  A1.x = -9.0; A1.y = 0.0;
754  A2.x = -4.0; A2.y = -5.0;
755  A3.x = 1.0; A3.y = 0.0;
756  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
757  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
758  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
759 
761  A1.x = -1.0; A1.y = 0.0;
762  A2.x = 4.0; A2.y = -5.0;
763  A3.x = 9.0; A3.y = 0.0;
764  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
765  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
766  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
767 
769  A1.x = 1.0; A1.y = 0.0;
770  A2.x = 6.0; A2.y = -5.0;
771  A3.x = 11.0; A3.y = 0.0;
772  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
773  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
774  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
775 
777  A1.x = 11.0; A1.y = 0.0;
778  A2.x = 16.0; A2.y = -5.0;
779  A3.x = 21.0; A3.y = 0.0;
780  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
781  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
782  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
783 
784 
786  A1.x = -15.0; A1.y = -6.0;
787  A2.x = -10.0; A2.y = -1.0;
788  A3.x = - 5.0; A3.y = -6.0;
789  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
790  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
791  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1.0, 0.000001);
792 
794  A1.x = -5.0; A1.y = 0.0;
795  A2.x = 0.0; A2.y = 5.0;
796  A3.x = 5.0; A3.y = 0.0;
797  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
798  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
799  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 5.0, 0.000001);
800 
802  A1.x = -5.0; A1.y = 0.0;
803  A2.x = 0.0; A2.y = -5.0;
804  A3.x = 5.0; A3.y = 0.0;
805  rv = lw_dist2d_arc_arc(&A1, &A2, &A3, &B1, &B2, &B3, &dl);
806  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
807  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 5.0, 0.000001);
808 
809 
810 }
811 
812 static void
814 {
815 /* double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3) */
816 
817  POINT2D A1, A2, A3;
818  double d;
819 
820  /* Unit semicircle at 0,0 */
821  A1.x = -1; A1.y = 0;
822  A2.x = 0 ; A2.y = 1;
823  A3.x = 1 ; A3.y = 0;
824 
825  /* Arc above the unit semicircle */
826  d = lw_arc_length(&A1, &A2, &A3);
827  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
828  d = lw_arc_length(&A3, &A2, &A1);
829  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
830 
831  /* Unit semicircle at 0,0 */
832  A1.x = 0; A1.y = 1;
833  A2.x = 1; A2.y = 0;
834  A3.x = 0; A3.y = -1;
835 
836  /* Arc to right of the unit semicircle */
837  d = lw_arc_length(&A1, &A2, &A3);
838  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
839  d = lw_arc_length(&A3, &A2, &A1);
840  CU_ASSERT_DOUBLE_EQUAL(d, M_PI, 0.000001);
841 
842  /* Unit 3/4 circle at 0,0 */
843  A1.x = -1; A1.y = 0;
844  A2.x = 1; A2.y = 0;
845  A3.x = 0; A3.y = -1;
846 
847  /* Arc to right of the unit semicircle */
848  d = lw_arc_length(&A1, &A2, &A3);
849  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
850  d = lw_arc_length(&A3, &A2, &A1);
851  CU_ASSERT_DOUBLE_EQUAL(d, 3*M_PI_2, 0.000001);
852 }
853 
854 static void
856 {
857  /* lw_dist2d_pt_ptarrayarc(const POINT2D *p, const POINTARRAY *pa, DISTPTS *dl) */
858  DISTPTS dl;
859  int rv;
860  LWLINE *lwline;
861  POINT2D P;
862 
863  /* Unit semi-circle above X axis */
864  lwline = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
865 
866  /* Point at origin */
867  P.x = P.y = 0;
869  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
870  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
871  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
872 
873  /* Point above arc on Y axis */
874  P.y = 2;
876  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
877  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
878  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
879 
880  /* Point 45 degrees off arc, 2 radii from center */
881  P.y = P.x = 2 * cos(M_PI_4);
883  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
884  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
885  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
886 
887  /* Four unit semi-circles surrounding the 2x2 box around origin */
888  lwline_free(lwline);
889  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)"));
890 
891  /* Point at origin */
892  P.x = P.y = 0;
894  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
895  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
896  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0), 0.000001);
897 
898  /* Point on box edge */
899  P.x = -1; P.y = 0;
901  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
902  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
903  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
904 
905  /* Point within a semicircle lobe */
906  P.x = -1.5; P.y = 0;
908  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
909  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
910  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
911 
912  /* Point outside a semicircle lobe */
913  P.x = -2.5; P.y = 0;
915  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
916  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
917  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
918 
919  /* Point outside a semicircle lobe */
920  P.y = -2.5; P.x = 0;
922  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
923  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
924  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
925 
926  /* Point outside a semicircle lobe */
927  P.y = 2; P.x = 1;
929  rv = lw_dist2d_pt_ptarrayarc(&P, lwline->points, &dl);
930  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
931  CU_ASSERT_DOUBLE_EQUAL(dl.distance, sqrt(2.0)-1.0, 0.000001);
932 
933  /* Clean up */
934  lwline_free(lwline);
935 }
936 
937 static void
939 {
940  /* int lw_dist2d_ptarray_ptarrayarc(const POINTARRAY *pa, const POINTARRAY *pb, DISTPTS *dl) */
941  DISTPTS dl;
942  int rv;
943  LWLINE *lwline1;
944  LWLINE *lwline2;
945 
946  /* Unit semi-circle above X axis */
947  lwline1 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-1 0, 0 1, 1 0)"));
948 
949  /* Line above top of semi-circle */
950  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 2, -1 2, 1 2, 2 2)"));
952  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
953  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
954  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
955 
956  /* Reversed arguments, should fail */
959  rv = lw_dist2d_ptarray_ptarrayarc(lwline1->points, lwline2->points, &dl);
960  //printf("%s\n", cu_error_msg);
961  CU_ASSERT_EQUAL( rv, LW_FAILURE );
962  CU_ASSERT_STRING_EQUAL("lw_dist2d_ptarray_ptarrayarc called with non-arc input", cu_error_msg);
963 
964  lwline_free(lwline2);
965 
966  /* Line along side of semi-circle */
967  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2 -3, -2 -2, -2 2, -2 3)"));
969  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
970  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
971  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 1, 0.000001);
972 
973  /* Four unit semi-circles surrounding the 2x2 box around origin */
974  lwline_free(lwline1);
975  lwline_free(lwline2);
976  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)"));
977  lwline2 = lwgeom_as_lwline(lwgeom_from_text("LINESTRING(-2.5 -3, -2.5 -2, -2.5 2, -2.5 3)"));
978  rv = lw_dist2d_ptarray_ptarrayarc(lwline2->points, lwline1->points, &dl);
979  CU_ASSERT_EQUAL( rv, LW_SUCCESS );
980  CU_ASSERT_DOUBLE_EQUAL(dl.distance, 0.5, 0.000001);
981 
982  lwline_free(lwline2);
983  lwline_free(lwline1);
984 }
985 
986 static void
988 {
989  LWGEOM *g1, *g2;
990  double m, dist;
991 
992  /* Invalid input, lack of dimensions */
993 
994  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
995  g2 = lwgeom_from_wkt("LINESTRING (0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
996  m = lwgeom_tcpa(g1, g2, NULL);
997  lwgeom_free(g1);
998  lwgeom_free(g2);
999  ASSERT_DOUBLE_EQUAL(m, -1.0);
1000  CU_ASSERT_STRING_EQUAL(
1001  "Both input geometries must have a measure dimension",
1002  cu_error_msg
1003  );
1004 
1005  /* Invalid input, not linestrings */
1006 
1007  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1008  g2 = lwgeom_from_wkt("POINT M (0 0 2)", LW_PARSER_CHECK_NONE);
1009  m = lwgeom_tcpa(g1, g2, NULL);
1010  lwgeom_free(g1);
1011  lwgeom_free(g2);
1012  ASSERT_DOUBLE_EQUAL(m, -1.0);
1013  CU_ASSERT_STRING_EQUAL(
1014  "Both input geometries must be linestrings",
1015  cu_error_msg
1016  );
1017 
1018  /* Invalid input, too short linestring */
1019 
1020  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1021  g2 = lwgeom_from_wkt("LINESTRING M(2 0 1)", LW_PARSER_CHECK_NONE);
1022  dist = -77;
1023  m = lwgeom_tcpa(g1, g2, &dist);
1024  lwgeom_free(g1);
1025  lwgeom_free(g2);
1026  ASSERT_DOUBLE_EQUAL(dist, -77.0); /* not touched */
1027  ASSERT_DOUBLE_EQUAL(m, -1.0);
1028  CU_ASSERT_STRING_EQUAL(
1029  "Both input lines must have at least 2 points", /* should be accepted ? */
1030  cu_error_msg
1031  );
1032 
1033  /* Invalid input, empty linestring */
1034 
1035  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1036  g2 = lwgeom_from_wkt("LINESTRING M EMPTY", LW_PARSER_CHECK_NONE);
1037  m = lwgeom_tcpa(g1, g2, NULL);
1038  lwgeom_free(g1);
1039  lwgeom_free(g2);
1040  ASSERT_DOUBLE_EQUAL(m, -1.0);
1041  CU_ASSERT_STRING_EQUAL(
1042  "Both input lines must have at least 2 points",
1043  cu_error_msg
1044  );
1045 
1046  /* Timeranges do not overlap */
1047 
1048  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1049  g2 = lwgeom_from_wkt("LINESTRING M(0 0 2, 0 0 5)", LW_PARSER_CHECK_NONE);
1050  m = lwgeom_tcpa(g1, g2, NULL);
1051  lwgeom_free(g1);
1052  lwgeom_free(g2);
1053  ASSERT_DOUBLE_EQUAL(m, -2.0); /* means timeranges do not overlap */
1054 
1055  /* One of the tracks is still, the other passes to that point */
1056 
1057  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1058  g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
1059  dist = -1;
1060  m = lwgeom_tcpa(g1, g2, &dist);
1061  ASSERT_DOUBLE_EQUAL(m, 1.0);
1062  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1063  CU_ASSERT( lwgeom_cpa_within(g1, g2, 0.0) == LW_TRUE );
1064  lwgeom_free(g1);
1065  lwgeom_free(g2);
1066 
1067  /* One of the tracks is still, the other passes at 10 meters from point */
1068 
1069  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 5)", LW_PARSER_CHECK_NONE);
1070  g2 = lwgeom_from_wkt("LINESTRING M(-10 10 1, 10 10 5)", LW_PARSER_CHECK_NONE);
1071  dist = -1;
1072  m = lwgeom_tcpa(g1, g2, &dist);
1073  ASSERT_DOUBLE_EQUAL(m, 3.0);
1074  ASSERT_DOUBLE_EQUAL(dist, 10.0);
1075  CU_ASSERT( lwgeom_cpa_within(g1, g2, 11.0) == LW_TRUE );
1076  CU_ASSERT( lwgeom_cpa_within(g1, g2, 10.0) == LW_TRUE );
1077  CU_ASSERT( lwgeom_cpa_within(g1, g2, 9.0) == LW_FALSE );
1078  lwgeom_free(g1);
1079  lwgeom_free(g2);
1080 
1081  /* Equal tracks, 2d */
1082 
1083  g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1084  g2 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1085  dist = -1;
1086  m = lwgeom_tcpa(g1, g2, &dist);
1087  lwgeom_free(g1);
1088  lwgeom_free(g2);
1089  ASSERT_DOUBLE_EQUAL(m, 10.0);
1090  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1091 
1092  /* Reversed tracks, 2d */
1093 
1094  g1 = lwgeom_from_wkt("LINESTRING M(0 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1095  g2 = lwgeom_from_wkt("LINESTRING M(10 0 10, 0 0 20)", LW_PARSER_CHECK_NONE);
1096  dist = -1;
1097  m = lwgeom_tcpa(g1, g2, &dist);
1098  lwgeom_free(g1);
1099  lwgeom_free(g2);
1100  ASSERT_DOUBLE_EQUAL(m, 15.0);
1101  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1102 
1103  /* Parallel tracks, same speed, 2d */
1104 
1105  g1 = lwgeom_from_wkt("LINESTRING M(2 0 10, 12 0 20)", LW_PARSER_CHECK_NONE);
1106  g2 = lwgeom_from_wkt("LINESTRING M(13 0 10, 23 0 20)", LW_PARSER_CHECK_NONE);
1107  dist = -1;
1108  m = lwgeom_tcpa(g1, g2, &dist);
1109  lwgeom_free(g1);
1110  lwgeom_free(g2);
1111  ASSERT_DOUBLE_EQUAL(m, 10.0);
1112  ASSERT_DOUBLE_EQUAL(dist, 11.0);
1113 
1114  /* Parallel tracks, different speed (g2 gets closer as time passes), 2d */
1115 
1116  g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1117  g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 9 0 20)", LW_PARSER_CHECK_NONE);
1118  dist = -1;
1119  m = lwgeom_tcpa(g1, g2, &dist);
1120  lwgeom_free(g1);
1121  lwgeom_free(g2);
1122  ASSERT_DOUBLE_EQUAL(m, 20.0);
1123  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1124 
1125  /* Parallel tracks, different speed (g2 left behind as time passes), 2d */
1126 
1127  g1 = lwgeom_from_wkt("LINESTRING M(4 0 10, 10 0 20)", LW_PARSER_CHECK_NONE);
1128  g2 = lwgeom_from_wkt("LINESTRING M(2 0 10, 6 0 20)", LW_PARSER_CHECK_NONE);
1129  dist = -1;
1130  m = lwgeom_tcpa(g1, g2, &dist);
1131  lwgeom_free(g1);
1132  lwgeom_free(g2);
1133  ASSERT_DOUBLE_EQUAL(m, 10.0);
1134  ASSERT_DOUBLE_EQUAL(dist, 2.0);
1135 
1136  /* Tracks, colliding, 2d */
1137 
1138  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
1139  g2 = lwgeom_from_wkt("LINESTRING M(5 -8 0, 5 8 10)", LW_PARSER_CHECK_NONE);
1140  dist = -1;
1141  m = lwgeom_tcpa(g1, g2, &dist);
1142  lwgeom_free(g1);
1143  lwgeom_free(g2);
1144  ASSERT_DOUBLE_EQUAL(m, 5.0);
1145  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1146 
1147  /* Tracks crossing, NOT colliding, 2d */
1148 
1149  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 10 0 10)", LW_PARSER_CHECK_NONE);
1150  g2 = lwgeom_from_wkt("LINESTRING M(8 -5 0, 8 5 10)", LW_PARSER_CHECK_NONE);
1151  dist = -1;
1152  m = lwgeom_tcpa(g1, g2, &dist);
1153  lwgeom_free(g1);
1154  lwgeom_free(g2);
1155  ASSERT_DOUBLE_EQUAL(m, 6.5);
1156  ASSERT_DOUBLE_EQUAL(rint(dist*100), 212.0);
1157 
1158  /* Same origin, different direction, 2d */
1159 
1160  g1 = lwgeom_from_wkt("LINESTRING M(0 0 1, 10 0 10)", LW_PARSER_CHECK_NONE);
1161  g2 = lwgeom_from_wkt("LINESTRING M(0 0 1, -100 0 10)", LW_PARSER_CHECK_NONE);
1162  dist = -1;
1163  m = lwgeom_tcpa(g1, g2, &dist);
1164  lwgeom_free(g1);
1165  lwgeom_free(g2);
1166  ASSERT_DOUBLE_EQUAL(m, 1.0);
1167  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1168 
1169  /* Same ending, different direction, 2d */
1170 
1171  g1 = lwgeom_from_wkt("LINESTRING M(10 0 1, 0 0 10)", LW_PARSER_CHECK_NONE);
1172  g2 = lwgeom_from_wkt("LINESTRING M(0 -100 1, 0 0 10)", LW_PARSER_CHECK_NONE);
1173  dist = -1;
1174  m = lwgeom_tcpa(g1, g2, &dist);
1175  lwgeom_free(g1);
1176  lwgeom_free(g2);
1177  ASSERT_DOUBLE_EQUAL(m, 10.0);
1178  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1179 
1180  /* Converging tracks, 3d */
1181 
1182  g1 = lwgeom_from_wkt("LINESTRING ZM(0 0 0 10, 10 0 0 20)", LW_PARSER_CHECK_NONE);
1183  g2 = lwgeom_from_wkt("LINESTRING ZM(0 0 8 10, 10 0 5 20)", LW_PARSER_CHECK_NONE);
1184  dist = -1;
1185  m = lwgeom_tcpa(g1, g2, &dist);
1186  lwgeom_free(g1);
1187  lwgeom_free(g2);
1188  ASSERT_DOUBLE_EQUAL(m, 20.0);
1189  ASSERT_DOUBLE_EQUAL(dist, 5.0);
1190 
1191  /* G1 stops at t=1 until t=4 to let G2 pass by, then continues */
1192  /* G2 passes at 1 meter from G1 t=3 */
1193 
1194  g1 = lwgeom_from_wkt("LINESTRING M(0 0 0, 0 1 1, 0 1 4, 0 10 13)", LW_PARSER_CHECK_NONE);
1195  g2 = lwgeom_from_wkt("LINESTRING M(-10 2 0, 0 2 3, 12 2 13)", LW_PARSER_CHECK_NONE);
1196  dist = -1;
1197  m = lwgeom_tcpa(g1, g2, &dist);
1198  lwgeom_free(g1);
1199  lwgeom_free(g2);
1200  ASSERT_DOUBLE_EQUAL(m, 3.0);
1201  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1202 
1203  /* Test for https://trac.osgeo.org/postgis/ticket/3136 */
1204 
1205  g1 = lwgeom_from_wkt("LINESTRING M (0 0 1432291464,2 0 1432291536) ", LW_PARSER_CHECK_NONE);
1206  g2 = lwgeom_from_wkt("LINESTRING M (0 0 1432291464,1 0 1432291466.25,2 0 1432291500)", LW_PARSER_CHECK_NONE);
1207  dist = -1;
1208  m = lwgeom_tcpa(g1, g2, &dist);
1209  lwgeom_free(g1);
1210  lwgeom_free(g2);
1211  ASSERT_DOUBLE_EQUAL(m, 1432291464.0);
1212  ASSERT_DOUBLE_EQUAL(dist, 0.0);
1213 
1214  /* Tracks share a single point in time */
1215 
1216  g1 = lwgeom_from_wkt("LINESTRINGM(0 0 0, 1 0 2)", LW_PARSER_CHECK_NONE);
1217  g2 = lwgeom_from_wkt("LINESTRINGM(0 0 2, 1 0 3)", LW_PARSER_CHECK_NONE);
1218  dist = -1;
1219  m = lwgeom_tcpa(g1, g2, &dist);
1220  lwgeom_free(g1);
1221  lwgeom_free(g2);
1222  ASSERT_DOUBLE_EQUAL(m, 2.0);
1223  ASSERT_DOUBLE_EQUAL(dist, 1.0);
1224 
1225 }
1226 
1227 static void
1229 {
1230  LWGEOM *g;
1231  int ret;
1232 
1233  g = lwgeom_from_wkt("POINT M(0 0 1)", LW_PARSER_CHECK_NONE);
1234  ret = lwgeom_is_trajectory(g);
1235  lwgeom_free(g);
1236  ASSERT_INT_EQUAL(ret, LW_FALSE); /* not a linestring */
1237 
1238  g = lwgeom_from_wkt("LINESTRING Z(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1239  ret = lwgeom_is_trajectory(g);
1240  lwgeom_free(g);
1241  ASSERT_INT_EQUAL(ret, LW_FALSE); /* no measure */
1242 
1243  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 1)", LW_PARSER_CHECK_NONE);
1244  ret = lwgeom_is_trajectory(g);
1245  lwgeom_free(g);
1246  ASSERT_INT_EQUAL(ret, LW_FALSE); /* same measure in two points */
1247 
1248  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 0)", LW_PARSER_CHECK_NONE);
1249  ret = lwgeom_is_trajectory(g);
1250  lwgeom_free(g);
1251  ASSERT_INT_EQUAL(ret, LW_FALSE); /* backward measure */
1252 
1253  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 1 0 3, 2 2 2)", LW_PARSER_CHECK_NONE);
1254  ret = lwgeom_is_trajectory(g);
1255  lwgeom_free(g);
1256  ASSERT_INT_EQUAL(ret, LW_FALSE); /* backward measure */
1257 
1258  g = lwgeom_from_wkt("LINESTRING M(0 0 1, 0 0 2)", LW_PARSER_CHECK_NONE);
1259  ret = lwgeom_is_trajectory(g);
1260  lwgeom_free(g);
1261  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (still) */
1262 
1263  g = lwgeom_from_wkt("LINESTRING M EMPTY", LW_PARSER_CHECK_NONE);
1264  ret = lwgeom_is_trajectory(g);
1265  lwgeom_free(g);
1266  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (empty) */
1267 
1268  g = lwgeom_from_wkt("LINESTRING M (0 0 1)", LW_PARSER_CHECK_NONE);
1269  ret = lwgeom_is_trajectory(g);
1270  lwgeom_free(g);
1271  ASSERT_INT_EQUAL(ret, LW_TRUE); /* ok (corner case) */
1272 }
1273 
1274 /*
1275 ** Used by test harness to register the tests in this file.
1276 */
1277 void measures_suite_setup(void);
1279 {
1280  CU_pSuite suite = CU_add_suite("measures", NULL, NULL);
1292  PG_ADD_TEST(suite, test_lwgeom_tcpa);
1294 }
static void test_lw_dist2d_seg_arc(void)
Definition: cu_measures.c:556
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:2034
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:655
Datum boundary(PG_FUNCTION_ARGS)
void lwfree(void *mem)
Definition: lwutil.c:214
double lwgeom_mindistance2d_tolerance(const LWGEOM *lw1, const LWGEOM *lw2, double tolerance)
Function handling min distance calculations and dwithin calculations.
Definition: measures.c:199
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:469
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1050
#define LW_SUCCESS
Definition: liblwgeom.h:65
void lwline_free(LWLINE *line)
Definition: lwline.c:63
static void do_test_mindistance2d_tolerance(char *in1, char *in2, double expected_res, int line)
Definition: cu_measures.c:34
static LWGEOM * lwgeom_from_text(const char *str)
Definition: cu_measures.c:24
LWGEOM * lwgeom_from_wkt(const char *wkt, const char check)
Definition: lwin_wkt.c:890
void lwgeom_request_interrupt(void)
Request interruption of any running code.
Definition: lwgeom_api.c:824
static void test_lw_dist2d_arc_arc(void)
Definition: cu_measures.c:636
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:1479
RECT_NODE * rect_tree_new(const POINTARRAY *pa)
Build a tree of nodes from a point array, one node per edge, and each with an associated measure rang...
Definition: lwtree.c:151
#define LW_FAILURE
Definition: liblwgeom.h:64
static void test_lw_dist2d_pt_arc(void)
Definition: cu_measures.c:483
static void test_lwgeom_segmentize2d(void)
Definition: cu_measures.c:404
static void test_lw_dist2d_pt_ptarrayarc(void)
Definition: cu_measures.c:855
int lwgeom_cpa_within(const LWGEOM *g1, const LWGEOM *g2, double maxdist)
Is the closest point of approach within a distance ?
#define LW_PARSER_CHECK_NONE
Definition: liblwgeom.h:1869
double x
Definition: liblwgeom.h:312
int rect_tree_contains_point(const RECT_NODE *node, const POINT2D *pt, int *on_boundary)
Definition: lwtree.c:34
static void test_lwgeom_tcpa(void)
Definition: cu_measures.c:987
double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)
Returns the length of a circular arc segment.
Definition: lwalgorithm.c:104
#define WKT_ISO
Definition: liblwgeom.h:1939
#define DIST_MIN
Definition: measures.h:17
#define LW_FALSE
Definition: liblwgeom.h:62
int lwgeom_parse_wkt(LWGEOM_PARSER_RESULT *parser_result, char *wktstr, int parse_flags)
Parse a WKT geometry string into an LWGEOM structure.
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:79
void cu_error_msg_reset()
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:61
Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM...
Definition: liblwgeom.h:1876
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:1169
int lw_dist2d_pt_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, DISTPTS *dl)
Definition: measures.c:1421
#define PG_ADD_TEST(suite, testfunc)
#define ASSERT_INT_EQUAL(o, e)
POINTARRAY ** rings
Definition: liblwgeom.h:441
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:1065
double y
Definition: liblwgeom.h:312
Datum distance(PG_FUNCTION_ARGS)
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:89
static void test_rect_tree_contains_point(void)
Definition: cu_measures.c:196
#define DIST2DTEST(str1, str2, res)
Definition: cu_measures.c:32
void measures_suite_setup(void)
Definition: cu_measures.c:1278
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.
Note that p1 and p2 are pointers into an independent POINTARRAY, do not free them.
Definition: lwtree.h:4
double distance
Definition: measures.h:24
static void test_lwgeom_locate_along(void)
Definition: cu_measures.c:456
LWGEOM * lwgeom_segmentize2d(LWGEOM *line, double dist)
Definition: lwgeom.c:668
int rect_tree_intersects_tree(const RECT_NODE *n1, const RECT_NODE *n2)
Definition: lwtree.c:55
static void test_mindistance2d_tolerance(void)
Definition: cu_measures.c:72
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:1274
static void test_lw_dist2d_ptarray_ptarrayarc(void)
Definition: cu_measures.c:938
#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:18
Structure used in distance-calculations.
Definition: measures.h:22
void lw_dist2d_distpts_init(DISTPTS *dl, int mode)
Definition: measures.c:53
char cu_error_msg[MAX_CUNIT_ERROR_LENGTH+1]
static void test_lwgeom_is_trajectory(void)
Definition: cu_measures.c:1228
static void test_rect_tree_intersects_tree(void)
Definition: cu_measures.c:346
static void test_lw_arc_length(void)
Definition: cu_measures.c:813
POINTARRAY * points
Definition: liblwgeom.h:406