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