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