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