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